使用R中st_join关联GIS普查区块数据时返回NA值的问题求助

我在R中处理两个数据集：

• 包含经度（lon）、纬度（lat）的地址列表add
• 韩国普查区块边界GIS文件census_id.shp（含BASE_DATE、ADM_CD、TOT_REG_CD字段）

目标是根据经纬度将GIS文件中的字段匹配到地址数据中，使用sf包编写的代码如下：

library(sf)
add
head(add)

census_BL_boundary <- st_read("census_id.shp")
census_BL_boundary

add_sf <- st_as_sf(add, coords = c("lon", "lat"), crs = st_crs(census_BL_boundary))
joined_data <- st_join(add_sf, census_BL_boundary)

执行后joined_data中BASE_DATE、ADM_CD、TOT_REG_CD均为NA值，相关数据示例：

• 普查区块GIS数据：

> census_BL_boundary
Simple feature collection with 104292 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 746111 ymin: 1458603 xmax: 1387949 ymax: 2068444
Projected CRS: KGD2002 / Unified CS
First 10 features:
   BASE_DATE   ADM_CD     TOT_REG_CD                       geometry
1   20220630 29010110 29010110010001 MULTIPOLYGON (((982172 1846...
...

• 地址数据：

> head(add)
           lon      lat
    1 126.9904 37.57180
    2 127.0153 37.57254
...

• 关联结果：

> joined_data
Simple feature collection with 2025 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 126.5687 ymin: 36.98421 xmax: 127.7803 ymax: 38.09119
Projected CRS: KGD2002 / Unified CS
First 10 features:
   BASE_DATE ADM_CD TOT_REG_CD                  geometry
1       <NA>   <NA>       <NA>  POINT (126.9904 37.5718)
...

核心问题是坐标系不匹配：

• 地址数据的经纬度属于地理坐标系（WGS84，EPSG:4326），单位是度
• 普查区块GIS文件使用的是投影坐标系（KGD2002 / Unified CS，EPSG:5179），单位是米

之前的代码直接将经纬度点的CRS指定为投影坐标系，相当于把WGS84的经纬度值当作KGD2002的坐标使用，导致点的位置完全错误，不在任何普查区块多边形范围内，因此匹配结果全为NA。

需要先将地址数据转换为正确的地理坐标系，再转换为与GIS文件一致的投影坐标系，最后执行空间连接：

library(sf)

# 读取地址数据（假设add已加载）
head(add)

# 读取普查区块边界
census_BL_boundary <- st_read("census_id.shp")
census_BL_boundary

# 1. 将地址数据转为sf对象，指定正确的地理坐标系WGS84（EPSG:4326）
add_sf <- st_as_sf(add, coords = c("lon", "lat"), crs = 4326)

# 2. 将地址数据的坐标系转换为与普查边界一致的投影坐标系
add_sf_transformed <- st_transform(add_sf, crs = st_crs(census_BL_boundary))

# 3. 执行空间连接，使用st_within确保点在多边形内才匹配，left连接保留所有地址点
joined_data <- st_join(add_sf_transformed, census_BL_boundary, join = st_within)

# 查看匹配结果
head(joined_data)

1. 检查转换后的地址点坐标范围，应该和普查区块的边界范围（xmin:746111_{xmax:1387949，ymin:1458603}ymax:2068444）一致
2. 查看joined_data中的字段是否不再为NA，确认匹配成功

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