R实现地图绘制

# ---- Step 1: Install and Load Necessary Packages ----
# Ensure these packages are installed
# install.packages(c("sf", "ggplot2", "dplyr", "readr", "purrr", "ggspatial", "grid", "officer", "rvg"))

library(sf)
library(ggplot2)
library(dplyr)
library(readr)     # For read_csv
library(purrr)     # For map functional programming
library(ggspatial)
library(grid)
# library(showtext) # Removed showtext dependency
library(officer)   # Optional
library(rvg)       # Optional
# library(readxl)    # Optional

# ---- Step 2: Set Font ----
# Specify the desired font family directly in ggplot elements
# Ensure "Times New Roman" is available on your system
base_font_family <- "Times New Roman"

# ---- Step 3: Read geojson Files ----
cities_data <- st_read("China Map/中国_市.geojson")
districts_data <- st_read("China Map/中国_县.geojson")

# ---- Step 4: Define Target Cities and Districts ----
target_areas_cities_cn <- c("成都市", "绵阳市", "宜宾市", "德阳市", "南充市", "达州市", "泸州市", "乐山市", "自贡市", "眉山市", "遂宁市", "内江市", "广安市", "雅安市", "资阳市")
target_areas_districts_cn <- c("渝北区", "渝中区", "九龙坡区", "沙坪坝区", "南岸区", "江津区", "北碚区", "巴南区", "涪陵区", "璧山区", "万州区", "永川区", "长寿区", "合川区", "大渡口区", "铜梁区", "大足区", "綦江区", "荣창区", "开州区", "潼南区", "梁平区", "南川区", "云阳县", "忠县", "垫江县", "黔江区", "丰都县", "石柱土家族自治县", "彭水苗族土家族自治县", "武隆区", "城口县", "秀山土家族苗族自治县", "酉阳土家族苗族自治县", "江北区") # Ensure Jiangbei is here
target_areas_cities_en <- c("Chengdu", "Mianyang", "Yibin", "Deyang", "Nanchong", "Dazhou", "Luzhou", "Leshan", "Zigong", "Meishan", "Suining", "Neijiang", "Guang'an", "Ya'an", "Ziyang")
target_areas_districts_en <- c("Yubei", "Yuzhong", "Jiulongpo", "Shapingba", "Nan'an", "Jiangjin", "Beibei", "Banan", "Fuling", "Bishan", "Wanzhou", "Yongchuan", "Changshou", "Hechuan", "Dadukou", "Tongliang", "Dazu", "Qijiang", "Rongchang", "Kaizhou", "Tongnan", "Liangping", "Nanchuan", "Yunyang", "Zhongxian", "Dianjiang", "Qianjiang", "Fengdu", "Shizhu", "Pengshui", "Wulong", "Chengkou", "Xiushan", "Youyang", "Jiangbei") # Add Jiangbei EN
gb_code_for_jiangbei <- '156500105'

# ---- Step 5: Filter Geographic Data ----
sichuan_cities <- cities_data %>% filter(name %in% target_areas_cities_cn)
chongqing_districts <- districts_data %>% filter((name %in% target_areas_districts_cn & name != "江北区") | (name == "江北区" & gb == gb_code_for_jiangbei))
sichuan_chongqing <- rbind(sichuan_cities, chongqing_districts) %>% distinct(name, .keep_all = TRUE) %>% mutate(name = trimws(name)) # Trim names early

