R中合并Shapefile后用coord_sf修正CRS并缩放至研究区的方法

第一步：正确转换坐标系 #

你的原始Shapefile是EPSG:4326，需要转换到带15846转换参数的EPSG:3355。操作如下：

library(sf)
library(ggplot2)

# 读取并合并三个Shapefile
shp1 <- st_read("你的文件1.shp")
shp2 <- st_read("你的文件2.shp")
shp3 <- st_read("你的文件3.shp")
combined_shp <- rbind(shp1, shp2, shp3)

# 定义带15846转换参数的目标CRS（若默认EPSG:3355未包含该参数，手动指定PROJ字符串）
target_crs <- st_crs("+proj=utm +zone=36 +ellps=GRS80 +towgs84=26,-152,-194,0,0,0,0 +units=m +no_defs")
# 转换坐标系
combined_shp_proj <- st_transform(combined_shp, target_crs)

第二步：缩放到目标研究区 #

有两种高效的方式实现区域缩放：

方式1：用coord_sf指定坐标范围 #

根据当前图层的CRS匹配设置范围：

• 若使用EPSG:4326（经纬度）：

ggplot() +
  geom_sf(data = combined_shp) +
  coord_sf(crs = st_crs(4326), 
           xlim = c(32, 34),  # 东经32-34度
           ylim = c(26, 28))  # 北纬26-28度

• 若使用转换后的EPSG:3355（米为单位）：

# 把经纬度范围转换为目标投影坐标
bbox_wgs84 <- st_bbox(c(xmin=32, xmax=34, ymin=26, ymax=28), crs=st_crs(4326))
bbox_proj <- st_transform(st_as_sfc(bbox_wgs84), target_crs) %>% st_bbox()

ggplot() +
  geom_sf(data = combined_shp_proj) +
  coord_sf(crs = target_crs,
           xlim = c(bbox_proj[1], bbox_proj[3]),
           ylim = c(bbox_proj[2], bbox_proj[4]))

方式2：先裁剪图层再绘图 #

先裁剪数据到目标区域，再绘图会更高效：

# 定义目标研究区边界（WGS84）
study_area <- st_bbox(c(xmin=32, xmax=34, ymin=26, ymax=28), crs=st_crs(4326)) %>% st_as_sfc()
# 裁剪图层
cropped_shp <- st_crop(combined_shp_proj, study_area)

# 自动适配裁剪后的区域绘图
ggplot() +
  geom_sf(data = cropped_shp) +
  coord_sf(crs = target_crs)

关键提示 #

• geom_sf不会自动统一CRS，所有图层需保持CRS一致，或通过coord_sf的crs参数统一转换后绘图。
• 切换CRS后仍显示世界地图，是因为未设置xlim/ylim，默认会展示全部数据范围。

内容的提问来源于stack exchange，提问作者Alice Hobbs