R绘制中国地图

Ggpaddy

已于 2024-10-17 23:15:32 修改

阅读量542

点赞数 5

CC 4.0 BY-SA版权

文章标签： java 前端数据库

于 2024-10-17 22:08:22 首次发布

本文链接：https://blog.csdn.net/Ggpaddy/article/details/143029506

# 安装并加载必要的包
# install.packages(c("sf", "ggspatial", "ggplot2", "cowplot", "colorspace", "readxl", "dplyr"))
library(sf)
library(ggspatial)
library(ggplot2)
library(cowplot)
library(colorspace)
library(readxl)
library(dplyr)

# 读取Excel数据
excel_data <- read_excel("D:/R画地图/极端高温merra2.xlsx")

# 检查Excel数据的列名
print(colnames(excel_data))

# 读取Shapefile
shapefile_data <- st_read("D:/R画地图/中国地图最新-审图号GS(2019)1822号-shp格式/中国地图最新-审图号GS(2019)1822号-shp格式/中国地图最新-审图号GS(2019)1822号-shp格式/市（等积投影）.shp")

# 检查Shapefile数据的列名
print(colnames(shapefile_data))

# 合并数据，使用正确的列名
merged_data <- left_join(shapefile_data, excel_data, by = c("市代码" = "id"))

# 修改字段名称，确保字段名不超过10个字符
merged_data <- merged_data %>%
  rename(
    P_Cd = `省代码`,  # 省代码
    P = `省`,         # 省
    C_Cd = `市代码`, # 市代码
    C = `市`,       # 市
    T = `类型`       # 类型
  )

# 检查并创建保存Shapefile的目录
output_directory <- "D:/R画地图/分年Shapefile/"
if (!dir.exists(output_directory)) {
  dir.create(output_directory, recursive = TRUE)
}

# 检查字符编码
check_encoding <- function(data) {
  lapply(data, function(column) {
    if (is.character(column)) {
      bad_chars <- column[!iconv(column, "UTF-8", "UTF-8", sub = "byte") %in% column]
      return(bad_chars)
    }
  })
}

# 打印编码问题
bad_chars <- check_encoding(merged_data)
print(bad_chars)

# 按年份分组并保存每个年份的Shapefile
for (year in 1990:2023) {
  year_data <- merged_data %>% filter(Year == year)
  
  # 检查是否有数据
  if (nrow(year_data) > 0) {
    # 确保字段名不超过10个字符
    colnames(year_data) <- gsub("(.{10}).*", "\\1", colnames(year_data))  # 限制字段名为10个字符
    
    # 保存Shapefile，指定覆盖现有文件
    st_write(year_data, paste0(output_directory, "中国地图_", year, ".shp"), append = FALSE)
  }
}










"D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\九段线.shp"
"D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\国界.shp"
library(sf)
library(ggplot2)
library(dplyr)
library(RColorBrewer)
library(ggspatial)
library(grid)

# 定义投影坐标系（修改为 Albers 等面积圆锥投影）
crs_projected <- "+proj=aea +lat_1=25 +lat_2=47 +lon_0=110 +lat_0=0 +ellps=WGS84 +units=m +no_defs"


# 读取2020年市级Shapefile数据
year_2020_data <- st_read("D:/R画地图/分年Shapefile/中国地图_2020.shp")

# 创建Total Strength的分类变量
year_2020_data <- year_2020_data %>%
  mutate(class = cut(TotalStren,
                     breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40),
                     include.lowest = TRUE))

# 读取国界的Shapefile
border <- st_read("D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\国界.shp")

# 读取九段线的Shapefile
nine_dash_line <- st_read("D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\九段线.shp")

# **修复无效的几何图形**
year_2020_data <- st_make_valid(year_2020_data)
border <- st_make_valid(border)
nine_dash_line <- st_make_valid(nine_dash_line)

# 将主地图数据转换为投影坐标系（使用新的 Albers 投影）
year_2020_data_proj <- st_transform(year_2020_data, crs = crs_projected)
border_proj <- st_transform(border, crs = crs_projected)
nine_dash_line_proj <- st_transform(nine_dash_line, crs = crs_projected)

