在R中触摸边界时测试几何是否包括其他

Question

我是地理数据的新手，所以也许这是一个愚蠢的问题。 我有两个多边形，我想做的是测试每个“病房”多边形包含在哪个“选区”多边形中。基本上，查看每个选区所在的选区。

但是，我遇到了来自 sf 的函数 st_within() 的问题，因为它仅在第一个几何体完全在第二个几何体中时才返回 TRUE。当较小的单位接触较大单位的边界时，这是一个问题。这是一个例子

我有一个称为“选区”的多面体和一个称为“病房”的多面体。你可以在这里找到它们：https://www.dropbox.com/sh/xu08wm79rym00zz/AADFDpyPe0EuDSDSY-6SvqpYa?dl=0

选区中有 2 个对象：

plot(constituencies$geometry)

然后，我尝试加入他们的行列，了解选区是否在选区

wards <- st_transform(wards, crs = st_crs(constituencies))
test98 <- st_join(wards, constituencies, join=st_within)

结果会在选区名称和选区所在的标签上吐出许多 NA。

    head(test98)
Simple feature collection with 6 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 531936.6 ymin: 180569.5 xmax: 533595.1 ymax: 182082.5
projected CRS:  OSGB 1936 / British National Grid
  WD98CD       WD98NM                             name label                       geometry
1 00AAFA   Aldersgate                             <NA>  <NA> MULTIPOLYGON (((532104.9 18...
2 00AAFB      Aldgate                             <NA>  <NA> MULTIPOLYGON (((533319.2 18...
3 00AAFC    Bassishaw Cities of London and Westminster   107 MULTIPOLYGON (((532587.3 18...
4 00AAFD Billingsgate Cities of London and Westminster   107 MULTIPOLYGON (((533167.9 18...
5 00AAFE  Bishopsgate                             <NA>  <NA> MULTIPOLYGON (((533410.7 18...
6 00AAFF Bread Street Cities of London and Westminster   107 MULTIPOLYGON (((532300.3 18...

当我绘制它们以查看哪些名称中有 NA（红色）时，您可以看到它是触及选区边界的那些：

plot(constituencies_short$geometry)
plot(test98$geometry[!is.na(test98$name),], col = "green", add = T)
plot(test98$geometry[is.na(test98$name),], col = "red", add = T)

我假设正在发生以下两种情况之一：1) 接触边界不能算作完全在几何范围内，或 2) 选区和选区的边界可能不完全匹配，因此选区并不完全在几何范围内。选区。

我的问题是：有没有一种方法可以完美匹配边界，然后测试哪些选区在哪些选区之内，或者有没有一种方法可以测试哪些选区不是全部，而是选区的大部分边界？

非常感谢！

Answer 1

您的边界不完全对齐。让我们看看当我们将第一个选区与两个选区相交时会发生什么：

> st_intersection(wards$geom[1], constituencies$geom)
Geometry set for 2 features 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 531936.6 ymin: 181262.6 xmax: 532309 ymax: 182012
projected CRS:  OSGB 1936 / British National Grid
POLYGON ((532104.9 182011.9, 532126.3 181948, 5...
MULTIPOLYGON (((532022.4 181893.5, 532022.4 181...

我们得到两个输出特征，因为选区在两个选区都有一些输出特征。多少钱？让我们看看：

> st_area(st_intersection(wards$geom[1], constituencies$geom))
Units: [m^2]
[1] 1.298829e+05 9.542473e-01

所以一个是 129882m^2，另一个是 0.95m^2。您的边界错位，因此“错误”选区中有 1 平方米的病房！

一种解决方法是根据病房面积的值或比例计算完整交叉集的面积和阈值。

您可以使用st_intersects 获取所有重叠的列表：

> st_intersects(wards, constituencies)

Sparse geometry binary predicate list of length 25, where the predicate was `intersects'
first 10 elements:
 1: 1, 2
 2: 1
 3: 1
 4: 1
 5: 1, 2
 6: 1

该列表中具有一个元素的任何元素都是与一个且仅一个选区重叠的选区，否则您需要计算该选区与这些选区重叠的面积比例，以计算出该选区主要结束的选区。

Answer 2

由于我们只处理细微的差异，一种方法是在加入之前在病房几何形状周围创建一个轻微的负/插入缓冲区。

准备

library(sf)
library(tmap)

wards_transformed <- st_transform(wards, st_crs(constituencies))

无缓冲：

wards_transformed %>%
  st_join(constituencies, st_within)

带负缓冲：

wards_transformed %>%
  st_buffer(-10) %>%
  st_join(constituencies, st_within)