# ---- Step 6: Calculate Average Node Strength (2020-2024) ----
avg_node_strength_data <- NULL # Initialize
tryCatch({
  years <- 2020:2024
  base_path <- "Python/Output Data" # Base path for data folders
  source_col_name <- "Source"
  target_col_name <- "Target"
  prob_col_name   <- "Weight"
  
  # Function to read and process data for a single year
  read_edge_data <- function(year, base_path, src_col, trg_col, prob_col) {
    file_path <- file.path(base_path, as.character(year), "边数据.csv")
    if (!file.exists(file_path)) {
      warning("File not found for year ", year, ": ", file_path)
      return(NULL) # Return NULL if file doesn't exist
    }
    tryCatch({
      df <- read_csv(file_path, show_col_types = FALSE)
      # Check for required columns
      if (!all(c(src_col, trg_col, prob_col) %in% names(df))) {
        warning("Missing required columns in file for year ", year, ": ", file_path)
        return(NULL)
      }
      df %>%
        rename(source = !!sym(src_col), target = !!sym(trg_col), prob = !!sym(prob_col)) %>%
        mutate(source = trimws(source), target = trimws(target)) %>%
        # Ensure probability is numeric, coerce non-numeric to NA
        mutate(prob = as.numeric(prob)) %>%
        filter(!is.na(prob)) %>% # Remove rows where probability couldn't be read as numeric
        select(source, target, prob) %>%
        mutate(year = year) %>% # Add year column
        # Create a unique pair ID, sorting source/target to handle directionality
        mutate(pair_id = map2_chr(source, target, ~ paste(sort(c(.x, .y)), collapse = "-")))
    }, error = function(e) {
      warning("Error reading or processing file for year ", year, ": ", file_path, " - ", e$message)
      return(NULL)
    })
  }
  
  # Read all yearly data into a single data frame
  all_years_edges <- map_dfr(years, ~read_edge_data(.x, base_path, source_col_name, target_col_name, prob_col_name))
  
  if (nrow(all_years_edges) == 0) {
    stop("No valid edge data could be read from any year.")
  }
  
  message("Read data for ", length(unique(all_years_edges$year)), " years.")
  
  # Calculate the average probability for each unique edge pair across all years
  # Ensure source and target are consistent for each pair_id before averaging
  consistent_pairs <- all_years_edges %>%
    select(pair_id, source, target) %>%
    distinct(pair_id, .keep_all = TRUE) # Get one consistent source/target for each pair_id
  
  average_edges <- all_years_edges %>%
    group_by(pair_id) %>%
    # Calculate the mean probability for each pair, handling potential NA values if a pair doesn't exist in a year
    summarise(average_prob = mean(prob, na.rm = TRUE), .groups = 'drop') %>%
    # Replace NaN with 0 (can happen if a pair exists but only has NA probabilities, though filter step should prevent this)
    mutate(average_prob = ifelse(is.nan(average_prob), 0, average_prob)) %>%
    filter(average_prob > 0) %>% # Keep only edges with positive average strength
    # Join back the consistent source/target names
    left_join(consistent_pairs, by = "pair_id") %>%
    select(source, target, average_prob) # Keep relevant columns
  
  
  if (nrow(average_edges) == 0) {
    stop("No edges with positive average probability found.")
  }
  
  message("Calculated average probability for ", nrow(average_edges), " unique edges.")
  
  # Calculate Average Node Strength
  # Sum the average strength for each node (appearing as source or target)
  strength_source <- average_edges %>% select(node = source, strength = average_prob)
  strength_target <- average_edges %>% select(node = target, strength = average_prob)
  
  avg_node_strength_data <- bind_rows(strength_source, strength_target) %>%
    group_by(node) %>%
    summarise(avg_node_strength = sum(strength), .groups = 'drop') %>%
    rename(name = node) # Rename column to 'name' for joining
  
  message("Successfully calculated average node strength (2020-2024) for ", nrow(avg_node_strength_data), " nodes.")
  
}, error = function(e) {
  warning("Error calculating average node strength (2020-2024): ", e$message, "\nMap fill will be missing.")
  avg_node_strength_data <- NULL # Ensure it's NULL on error
})

# ---- Step 7: Merge Average Node Strength to Geographical Data ----
if (!is.null(avg_node_strength_data)) {
  sichuan_chongqing <- sichuan_chongqing %>%
    left_join(avg_node_strength_data, by = "name")
} else {
  # Add a dummy column if calculation failed, to avoid errors later
  sichuan_chongqing$avg_node_strength <- NA_real_ # Use numeric NA
  warning("Using NA for node strength due to calculation error.")
}


# ---- Step 8: Calculate Centroids and Prepare Labels ----
# Use st_point_on_surface for potentially better placement inside polygons
sichuan_chongqing_centroids <- suppressWarnings(st_point_on_surface(sichuan_chongqing))

