太平洋中心NaturalEarth世界地图sf对象的st_crop裁剪问题

错误原因 #

你裁剪失败的核心是st_crop传入的坐标单位与待裁剪图层的坐标单位不匹配：你传入的xmin=30等参数是WGS84坐标系下的经纬度值（单位为度），但world3已经转换为中心经度133°的eqc投影，该投影的坐标单位为米，30米的范围自然会导致裁剪结果异常偏小。

解决方案 #

以下提供两种可直接运行的实现方案：

方案1：转换裁剪边界后裁剪投影图层

先将你定义的经纬度裁剪边界转换为目标投影下的坐标，再对投影后的图层执行裁剪，适合后续分析需要使用投影后数据的场景：

library(rnaturalearth)
library(sf)
library(ggplot2)

# 原有基础代码
worldMap <- ne_countries(scale = "medium", returnclass = "sf") %>%
  st_make_valid()
target_crs <- st_crs("+proj=eqc +x_0=0 +y_0=0 +lat_0=0 +lon_0=133")
offset <- 180 - 133
polygon <- st_polygon(x = list(rbind(
  c(-0.0001 - offset, 90),
  c(0 - offset, 90),
  c(0 - offset, -90),
  c(-0.0001 - offset, -90),
  c(-0.0001 - offset, 90)
))) %>%
  st_sfc() %>%
  st_set_crs(4326)
world2 <- worldMap %>% st_difference(polygon)
world3 <- world2 %>% st_transform(crs = target_crs)

# 新增裁剪逻辑
# 1. 生成WGS84坐标系下的裁剪边界
crop_bounds_lonlat <- st_bbox(c(xmin= 30, ymin = -40, xmax = 230, ymax = 40), crs = st_crs(4326)) %>%
  st_as_sfc()
# 2. 将边界转换为目标投影的坐标
crop_bounds_proj <- st_transform(crop_bounds_lonlat, crs = target_crs)
# 3. 执行裁剪
world4 <- st_crop(world3, crop_bounds_proj)

# 绘图验证
ggplot(data = world4, aes(group = admin)) +
  geom_sf(fill = "grey")

方案2：先裁剪原始经纬度图层再投影

先在WGS84坐标系下完成经纬度范围裁剪，再执行原有的子午线切割、投影逻辑，操作更直观不易出错：

library(rnaturalearth)
library(sf)
library(ggplot2)

worldMap <- ne_countries(scale = "medium", returnclass = "sf") %>%
  st_make_valid()
# 先在WGS84坐标系下按经纬度裁剪
crop_bounds_lonlat <- st_bbox(c(xmin= 30, ymin = -40, xmax = 230, ymax = 40), crs = st_crs(4326))
world_cropped_lonlat <- st_crop(worldMap, crop_bounds_lonlat)

# 执行原有切割、投影逻辑
target_crs <- st_crs("+proj=eqc +x_0=0 +y_0=0 +lat_0=0 +lon_0=133")
offset <- 180 - 133
polygon <- st_polygon(x = list(rbind(
  c(-0.0001 - offset, 90),
  c(0 - offset, 90),
  c(0 - offset, -90),
  c(-0.0001 - offset, -90),
  c(-0.0001 - offset, 90)
))) %>%
  st_sfc() %>%
  st_set_crs(4326)
world2 <- world_cropped_lonlat %>% st_difference(polygon)
world3 <- world2 %>% st_transform(crs = target_crs)

# 绘图验证
ggplot(data = world3, aes(group = admin)) +
  geom_sf(fill = "grey")

内容的提问来源于stack exchange，提问作者Gemma Galbraith