在R中使用多纬度边界而非边界框裁剪多边形Shapefile

我正在R中开展项目，需将8个含年度数据的多边形Shapefile按纬度划分为15个区域（仅关注纬度，不考虑经度，因研究仅涉及浅水区），以分析区域面积的年度变化。

尝试用15个边界框实现分割，但遇到两个问题：

• 仅找到单个边界框生成新Shapefile的代码，未找到多边界框的实现方法；
• 不确定边界框是否为合适方案，能否仅使用纬度边界？

尝试生成边界框时发现，尽管多边形投影与地图投影均正确，但边界框在ArcGIS中显示为直线，未贴合纬度曲线（可能因仅选取单个纬度点）。这种情况是否会导致裁剪掉所需区域或包含无关区域？

想了解是否有替代方法，仅通过纬度分界点在R中分割Shapefile为新多边形（例如，某区域仅包含933707.9m N至1038634.6m N之间的Shapefile区域），无需经度参数。这种方法能否解决纬度曲线问题？或是ArcGIS中的显示不影响边界框的实际功能？

附上从CSV文件读取经纬度，通过循环生成边界框的代码：

library(sp)
library(rgdal)
library(raster)
for(i in 1:nrow(KelpBounds)){
  
  XCoords <- c(KelpBounds[i,6],KelpBounds[i,8], KelpBounds[i,9], KelpBounds [i,7])
  XCoords
  
  YCoords <- c(KelpBounds[i,4],KelpBounds[i,5], KelpBounds[i,2], KelpBounds [i,3])
  YCoords
  
  xym <- cbind(XCoords, YCoords)
  xym
  
  
  p = Polygon(xym)
  ps = Polygons(list(p),1)
  sps = SpatialPolygons(list(ps))
  plot(sps)
  
  proj4string(sps) = CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0
                       +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
  
  proj4string(sps)
  
  data = data.frame(f=99.9)
  spdf = SpatialPolygonsDataFrame(sps,data)
  spdf
  summary(spdf)
  
 ## writeOGR(spdf, dsn='NewRegions', layer= paste0("Region_",i,"_Matrix") 
       ##   driver= "ESRI Shapefile") 
  
  dsn <- layer <- gsub(".csv","",i)
  writeOGR(spdf, dsn, layer, driver="ESRI Shapefile")
  
}

一、边界框显示与功能的说明 #

ArcGIS中边界框显示为直线是投影后的视觉效果——你用的AEAlbers等积投影属于平面投影，地理坐标系中的曲线纬线会被转换为直线，但裁剪逻辑是基于投影坐标值的，和视觉显示无关。只要边界框的Y坐标范围是目标纬度对应的投影值，裁剪结果就不会出现多裁或少裁的问题，实际功能不受视觉表现影响。

二、多边界框批量处理的优化方法 #

你的循环代码可以调整为批量生成边界框并直接裁剪原始数据，解决多区域处理的问题：

library(sp)
library(rgdal)
library(raster)

# 1. 读取并合并所有原始Shapefile（假设文件存放在同一文件夹）
raw_shp_paths <- list.files(path = "你的原始Shapefile文件夹路径", pattern = "\\.shp$", full.names = TRUE)
raw_spdf <- do.call(rbind, lapply(raw_shp_paths, readOGR))

# 2. 批量生成边界框并裁剪输出
output_dir <- "NewRegions"
if(!dir.exists(output_dir)) dir.create(output_dir)

for(i in 1:nrow(KelpBounds)){
  # 提取当前区域的Y范围（纬度对应的投影坐标）
  y_vals <- KelpBounds[i,c(2,3,4,5)]
  y_min <- min(y_vals)
  y_max <- max(y_vals)
  # 取原始数据的全经度范围（无需考虑经度限制）
  x_min <- bbox(raw_spdf)[1,1]
  x_max <- bbox(raw_spdf)[1,2]
  
  # 构造边界框多边形
  xym <- cbind(c(x_min, x_max, x_max, x_min), c(y_min, y_min, y_max, y_max))
  p <- Polygon(xym)
  ps <- Polygons(list(p), ID = as.character(i))
  sps <- SpatialPolygons(list(ps))
  proj4string(sps) <- proj4string(raw_spdf)
  spdf_box <- SpatialPolygonsDataFrame(sps, data = data.frame(ID = i))
  
  # 裁剪原始数据到当前边界框
  clipped_spdf <- raster::intersect(raw_spdf, spdf_box)
  
  # 输出裁剪后的Shapefile
  writeOGR(clipped_spdf, 
           dsn = output_dir, 
           layer = paste0("Region_", i), 
           driver = "ESRI Shapefile",
           overwrite_layer = TRUE)
}

三、仅用纬度边界分割的替代方案 #

如果不想手动构造边界框，推荐使用sf包（sp包已逐步被替代）直接按Y坐标范围裁剪，更简洁高效：

library(sf)

# 读取原始Shapefile
raw_sf <- st_read("你的原始Shapefile路径/文件名.shp")

# 定义15个纬度区间的分界点（示例：包含16个值，对应15个区间）
lat_bounds <- c(933707.9, 1038634.6, ...) # 补充完整所有分界点

# 循环生成各区域Shapefile
output_dir <- "NewRegions"
if(!dir.exists(output_dir)) dir.create(output_dir)

for(i in 1:(length(lat_bounds)-1)){
  y_min <- lat_bounds[i]
  y_max <- lat_bounds[i+1]
  # 直接按Y范围裁剪，自动覆盖全经度
  clipped_sf <- raw_sf %>%
    st_crop(xmin = st_bbox(raw_sf)[1], xmax = st_bbox(raw_sf)[2],
            ymin = y_min, ymax = y_max)
  # 输出文件
  st_write(clipped_sf, paste0(output_dir, "/Region_", i, ".shp"), delete_layer = TRUE)
}

这个方法完全满足“仅用纬度分界”的需求，无需手动构建边界框，sf包对投影的处理更稳定，代码可读性也更高。

内容的提问来源于stack exchange，提问作者Delaney Chabot