# 定义南海插图的地理范围（使用地理坐标系）
south_china_sea_bbox <- st_bbox(c(xmin = 105, xmax = 126, ymin = 3, ymax = 28), crs = st_crs(4326))

# 裁剪数据到南海插图的地理范围（使用地理坐标系）
year_2020_data_geo <- st_transform(year_2020_data, crs = 4326)
border_geo <- st_transform(border, crs = 4326)
nine_dash_line_geo <- st_transform(nine_dash_line, crs = 4326)

south_china_sea_bbox_sfc <- st_as_sfc(south_china_sea_bbox)

# **再次修复转换后的数据的无效几何**
year_2020_data_geo <- st_make_valid(year_2020_data_geo)
border_geo <- st_make_valid(border_geo)

# **执行空间操作前，设置 s2 选项**
sf::sf_use_s2(FALSE)

# 裁剪数据
year_2020_data_scs <- st_intersection(year_2020_data_geo, south_china_sea_bbox_sfc)
border_scs <- st_intersection(border_geo, south_china_sea_bbox_sfc)

# 绘制南海插图（使用地理坐标系）
nine_map <- ggplot() +
  geom_sf(data = year_2020_data_scs, aes(fill = class), color = NA) +
  geom_sf(data = border_scs, fill = NA, color = "#241E3E", size = 0.3) +
  geom_sf(data = nine_dash_line_geo, color = "#241E3E", size = 0.5) +
  coord_sf(
    xlim = c(105, 122),
    ylim = c(3, 28),
    crs = st_crs(4326),
    expand = FALSE
  ) +
  scale_fill_manual(values = rev(brewer.pal(8, "Spectral")), name = NULL, na.value = "grey95") +
  theme_void() +
  theme(
    legend.position = "none",
    panel.border = element_rect(fill = NA, color = "grey10", linetype = 1, size = 0.8),
    plot.margin = unit(c(0, 0, 0, 0), "mm")
  )

# 将南海插图转换为 grob 对象
nine_map_grob <- ggplotGrob(nine_map)

# 绘制主地图，并添加南海插图
china_map <- ggplot() +
  geom_sf(data = year_2020_data_proj, aes(fill = class), color = "#241E3E", size = 0.1) +
  geom_sf(data = border_proj, fill = NA, color = "#241E3E", size = 0.3) +
  geom_sf(data = nine_dash_line_proj, color = "#241E3E", size = 0.5) +
  coord_sf(crs = crs_projected, default_crs = st_crs(4326)) +
  scale_x_continuous(expand = c(0, 0), limits = c(72, 142), breaks = seq(70, 140, 10)) +
  scale_y_continuous(expand = c(0, 0), limits = c(17, 55.5), breaks = seq(10, 60, 10)) +
  # 调整图例为两列
  scale_fill_manual(values = rev(brewer.pal(8, "Spectral")), name = "Total Strength", na.value = "grey95") +
  guides(fill = guide_legend(ncol = 2)) +  # 设置图例为两列
  theme_bw() +
  theme(
    axis.text = element_text(family = "serif", color = "black"),
    axis.title = element_blank(),
    legend.position = c(0.2, 0.3),
    legend.justification = c(1, 1),
    legend.key.width = unit(0.9, "cm"),
    legend.key.height = unit(0.8, "cm"),
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 10),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  # 调整比例尺的位置（删除指北针）
  annotation_scale(
    location = "bl",        # 将比例尺移动到左下角
    pad_x = unit(1, "cm"),  # 水平边距
    pad_y = unit(1, "cm")   # 垂直边距
  ) +
  # 删除指北针（已删除该行代码）
  # annotation_north_arrow(...) +
  annotation_custom(nine_map_grob, xmin = 122, xmax = 134, ymin = 16, ymax = 30)  # 调整插图的位置和大小

# 绘制并保存最终的地图
print(china_map)
ggsave("Total_Strength_Map_2020.png", plot = china_map, width = 15, height = 9, dpi = 300)

# 安装并加载必要的包
# install.packages(c("sf", "ggspatial", "ggplot2", "cowplot", "colorspace", "readxl", "dplyr"))
library(sf)
library(ggspatial)
library(ggplot2)
library(cowplot)
library(colorspace)
library(readxl)
library(dplyr)

