如何用R的Terra包制作带属性的等面积中心对齐网格并导出为Shapefile

步骤说明 #

• 转换地理数据到等面积投影坐标系，确保网格面积精准匹配400平方英里要求
• 基于指定经纬度中心点生成对齐的网格，再裁剪至威斯康星州边界内
• 为单元格添加属性后导出为Shapefile

完整代码 #

library(terra)

# 1. 加载威斯康星州边界数据
wi_shape <- vect('C:\\Users\\ruben\\Downloads\\Wisconsin_State_Boundary_24K\\Wisconsin_State_Boundary_24K.shp')

# 2. 切换到等面积投影（威斯康星州专属EPSG:3071，单位为米）
target_crs <- "EPSG:3071"
wi_proj <- project(wi_shape, target_crs)

# 3. 定义指定中心点（替换为你的目标经纬度），转换到目标投影
center_lonlat <- c(-89.4, 43.0) # 示例中心点：威斯康星州大致地理中心
center_proj <- project(vect(matrix(center_lonlat, ncol=2), crs="EPSG:4326"), target_crs)
center_coords <- crds(center_proj)

# 4. 计算单元格尺寸：400平方英里=20英里边长，转换为米
cell_size_m <- 20 * 1609.344

# 5. 计算网格偏移量，让指定点成为某一网格的中心
offset_x <- center_coords[1] - (cell_size_m / 2)
offset_y <- center_coords[2] + (cell_size_m / 2) # 投影坐标系中y值向上递增

# 6. 生成覆盖全州的网格（扩展边界避免边缘截断）
wi_ext <- ext(wi_proj)
wi_ext_expand <- wi_ext + cell_size_m * 2
wi_grid <- make_grid(wi_proj, cellsize=cell_size_m, offset=c(offset_x, offset_y), square=TRUE)

# 7. 裁剪网格到州边界内
wi_grid_clipped <- intersect(wi_grid, wi_proj)

# 8. 添加单元格属性（ID、中心点经纬度等）
wi_grid_clipped$cell_id <- 1:nrow(wi_grid_clipped)
grid_centers <- centroids(wi_grid_clipped)
wi_grid_clipped$center_lon <- project(grid_centers, "EPSG:4326")[,1]
wi_grid_clipped$center_lat <- project(grid_centers, "EPSG:4326")[,2]

# 9. 导出为Shapefile
writeVector(wi_grid_clipped, "wisconsin_400sqmi_grid.shp", overwrite=TRUE)

# 可选：可视化验证
plot(wi_proj)
plot(wi_grid_clipped, add=TRUE, border="blue", lwd=0.5)
plot(center_proj, add=TRUE, col="red", pch=19)

关键细节说明 #

• 投影选择：EPSG:3071是威斯康星州本地等面积投影，能严格保证网格面积准确性；若需使用英里单位投影，需对应调整坐标系与单元格尺寸计算。
• 网格对齐：通过offset参数精准调整网格起始位置，确保指定中心点落在某一网格的几何中心。
• 裁剪逻辑：intersect会保留所有与州边界相交的单元格；若需仅保留完全位于州内的单元格，可改用crop后结合is.inside筛选。
• 属性扩展：可根据需求添加更多属性（如单元格面积、边界坐标等）。

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