如何判断.shp文件是否含平面几何及非平面坐标平面化方法

问题复现 #

调用st_intersection()计算两个多边形交集的代码如下：

poly.3 <- st_intersection(poly.1, poly.2)

运行后返回警告信息：

Warning: attribute variables are assumed to be spatially constant throughout all geometries
although coordinates are longitude/latitude, st_intersection assumes that they are planar

已提前执行poly.2 <- st_transform(poly.2, st_crs(poly.1))完成两个图层的空间参考对齐。 #

问题1：判断坐标是否为平面坐标的方法 #

你当前通过st_crs(poly.1)返回的坐标系不是平面坐标，判断依据有3个：

• 看CRS类型：返回的WKT标识开头为GEOGCRS，代表地理坐标系（经纬度坐标）；平面投影坐标系的WKT开头会显示为PROJCRS。
• 看坐标单位：当前坐标轴的单位是ANGLEUNIT["degree",0.0174532925199433]，即角度单位“度”；平面坐标的坐标轴单位为长度单位（米/英尺等）。
• 看EPSG编码：当前坐标系EPSG编号为4326，对应WGS84地理坐标系，是全球通用的经纬度坐标，不属于平面投影坐标。

注：第一条关于属性变量的警告属于常规提示，只要你的矢量数据属性值在单个几何范围内是空间恒定的，该提示不会影响计算结果，不需要额外处理。

问题2：非平面坐标转平面坐标的操作方法 #

sf_use_s2(FALSE)只是强制关闭球面几何计算逻辑，本质是将经纬度坐标直接当作平面坐标计算，会带来明显的计算误差，不推荐使用。正确操作流程如下：

1. 选择适配研究区域的平面投影坐标系：
   • 大尺度（如全国范围）分析：可选择对应区域的等积投影，例如中国区域可选CGCS2000 Albers等积投影（EPSG:4527）
   • 小尺度（如单个省市/区县）分析：选择研究区对应的UTM分带投影、高斯-克吕格投影即可，优先选择单位为米的投影，保证几何计算精度
2. 执行坐标转换与交集计算，示例代码（以转EPSG:4527为例）：

   # 再次确认两个输入图层CRS一致
   stopifnot(st_crs(poly.1) == st_crs(poly.2))
   # 统一转换到选定的平面投影坐标系
   poly.1_proj <- st_transform(poly.1, crs = 4527)
   poly.2_proj <- st_transform(poly.2, crs = 4527)
   # 执行交集计算
   poly.3 <- st_intersection(poly.1_proj, poly.2_proj)
   

转换完成后可再次调用st_crs()校验，当WKT开头为PROJCRS、坐标轴单位为米时，即说明已转换为平面坐标，此时计算不会再出现经纬度假设为平面的警告，计算结果精度符合分析要求。 #

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