# Manual Adjustments (keep if needed)
suppressWarnings({
  # Find indices safely
  idx_luzhou <- which(sichuan_chongqing_centroids$name == "泸州市")
  idx_yubei <- which(sichuan_chongqing_centroids$name == "渝北区")
  idx_meishan <- which(sichuan_chongqing_centroids$name == "眉山市")
  
  if(length(idx_luzhou) > 0) sichuan_chongqing_centroids[idx_luzhou, ]$geometry <- st_geometry(sichuan_chongqing_centroids[idx_luzhou, ]) + c(-0.23, 0.23)
  if(length(idx_yubei) > 0) sichuan_chongqing_centroids[idx_yubei, ]$geometry <- st_geometry(sichuan_chongqing_centroids[idx_yubei, ]) + c(0.03, -0.03)
  if(length(idx_meishan) > 0) sichuan_chongqing_centroids[idx_meishan, ]$geometry <- st_geometry(sichuan_chongqing_centroids[idx_meishan, ]) + c(0, 0.1)
})

# ---- Step 9: Define and Process Labels (English with Numbered Core) ----
# Combine CN and EN names for mapping (ensure Jiangbei is included)
name_map_cn <- c(target_areas_cities_cn, target_areas_districts_cn, "江北区")
name_map_en <- c(target_areas_cities_en, target_areas_districts_en, "Jiangbei") # Make sure Jiangbei is here
name_mapping <- setNames(name_map_en, name_map_cn) %>% .[!duplicated(names(.))] # Create named vector, remove duplicate CN names if any

# Define the core Chongqing districts for numbering
chongqing_core_cn <- c("江北区", "沙坪坝区", "九龙坡区", "渝中区", "大渡口区", "南岸区")
chongqing_core_labels <- c("1", "2", "3", "4", "5", "6")
core_mapping <- setNames(chongqing_core_labels, chongqing_core_cn)

# Apply mapping: use number if core, else use English name from mapping, else original name
sichuan_chongqing_centroids <- sichuan_chongqing_centroids %>%
  mutate(label = case_when(
    name %in% chongqing_core_cn ~ recode(name, !!!core_mapping),
    name %in% names(name_mapping) ~ recode(name, !!!name_mapping),
    TRUE ~ name # Fallback to original name if no mapping found
  ))

# ---- Step 10: Calculate Shared Boundary ----
shared_boundary <- NULL
tryCatch({
  # Ensure input layers are valid before proceeding
  if(!all(st_is_valid(sichuan_cities))) sichuan_cities <- st_make_valid(sichuan_cities)
  if(!all(st_is_valid(chongqing_districts))) chongqing_districts <- st_make_valid(chongqing_districts)
  
  sichuan_cities_proj <- st_transform(sichuan_cities, crs = 3857) # Use a projected CRS for geometric operations
  chongqing_districts_proj <- st_transform(chongqing_districts, crs = 3857)
  
  sichuan_boundary_proj <- suppressMessages(suppressWarnings(st_union(sichuan_cities_proj) %>% st_boundary()))
  chongqing_boundary_proj <- suppressMessages(suppressWarnings(st_union(chongqing_districts_proj) %>% st_boundary()))
  
  # Perform intersection only if boundaries are valid and not empty
  if(!st_is_empty(sichuan_boundary_proj) && !st_is_empty(chongqing_boundary_proj) && st_is_valid(sichuan_boundary_proj) && st_is_valid(chongqing_boundary_proj)){
    shared_boundary_proj <- suppressMessages(suppressWarnings(st_intersection(sichuan_boundary_proj, chongqing_boundary_proj)))
    if (length(shared_boundary_proj) > 0 && !st_is_empty(shared_boundary_proj)) {
      shared_boundary <- st_transform(shared_boundary_proj, crs = 4326) # Transform back to WGS84
      message("Shared boundary calculated.")
    } else {
      message("No shared boundary found or intersection empty.")
    }
  } else {
    message("Could not calculate shared boundary due to invalid or empty input boundaries.")
  }
}, error = function(e){
  warning("Could not calculate shared boundary: ", e$message)
})

# ---- Step 11: Ensure Correct CRS for Plotting ----
sichuan_chongqing <- st_transform(sichuan_chongqing, crs = 4326)
sichuan_chongqing_centroids <- st_transform(sichuan_chongqing_centroids, crs = 4326)
# shared_boundary is already transformed if calculated successfully