# 读取Excel数据
excel_data <- read_excel("D:/R画地图/极端高温merra2.xlsx")

# 检查Excel数据的列名
print(colnames(excel_data))

# 读取Shapefile
shapefile_data <- st_read("D:/R画地图/中国地图最新-审图号GS(2019)1822号-shp格式/中国地图最新-审图号GS(2019)1822号-shp格式/中国地图最新-审图号GS(2019)1822号-shp格式/市（等积投影）.shp")

# 检查Shapefile数据的列名
print(colnames(shapefile_data))

# 合并数据，使用正确的列名
merged_data <- left_join(shapefile_data, excel_data, by = c("市代码" = "id"))

# 修改字段名称，确保字段名不超过10个字符
merged_data <- merged_data %>%
  rename(
    P_Cd = `省代码`,  # 省代码
    P = `省`,         # 省
    C_Cd = `市代码`, # 市代码
    C = `市`,       # 市
    T = `类型`       # 类型
  )

# 检查并创建保存Shapefile的目录
output_directory <- "D:/R画地图/分年Shapefile/"
if (!dir.exists(output_directory)) {
  dir.create(output_directory, recursive = TRUE)
}

# 检查字符编码
check_encoding <- function(data) {
  lapply(data, function(column) {
    if (is.character(column)) {
      bad_chars <- column[!iconv(column, "UTF-8", "UTF-8", sub = "byte") %in% column]
      return(bad_chars)
    }
  })
}

# 打印编码问题
bad_chars <- check_encoding(merged_data)
print(bad_chars)

# 按年份分组并保存每个年份的Shapefile
for (year in 1990:2023) {
  year_data <- merged_data %>% filter(Year == year)
  
  # 检查是否有数据
  if (nrow(year_data) > 0) {
    # 确保字段名不超过10个字符
    colnames(year_data) <- gsub("(.{10}).*", "\\1", colnames(year_data))  # 限制字段名为10个字符
    
    # 保存Shapefile，指定覆盖现有文件
    st_write(year_data, paste0(output_directory, "中国地图_", year, ".shp"), append = FALSE)
  }
}










"D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\九段线.shp"
"D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\国界.shp"
library(sf)
library(ggplot2)
library(dplyr)
library(RColorBrewer)
library(ggspatial)
library(grid)
library(classInt)
# 定义投影坐标系（使用 Albers 等面积圆锥投影）
crs_projected <- "+proj=aea +lat_1=25 +lat_2=47 +lon_0=110 +lat_0=0 +ellps=WGS84 +units=m +no_defs"

# 设置年份范围
years <- 1990:2023

