R语言极投影下如何沿指定纬度线裁剪shapefile文件

云上老周

2026-8-29

问题说明

需要对shapefile执行裁剪时，让裁剪边界平滑贴合指定纬度线，例如将研究区底部边界设置为沿北纬70°纬线延伸。st_crop()默认采用当前坐标系下的矩形范围裁剪，直接使用会导致投影后边界无法贴合纬线，尝试组合st_segmentize()与st_intersection()未得到预期结果，复现代码如下：

library(ggplot2)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
library(rnaturalearth)

sf_use_s2(use_s2 = FALSE)
#> Spherical geometry (s2) switched off

ocean <-
  ne_download(
    category = "physical",
    type = "ocean",
    returnclass = "sf",
    scale = "large"
  ) 
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/RtmpUOEg4b", layer: "ne_10m_ocean"
#> with 1 features
#> It has 3 fields

bbox <- c(
  xmin = -90,
  xmax = -40,
  ymin = 70,
  ymax = 80
)

ocean |>
  st_crop(bbox) |>
  st_transform(6056) |>
  ggplot() +
  geom_sf(size = 0.2)
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

默认裁剪方式投影后效果
本示例于2022-06-21由reprex包（v2.0.1）创建生成

解决方法

问题核心原因：直接传入数值bbox给st_crop()时，函数会生成仅含4个顶点的矩形，矩形边在投影坐标系下为直线，仅靠四个顶点无法还原纬线的投影曲线，自然无法贴合指定纬度线。

按以下步骤操作即可实现沿纬线平滑裁剪：

在WGS84地理坐标系（EPSG:4326）下手动构造对应经纬度范围的闭合裁剪多边形
用st_segmentize()给裁剪多边形的边加密节点，保证投影后曲线平滑
用加密后的多边形和原始数据做相交裁剪，再将结果转换到目标投影

修正后代码：

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

sf_use_s2(use_s2 = FALSE)

# 读取原始海洋数据
ocean <-
  ne_download(
    category = "physical",
    type = "ocean",
    returnclass = "sf",
    scale = "large"
  ) 

# 定义裁剪经纬度范围
xmin <- -90
xmax <- -40
ymin <- 70
ymax <- 80

# 构造沿经纬线的闭合裁剪面，指定坐标系为WGS84
crop_poly <- st_polygon(list(matrix(c(
  xmin, ymin,
  xmax, ymin,
  xmax, ymax,
  xmin, ymax,
  xmin, ymin
), ncol = 2, byrow = TRUE))) |> 
  st_sfc(crs = 4326) |> 
  # 按0.1度间隔插入节点，足够保证投影后边界平滑
  st_segmentize(dfMaxLength = 0.1)

# 先在地理坐标系下完成裁剪，再转换到目标投影EPSG:6056
ocean_cropped <- ocean |> 
  st_intersection(crop_poly) |> 
  st_transform(6056)

# 绘图查看结果
ggplot(ocean_cropped) +
  geom_sf(size = 0.2)

注意事项：

不要先将数据转换到投影坐标系再裁剪，否则裁剪边界为投影坐标系下的矩形，无法匹配纬线
st_segmentize()的dfMaxLength参数可根据需求调整，数值越小节点越密，边界越平滑，对应计算量也会略有上升

内容的提问来源于stack exchange，提问作者Philippe Massicotte