# ---- Step 12: Plot the Map ----
map_plot <- ggplot() +
  # Base Map: Fill by Average Node Strength
  geom_sf(data = sichuan_chongqing,
          aes(fill = avg_node_strength), # Use calculated average node strength
          color = "gray50", # Boundary color for regions
          linewidth = 0.2) +
  
  # Shared Boundary (if calculated)
  {
    if (!is.null(shared_boundary) && inherits(shared_boundary, "sf") && nrow(shared_boundary) > 0 && !st_is_empty(shared_boundary))
      geom_sf(data = shared_boundary, color = "black", linewidth = 0.8)
  } +
  
  # Text Labels
  geom_sf_text(data = sichuan_chongqing_centroids,
               aes(label = label),
               size = 2.5, # Adjust size as needed
               color = "black",
               family = base_font_family, # Use Times New Roman
               check_overlap = FALSE) + # Or use ggrepel::geom_text_repel
  
  # North Arrow
  annotation_north_arrow(location = "tl", which_north = "true",
                         pad_x = unit(0.3, "cm"), pad_y = unit(0.9, "cm"),
                         height = unit(1.2, "cm"), width = unit(1.2, "cm"),
                         style = north_arrow_fancy_orienteering) +
  
  # Scale Bar
  annotation_scale(location = "bl", width_hint = 0.3, # Bottom-left
                   pad_x = unit(0.6, "cm"), pad_y = unit(0.6, "cm"),
                   text_family = base_font_family) + # Use Times New Roman
  
  # Color Scale for Node Strength
  scale_fill_gradient(name = "Average Node Strength\n(2020-2024)", # Legend Title
                      low = "#FDE7D7", high = "#F08B5D", # Orange gradient
                      na.value = "grey90", # Color for NA values
                      n.breaks = 5) + # Adjust number of breaks in legend
  
  # Labels and Theme
  labs(title = "Sichuan-Chongqing Area: Average Node Strength (2020-2024)") +
  theme_minimal(base_family = base_font_family) + # Use Times New Roman
  theme(
    plot.title = element_text(hjust = 0.5, size=14, face="bold"),
    panel.grid.major = element_line(color = "grey85", linetype = "dashed", linewidth = 0.3),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.text = element_text(color = "black", size = 9, family=base_font_family), # Font for axis text
    axis.ticks = element_line(color = "grey70"),
    legend.title = element_text(size = 9, face = "bold", family=base_font_family), # Font for legend title
    legend.text = element_text(size = 8, family=base_font_family), # Font for legend text
    legend.position = "right",
    plot.margin = unit(c(0.5, 2.5, 0.5, 0.5), "cm") # Adjust right margin for annotation
  )

# ---- Step 13: Add Custom Annotation for Numbered Labels ----
# Define annotation text (ensure names match numbered labels 1-6)
annotation_text <- "Numbered Districts:\n1: Jiangbei\n2: Shapingba\n3: Jiulongpo\n4: Yuzhong\n5: Dadukou\n6: Nan'an"

# Create grob with Times New Roman font
annotation_grob <- textGrob(
  label = annotation_text,
  hjust = 0, vjust = 0,
  gp = gpar(fontsize = 9, fontfamily = base_font_family, lineheight = 1.3) # Use Times New Roman
)

# Add annotation outside plot area
# Using annotation_custom with Inf/-Inf and clip="off"
final_plot <- map_plot +
  annotation_custom(
    grob = annotation_grob,
    xmin = Inf, xmax = Inf, ymin = -Inf, ymax = -Inf # Place outside data range
  ) +
  coord_sf(crs = 4326, clip = "off", expand = TRUE) # Allow drawing outside, expand=TRUE might help positioning annotation


# ---- Step 14: Print and Save ----
# Print to RStudio plotting window
print(final_plot)

# Save to PPT (Optional)
# Make sure 'final_plot' object is used
# try({
#   ppt_path <- "images/Sichuan_Chongqing_Avg_Node_Strength_Map_2020_2024.pptx"
#   ppt <- read_pptx()
#   ppt <- add_slide(ppt, layout = "Blank", master = "Office Theme")
#   # Ensure rvg is loaded if using dml
#   plot_dml <- rvg::dml(ggobj = final_plot)
#   ppt <- ph_with(ppt, value = plot_dml, location = ph_location_fullsize())
#   print(ppt, target = ppt_path)
#   message("Plot saved to: ", ppt_path)
# }, error = function(e){
#   warning("Could not save to PPT: ", e$message)
# })

# Save as PNG (Optional)
# ggsave("images/Sichuan_Chongqing_Avg_Node_Strength_Map_2020_2024.png", plot = final_plot, width = 11, height = 8, dpi = 300, bg = "white")