尽管具有相同的 crs，但加入 sf 对象仍然出错

Question

我想弄清楚为什么我无法将两个 sf 对象连接在一起。

我有一个坐标数据框：

# A tibble: 80,840 x 2
   ycoord xcoord
    <dbl>  <dbl>
 1   51.9   4.44
 2   52.1   4.31
 3   52.1   4.31
 4   52.1   4.31
 5   52.3   5.17
 6   51.9   4.45
 7   52.0   4.17
 8   52.0   5.64
 9   51.8   5.81
10   52.1   4.66

以及带有邻域代码（例如邮政编码）的邻域边界的 sf 对象。

Simple feature collection with 14175 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
Projected CRS: Amersfoort / RD New
First 10 features:
      bu_code                       geometry
1  BU00140000 MULTIPOLYGON (((233836.2 58...
2  BU00140001 MULTIPOLYGON (((233934 5819...
3  BU00140002 MULTIPOLYGON (((233998.8 58...
4  BU00140003 MULTIPOLYGON (((233187.2 58...
5  BU00140004 MULTIPOLYGON (((233379.8 58...
6  BU00140005 MULTIPOLYGON (((234098.5 58...
7  BU00140007 MULTIPOLYGON (((234261.6 58...
8  BU00140008 MULTIPOLYGON (((233337.5 58...
9  BU00140100 MULTIPOLYGON (((234798 5816...
10 BU00140101 MULTIPOLYGON (((234489 5816...

为了获得 bu_code 变量，我尝试了以下方法来创建一个 sf 对象，该对象具有与邻域代码 sf 对象相同的 crs。

locaties_sf <- 
  locaties_df %>% 
  st_as_sf(coords = c(x = "xcoord", y = "ycoord"),
           crs = st_crs(neigbourhood_boundries_sf))

现在好像有相同的crs数据：

Simple feature collection with 80840 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3.359644 ymin: 50.76168 xmax: 7.21876 ymax: 53.48695
Projected CRS: Amersfoort / RD New
# A tibble: 80,840 x 1
              geometry
 *         <POINT [m]>
 1 (4.442658 51.91203)
 2 (4.314241 52.07851)
 3 (4.313637 52.05612)
 4 (4.313137 52.05622)
 5 (5.174235 52.26729)
 6  (4.44545 51.87896)
 7 (4.167927 52.00326)
 8  (5.643643 52.0424)
 9 (5.814333 51.83075)
10 (4.657972 52.06908)
# … with 80,830 more rows

然后我尝试将两者结合起来：

st_join(locaties_sf,
        neigbourhood_boundries_sf,
        join = st_within)

但我只得到 NA。

如果我绘制两者，它们似乎仍然具有不同的坐标系。

这就是 boundriers_sf 街区的样子：

这就是 locaties_sf 的样子：

非常感谢大家的帮助。

P.S.: 抱歉，我的数据不可重现，但多边形对于 dput()

来说太大了

Answer 1

您的坐标不匹配 - 点似乎是纬度和经度（在球体上以度为单位），而多边形是投影 CRS（在平面上以米为单位）。

你所做的是你“声明”你的点是平面上的米 = 显然非常接近原点，因此离多边形不远（注意 Amersfoort / RD New 的预期坐标范围） .

如果我的猜测是正确的，您需要首先在未投影的 CRS 中声明点（WGS84 是一个很好的默认值），然后（并且只有这样）将它们转换为多边形的投影 CRS。

我建议像这样更改您的 st_as_sf 调用：

locaties_sf <- locaties_df %>% 
  st_as_sf(coords = c(x = "xcoord", y = "ycoord"),
           crs = 4326) %>% # WGS84 - a good default for unprojected coords
  st_tranform(28992)

尽管具有相同的 crs，但加入 sf 对象仍然出错

Joining sf objects go wrong despite having the same crs

r

sf

coordinate-transformation