# 遍历每个年份
for (year in years) {
  
  # 构建年份对应的 Shapefile 路径（不修改路径，仅替换年份）
  shapefile_path <- paste0("D:/R画地图/分年Shapefile/中国地图_", year, ".shp")
  
  # 读取对应年份的市级 Shapefile 数据
  year_data <- st_read(shapefile_path)
  
  # 创建 Total Strength 的分类变量
  class_intervals <- classIntervals(year_data$TotalStren, n = 7, style = "jenks")
  year_data <- year_data %>%
    mutate(class = cut(TotalStren,
                       breaks = class_intervals$brks,
                       include.lowest = TRUE))
  
  # **修复无效的几何图形**
  year_data <- st_make_valid(year_data)
  
  # 将主地图数据转换为投影坐标系
  year_data_proj <- st_transform(year_data, crs = crs_projected)
  
  # 读取国界的 Shapefile（路径不变）
  border <- st_read("D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\国界.shp")
  border <- st_make_valid(border)
  border_proj <- st_transform(border, crs = crs_projected)
  
  # 读取九段线的 Shapefile（路径不变）
  nine_dash_line <- st_read("D:/R画地图\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\中国地图最新-审图号GS(2019)1822号-shp格式\\九段线.shp")
  nine_dash_line <- st_make_valid(nine_dash_line)
  nine_dash_line_proj <- st_transform(nine_dash_line, crs = crs_projected)
  
  # 定义南海插图的地理范围（使用地理坐标系）
  south_china_sea_bbox <- st_bbox(c(xmin = 105, xmax = 126, ymin = 3, ymax = 28), crs = st_crs(4326))
  
  # 裁剪数据到南海插图的地理范围（使用地理坐标系）
  year_data_geo <- st_transform(year_data, crs = 4326)
  year_data_geo <- st_make_valid(year_data_geo)
  
  border_geo <- st_transform(border, crs = 4326)
  border_geo <- st_make_valid(border_geo)
  
  nine_dash_line_geo <- st_transform(nine_dash_line, crs = 4326)
  
  south_china_sea_bbox_sfc <- st_as_sfc(south_china_sea_bbox)
  
  # **执行空间操作前，设置 s2 选项**
  sf::sf_use_s2(FALSE)
  
  # 裁剪数据
  year_data_scs <- st_intersection(year_data_geo, south_china_sea_bbox_sfc)
  border_scs <- st_intersection(border_geo, south_china_sea_bbox_sfc)
  
  # 绘制南海插图（使用地理坐标系）
  nine_map <- ggplot() +
    geom_sf(data = year_data_scs, aes(fill = class), color = NA) +
    geom_sf(data = border_scs, fill = NA, color = "#241E3E", size = 0.3) +
    geom_sf(data = nine_dash_line_geo, color = "#241E3E", size = 0.5) +
    coord_sf(
      xlim = c(105, 122),
      ylim = c(3, 28),
      crs = st_crs(4326),
      expand = FALSE
    ) +
    scale_fill_manual(values = rev(brewer.pal(8, "Spectral")), name = NULL, na.value = "grey95") +
    theme_void() +
    theme(
      legend.position = "none",
      panel.border = element_rect(fill = NA, color = "grey10", linetype = 1, size = 0.8),
      plot.margin = unit(c(0, 0, 0, 0), "mm")
    )
  
  # 将南海插图转换为 grob 对象
  nine_map_grob <- ggplotGrob(nine_map)
  
  # 绘制主地图，并添加南海插图
  china_map <- ggplot() +
    geom_sf(data = year_data_proj, aes(fill = class), color = "#241E3E", size = 0.1) +
    geom_sf(data = border_proj, fill = NA, color = "#241E3E", size = 0.3) +
    geom_sf(data = nine_dash_line_proj, color = "#241E3E", size = 0.5) +
    coord_sf(crs = crs_projected, default_crs = st_crs(4326)) +
    scale_x_continuous(expand = c(0, 0), limits = c(72, 142), breaks = seq(70, 140, 10)) +
    scale_y_continuous(expand = c(0, 0), limits = c(17, 55.5), breaks = seq(10, 60, 10)) +
    # 调整图例为两列
    scale_fill_manual(values = rev(brewer.pal(8, "Spectral")), name = "Total Strength", na.value = "grey95") +
    guides(fill = guide_legend(ncol = 2)) +  # 设置图例为两列
    theme_bw() +
    theme(
      axis.text = element_text(family = "serif", color = "black"),
      axis.title = element_blank(),
      legend.position = c(0.2, 0.3),
      legend.justification = c(1, 1),
      legend.key.width = unit(0.9, "cm"),
      legend.key.height = unit(0.8, "cm"),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 10),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank()
    ) +
    # 调整比例尺的位置（删除指北针）
    annotation_scale(
      location = "bl",        # 将比例尺移动到左下角
      pad_x = unit(1, "cm"),  # 水平边距
      pad_y = unit(1, "cm")   # 垂直边距
    ) +
    # 删除指北针（已删除该行代码）
    # annotation_north_arrow(...) +
    annotation_custom(nine_map_grob, xmin = 122, xmax = 134, ymin = 16, ymax = 30)  # 调整插图的位置和大小
  
  # 绘制并保存最终的地图
  print(china_map)
  
  # 保存地图，文件名包含年份
  output_filename <- paste0("D:/R画地图/极端高温总强度/Total_Strength_Map_", year, ".png")
  ggsave(output_filename, plot = china_map, width = 15, height = 9, dpi = 300)
}