尽管有相同的 crs，但加入 sf 对象会出错答案

【问题标题】：Joining sf objects go wrong despite having the same crs尽管有相同的 crs，但加入 sf 对象会出错
【发布时间】：2021-11-26 11:29:27
【问题描述】：

我试图弄清楚为什么我不能让两个 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 对象具有相同crs 的sf 对象。

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。

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

这就是 neigbourhood boundriers_sf 的样子：

这就是 locaties_sf 的样子：

非常感谢所有帮助。

P.S.：很抱歉我的数据不可重现，但多面体对于dput()来说太大了

【问题讨论】：

您好@Tdebeus，您的“点”数据采用经度/纬度：它们是地理坐标（CRS WGS 84）。因此，您必须使用st_transform() 将它们转换为投影坐标系“Amersfoort / RD New”。

标签： r sf coordinate-transformation

【解决方案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)

【讨论】：