使用 r 中的 st_join 从 GIS 文件获取人口普查块 ID

Question

我正在 R 中处理两个数据集：带有纬度和经度值的地址列表，以及表示韩国人口普查区块边界的 GIS 文件。

GIS 文件包含每个人口普查区块 (TOT_REG_CD) 的代码以及其他信息。该文件名为 census_id.shp，如下所示：

> 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...
2   20220630 29010110 29010110010101 MULTIPOLYGON (((981831.7 18...
3   20220630 29010110 29010110020003 MULTIPOLYGON (((982701 1844...
4   20220630 29010110 29010110020002 MULTIPOLYGON (((982185.9 18...
5   20220630 29010110 29010110020001 MULTIPOLYGON (((982448.6 18...
6   20220630 29010110 29010110020007 MULTIPOLYGON (((982521.4 18...
7   20220630 29010110 29010110020006 MULTIPOLYGON (((982254 1845...
8   20220630 29010110 29010110020005 MULTIPOLYGON (((982027.1 18...
9   20220630 29010110 29010110020004 MULTIPOLYGON (((982666.7 18...
10  20220630 29010110 29010110020009 MULTIPOLYGON (((981869.9 18...

另一个文件包含一长串地址及其各自的纬度和经度。它看起来像这样：

   > head(add)
           lon      lat
    1 126.9904 37.57180
    2 127.0153 37.57254
    3 126.9860 37.56995
    4 126.9670 37.56739
    5 126.9710 37.57226
    6 126.9729 37.57483

我的目标是根据纬度和经度将 GIS 文件中的 TOT_REG_CD 和相关信息附加到地址数据中。

这是我尝试过的代码：

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) 仅显示 NA 值的问题。 我从上面的代码得到的是：

> 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)
2       <NA>   <NA>       <NA> POINT (127.0153 37.57254)
3       <NA>   <NA>       <NA>  POINT (126.986 37.56995)
4       <NA>   <NA>       <NA>  POINT (126.967 37.56739)
5       <NA>   <NA>       <NA>  POINT (126.971 37.57226)
6       <NA>   <NA>       <NA> POINT (126.9729 37.57483)
7       <NA>   <NA>       <NA> POINT (126.9829 37.58806)
8       <NA>   <NA>       <NA> POINT (126.9692 37.60682)
9       <NA>   <NA>       <NA> POINT (127.0203 37.57404)
10      <NA>   <NA>       <NA> POINT (126.9988 37.57916)

有什么帮助吗？

提前致谢

Answer 1

crs

st_as_sf()

的参数不会转换，而是定义输入的坐标参考系，

add

的经/纬度列，但这些值很可能是地理坐标（WGS84），绝对不是投影的 KGD2002。

尝试先使用

WGS84

导入您的点，然后才转换为目标 crs:

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