CRS一致的两个SF对象执行st_intersection报错的技术问询

问题背景 #

我有两个SF空间对象：

• veg_plots_2024

Reading layer `veg_plots_2024' from data source 
  `C:\XXXX\veg_plots_2024.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 940 features and 40 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 170071.9 ymin: 26088.89 xmax: 646894 ymax: 640827.2
Projected CRS: OSGB36 / British National Grid

• land_class

Reading layer `land_class' from data source 
  `C:\XXXX\land_class.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 45 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 54 ymin: 11 xmax: 656 ymax: 1218
Projected CRS: OSGB36 / British National Grid

尝试执行st_intersection(veg_plots_2024, land_class)做相交分析，以获取样方对应的土地类型，但反复报错：

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : 
  st_crs(x) == st_crs(y) is not TRUE

表面上两者CRS名称一致，但仅land_class出现问题，怀疑和边界框有关，尝试转换CRS反而引发更多错误，补充land_class的CRS WKT信息如下：

> cat(st_crs(land_class)$wkt)
PROJCRS["OSGB36 / British National Grid",
    BASEGEOGCRS["OSGB36",
        DATUM["Ordnance Survey of Great Britain 1936",
            ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6277]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",49,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-2,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.999601272,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",400,
            LENGTHUNIT["kilometre",1000],
            ID["EPSG",8806]],
        PARAMETER["False northing",-100,
            LENGTHUNIT["kilometre",1000],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["kilometre",1000,
                ID["EPSG",9036]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["kilometre",1000,
                ID["EPSG",9036]]]]

问题原因 #

对比两个对象的CRS细节就能发现问题：

• veg_plots_2024的坐标单位是米（边界框数值170071.9等符合英国国家网格的米级坐标规范）
• land_class的WKT显示坐标单位是千米，且假东、假北参数也以千米为单位，这导致两个对象的CRS仅名称相同，底层参数（单位）不匹配，所以st_crs(x) == st_crs(y)返回FALSE。
• land_class的边界框数值极小（xmin:54，xmax:656），正好是veg_plots_2024坐标除以1000后的范围，进一步验证了单位差异的问题。

解决方法 #

方法1：修正land_class的单位并转换坐标 #

先将land_class的坐标从千米转为米，同时修正CRS的单位信息：

# 方式一：直接缩放坐标后匹配CRS
land_class_m <- land_class %>%
  mutate(geometry = geometry * 1000) %>%
  st_set_crs(st_crs(veg_plots_2024))

# 方式二：先修正CRS的单位定义再转换
land_class_crs_fixed <- st_crs(land_class)
# 替换WKT中的千米单位为米
land_class_crs_fixed$wkt <- gsub("kilometre", "metre", land_class_crs_fixed$wkt)
land_class_crs_fixed$wkt <- gsub("LENGTHUNIT\\[\"metre\",1000", "LENGTHUNIT[\"metre\",1", land_class_crs_fixed$wkt)

land_class_m <- land_class %>%
  st_set_crs(land_class_crs_fixed) %>%
  st_transform(st_crs(veg_plots_2024))

方法2：直接指定正确的EPSG代码 #

英国国家网格的标准EPSG代码是27700（单位为米），直接用这个代码重新定义land_class的CRS：

# 先缩放坐标（千米转米），再设置正确的EPSG
land_class_m <- land_class %>%
  mutate(geometry = geometry * 1000) %>%
  st_set_crs(27700)

验证与执行 #

先确认CRS匹配：

st_crs(veg_plots_2024) == st_crs(land_class_m)

返回TRUE后，再运行相交分析：

intersection_result <- st_intersection(veg_plots_2024, land_class_m)

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