如何解决空间数据连接时出现的球面几何失败问题

Question

我有一个包含多个多边形的shapefile文件和一个带有坐标的数据框。我想将数据框中的每个坐标分配给shapefile中的一个多边形。因此，需要在数据框中添加一列，列名为多边形名称或ID。

这是数据的链接

library(sf)
library(readr)
shape <- read_sf("data/Provinces_v1_2017.shp")
data<- read_csv("data/data.csv")

但是当我尝试加入它们时，我总是会出现错误

pts = st_as_sf(data, coords = c("dec_lon", "dec_lat"), crs= 4326)

st_join(pts, shape)

我尝试了over()函数和其他技巧，比如st_make_valid()，但是总是会报错：

Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : Evaluation error: Found 30 features with invalid spherical geometry.

这是最近出现的问题（之前我的代码可以正常工作），但现在我无法使用sf包完成这个任务，总是遇到这个错误。我更新了库以查看是否有帮助，但无法解决。

非常感谢您的帮助。

- yuliaUU

Answer 1

您有两个选择：

1. 通过在脚本中使用sf::sf_use_s2(FALSE)关闭s2处理；理论上，行为应该恢复到1.0发布之前的状态。
2. 修复您的多边形对象的球面几何；这将取决于实际错误的性质。

我无法访问您的文件并进行确认，但这段代码过去曾经帮助过我：

yer_object$geometry <- yer_object$geometry %>%
  s2::s2_rebuild() %>%
  sf::st_as_sfc()

- Jindra Lacko

---

你是正确的! 我尝试了 sf::sf_use_s2(FALSE)，现在它可以工作了! - yuliaUU

非常荣幸能为您服务！ :) - Jindra Lacko

我也遇到了这个错误，让我很疯狂，因为在Python中使用相同的输入时操作是顺利运行的。@jiladata，您能否解释一下在1.0版本中发生了什么变化导致了意外的行为？谢谢。 - IvanP

嗨 @IvanP！ sf v1.0 的更改是从 GEOS 到 Google 的 s2，针对未投影坐标（即地理坐标，例如 EPSG 4326 中的经纬度）的后端引擎。GEOS 将投影坐标视为平面（即两个点位于无限长度的直线上），而 s2 更为“正确”（两个点位于周长为 40,075 公里的大圆上）。更改默认的后端引擎具有影响，因为无论是 GEOS 还是 s2 都在做出一些简化和不同的假设。请查看 https://r-spatial.github.io/sf/articles/sf7.html 以获取更多信息。 - Jindra Lacko

感谢您清晰的解释！我简要查看了sf v1.0的文档，并最终进入了您提供链接的同一页。 - IvanP

正如所料，因为我链接了{sf}维护者的官方通讯，:) - Jindra Lacko

Answer 2

我发现这个“无效的球面几何”问题经常出现。如果上面的s2::s2_rebuild()解决方案不起作用，通常适用于我的解决方案涉及投影和简化（稍微降低地图分辨率）。如果您的应用程序可以使用更低的分辨率，请尝试此方法。

library(tidyverse)
library(sf)
crs_N = 3995 #northern polar projection

# example of FAILING map - with bad spherical geometry.
m_RU <- rnaturalearthdata::countries50 %>% 
  st_as_sf() %>% 
  filter((admin %in% c("Russia") )) |> 
  st_as_s2()

在这个例子中，我选择了俄罗斯，因为它横跨日期变更线，这可能是其中一个挑战。我切换到北极极地投影，并将地图减少到10公里的分辨率（在这种情况下，5公里的分辨率不足够！）。

# with 2 extra lines the problem is gone
m_RU <- rnaturalearthdata::countries50 %>% 
  st_as_sf() %>% 
  filter((admin %in% c("Russia") )) |> 
  st_transform(crs = crs_N)   |>  
  st_simplify(dTolerance = 10000) |> # to get rid of duplicate vertex (reduce to 10km steps)
  st_as_s2()

- Scrope