st_crop 在以太平洋为中心的自然地球世界地图上

st_crop on pacific centred naturalearth world map

我正在尝试使用 naturalearth 地图数据作为 sf 对象来裁剪以太平洋为中心的世界地图。该地图是根据 和以下代码创建的:


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")

# define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world
# Centered in lon 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)


# modify world dataset to remove overlapping portions with world's polygons
world2 <- worldMap %>% st_difference(polygon)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# Transform
world3 <- world2 %>% st_transform(crs = target_crs)

ggplot(data = world3, aes(group = admin)) +
  geom_sf(fill = "grey")

我想裁剪到覆盖印度洋和太平洋的感兴趣区域(就像这张取自“地图”world2 投影的屏幕截图):

为了裁剪,我尝试了以下方法:


world4 <- st_crop(world3, c(xmin= 30, ymin = -40, xmax = 230, ymax = 40))


但是,上面的代码产生了不正确的剪裁,增量非常小 lat/long:

ggplot(data = world3, aes(group = admin)) +
  geom_sf(fill = "grey")

如果有人能帮我弄清楚如何使用纬度和经度裁剪来定义这张以太平洋为中心的地图的区域,我将非常感激。我知道还有其他方法可以获得以太平洋为中心的地图,但该项目的未来计划分析要求它是 rnaturalearth 地图。

您的代码存在的问题是您需要使用与 world3 对象相同的 CRS 来定义裁剪过滤器的 bbox。例如:

world4 <- st_crop(
  x = world3, 
  y = st_as_sfc(
    st_bbox(c(xmin= 30, ymin = -40, xmax = 230, ymax = 40), crs = 4326)
  ) %>% st_transform(target_crs)
)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
ggplot(data = world4, aes(group = admin)) +
  geom_sf(fill = "grey")