st_join 不够精确，无法计算一个点是否在多边形内部

Question

我有一个数据框，其中包含不同蛇的纬度/经度坐标，如下所示：

df1 <- data.frame(snake = c("snake_1", "snake_2"),
                  lat = c(40.68491, 40.8),
                  long = c(-83.34485, -83.3)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

> df1
Simple feature collection with 2 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -83.34485 ymin: 40.68491 xmax: -83.3 ymax: 40.8
Geodetic CRS:  WGS 84
    snake                   geometry
1 snake_1 POINT (-83.34485 40.68491)
2 snake_2         POINT (-83.3 40.8)

我有第二个 sf 多边形数据框，如下所示：

df2 <- data.frame(lat = c(40.68477,
                          40.68477,
                          40.9923,
                          40.9923,
                          40.68477),
                  long = c(-83.51883,
                           -83.11233,
                           -83.11233,
                           -83.51883,
                           -83.51883)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON") %>% 
  mutate(region = "region_1")

> df2
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -83.51883 ymin: 40.68477 xmax: -83.11233 ymax: 40.9923
Geodetic CRS:  WGS 84
                        geometry   region
1 POLYGON ((-83.51883 40.6847... region_1

如您所见，snake_1和snake_2都应该在多边形的边界内，尽管snake_1接近最小纬度。您可以在 ggplot 中看到这两点：

ggplot() +
  geom_sf(data = df2, fill = "lightblue") +
  geom_sf(data = df1, color = "red", size = 3) +
  theme_minimal()

我的目标是查看每个点是否在多边形内部。当我使用以下代码执行此操作时，生成的数据帧指示 Snake_1 不在边界框内。怎么会这样？该点肯定在多边形的边界内。除了 st_join 之外，我还应该使用其他函数吗？

test <- st_join(df1, df2, join = st_within)
> test
Simple feature collection with 2 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -83.34485 ymin: 40.68491 xmax: -83.3 ymax: 40.8
Geodetic CRS:  WGS 84
    snake   region                   geometry
1 snake_1     <NA> POINT (-83.34485 40.68491)
2 snake_2 region_1         POINT (-83.3 40.8)

Answer 1

在使用未投影的坐标系处理仅由 4 个角定义的多边形时，可能会有点棘手。例如，人们可能会假设

df2

多边形底部边缘的纬度值是恒定的（即底部部分与平行线平行），但情况并非总是如此：

# a linestring defined by geographic coordinates, CRS WGS 84
# latitude values for both ends are the same
line_1 <- sf::st_as_sfc("LINESTRING (-83.51883 40.68477, -83.11233 40.68477)", crs = "WGS84")

# but when adding more vertices to the linestring, our latitudes form a curve
# that "bends" behind snake_1 latitude
(lats <- sf::st_coordinates(sf::st_segmentize(line_1, 2000))[,"Y"])
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
#>  [1] 40.68477 40.68479 40.68481 40.68483 40.68485 40.68486 40.68488 40.68489
#>  [9] 40.68490 40.68491 40.68492 40.68493 40.68494 40.68494 40.68495 40.68495
#> [17] 40.68495 40.68495 40.68495 40.68494 40.68494 40.68493 40.68492 40.68491
#> [25] 40.68490 40.68489 40.68488 40.68486 40.68485 40.68483 40.68481 40.68479
#> [33] 40.68477
plot(lats)
abline(h = 40.68491, col = "red")
text(1,40.68491, "snake_1 latitude", adj = c(0, -1))

一种选择是向多边形添加更多顶点，但应注意何时以及如何执行此操作。在这里，我们将暂时删除 CRS，以便纬度/经度值将被视为平面坐标。

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(dplyr)
library(ggplot2)

df1 <- data.frame(snake = c("snake_1", "snake_2"),
                  lat = c(40.68491, 40.8),
                  long = c(-83.34485, -83.3)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

df2 <- data.frame(lat = c(40.68477,40.68477,40.9923,40.9923,40.68477),
                  long = c(-83.51883,-83.11233,-83.11233,-83.51883,-83.51883)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON") %>% 
  mutate(region = "region_1")

# segmentize region polygon, add point per every .1 deg
df2_seg <- st_set_crs(df2, NA) %>% 
  st_segmentize(0.1) %>% 
  st_set_crs(4326)

# test with updates polygons
st_join(df1, df2_seg, join = st_within)
#> Simple feature collection with 2 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -83.34485 ymin: 40.68491 xmax: -83.3 ymax: 40.8
#> Geodetic CRS:  WGS 84
#>     snake   region                   geometry
#> 1 snake_1 region_1 POINT (-83.34485 40.68491)
#> 2 snake_2 region_1         POINT (-83.3 40.8)

---

这是另一种可视化这一点的尝试。让我们放大盒子，将其向北移动并使用圆锥投影。这不是最好的视觉效果，但我们仍然可以看到 4 角多边形与具有更多顶点的多边形有何不同。

bb1 <- 
  st_bbox(c(xmin = -85, xmax = -82, ymax = 66, ymin = 65.8), crs = "WGS84") %>% 
  st_as_sfc()
# 4 corners
st_as_text(bb1)
#> [1] "POLYGON ((-85 65.8, -82 65.8, -82 66, -85 66, -85 65.8))"

bb2 <- 
  st_bbox(c(xmin = -85, xmax = -82, ymax = 66, ymin = 65.8)) %>% 
  st_as_sfc() %>% 
  st_segmentize(0.1) %>% 
  st_set_crs("WGS84")
# point per every .1 deg
st_as_text(bb2)
#> [1] "POLYGON ((-85 65.8, -84.9 65.8, -84.8 65.8, -84.7 65.8, -84.6 65.8, -84.5 65.8, -84.4 65.8, -84.3 65.8, -84.2 65.8, -84.1 65.8, -84 65.8, -83.9 65.8, -83.8 65.8, -83.7 65.8, -83.6 65.8, -83.5 65.8, -83.4 65.8, -83.3 65.8, -83.2 65.8, -83.1 65.8, -83 65.8, -82.9 65.8, -82.8 65.8, -82.7 65.8, -82.6 65.8, -82.5 65.8, -82.4 65.8, -82.3 65.8, -82.2 65.8, -82.1 65.8, -82 65.8, -82 65.9, -82 66, -82.1 66, -82.2 66, -82.3 66, -82.4 66, -82.5 66, -82.6 66, -82.7 66, -82.8 66, -82.9 66, -83 66, -83.1 66, -83.2 66, -83.3 66, -83.4 66, -83.5 66, -83.6 66, -83.7 66, -83.8 66, -83.9 66, -84 66, -84.1 66, -84.2 66, -84.3 66, -84.4 66, -84.5 66, -84.6 66, -84.7 66, -84.8 66, -84.9 66, -85 66, -85 65.9, -85 65.8))"

ggplot() +
  geom_sf(data = bb2, aes(fill = "seg")) +
  geom_sf(data = bb1, aes(fill = "4-point"), alpha = .8) +
  scale_fill_brewer(type = "qual") +
  coord_sf(crs = "+proj=lcca +lon_0=-83 +lat_0=70") +
  theme_minimal()

^{创建于 2024-03-16，使用 reprex v2.1.0}