r - 为什么从 OSM/Nominatim 获取的空间数据框无效？

Question

我想使用 nominatim R package 从 OpenStreetMap 获取可绘制的空间(边界/多边形)数据.

nominatim::osm_search_spatial 请求似乎成功运行，但生成的数据尽管被发现为“SpatialPolygonsDataFrame”，但显示为无效，我无法使用 plot 或 tmap。它似乎也缺少 CRS。

我不确定我是否遗漏了一两步，或者从 OSM 检索的数据是否确实有问题(似乎不太可能，但有可能)。

编辑:接收到的数据中的纬度和经度似乎有问题，因此已使用新的 reprex 编辑了此问题。事实上，看起来可能需要在 GIS 论坛上而不是在 R 中提问。

删除之前的 reprex 显示多余的步骤和结果(移至要点 here)——这就是 Eugene Chong 在他的第一个回复中的回应。

这是我新的、更集中的、显示问题的 reprex(隐藏了 OSM API key )。您可以看到多边形 @coords 从第一个点 -1,-1 跃升到第二个点的 53,53，因此绘制了对角线。

library(nominatim)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
#> Nominatim Usage Policy: http://wiki.openstreetmap.org/wiki/Nominatim_usage_policy
#> MapQuest Nominatim Terms of Use: http://info.mapquest.com/terms-of-use/
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(tmap)
library(tibble)

# get OSM search results for Ashfield district (UK)
ashfield <- nominatim::osm_search_spatial("Ashfield", limit = 1, key = $my_osm_api_key)
class(ashfield)
#> [1] "list"

# extract SPDF from list
ashfield <- ashfield[[1]]
class(ashfield)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"

# Convert to an SF object and try again
ashfield_sf <- sf::st_as_sf(ashfield)
class(ashfield_sf)
#> [1] "sf"         "data.frame"

# set CRS (thanks to Eugene Chong)
st_crs(ashfield_sf) <- 4326
tmap::qtm(ashfield_sf)
#> Warning: The shape ashfield_sf is invalid. See sf::st_is_valid

glimpse(ashfield@data)
#> Observations: 1
#> Variables: 15
#> $ place_id     <chr> "186877616"
#> $ licence      <chr> "Data © OpenStreetMap contributors, ODbL 1.0. https://...
#> $ osm_type     <chr> "relation"
#> $ osm_id       <chr> "154043"
#> $ lat          <dbl> 53.08977
#> $ lon          <dbl> -1.251877
#> $ display_name <chr> "Ashfield, Nottinghamshire, East Midlands, England, Un...
#> $ class        <chr> "boundary"
#> $ type         <chr> "administrative"
#> $ importance   <dbl> 0.2116014
#> $ icon         <chr> "http://ip-10-98-176-55.mq-us-east-1.ec2.aolcloud.net/...
#> $ bbox_left    <fct> 53.0080617
#> $ bbox_top     <fct> 53.1714343
#> $ bbox_right   <fct> -1.3445928
#> $ bbox_bottom  <fct> -1.1642542
head(ashfield@polygons[[1]]@Polygons[[1]]@coords)
#>           [,1]      [,2]
#> [1,] -1.344593 -1.344409
#> [2,] 53.063537 53.063260
#> [3,] 53.064985 53.063764
#> [4,] 53.065520 53.065521
#> [5,] 53.065553 53.065526
#> [6,] 53.065725 53.065656
ashfield_sf$geometry
#> Geometry set for 1 feature 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -1.344593 ymin: -1.344593 xmax: 53.17143 ymax: 53.17142
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#> POLYGON ((-1.344593 -1.344409, 53.06354 53.0632...

^{由 reprex package 创建于 2020-02-20 (v0.3.0)}

下一步:我将尝试从 nominatim 直接下载，而不是通过 R 包，看看我得到了什么。

Answer 1

使用 lwgeom 包中的 st_make_valid() 可以使警告消失:

添加 WGS84 作为 crs。

ashfield_sf <- sf::st_as_sf(ashfield)
st_crs(ashfield_sf) <- 4326

sf::st_is_valid(ashfield_sf, reason = TRUE)

[1] "Self-intersection[53.0709899483331 53.0709561246412]"

使用st_make_valid()

library(lwgeom)
ashfield_sf_2 <- st_make_valid(ashfield_sf)
sf::st_is_valid(ashfield_sf_2, reason = TRUE)

[1] "Valid Geometry"

不过，就几何本身而言，它仍然只是一条从俄罗斯延伸到非洲西海岸的线，而不是预期的 Ashfield 多边形。

至于为什么几何首先无效，我不确定，但自相交错误表明多边形在某个点与自身交叉(就像领结一样)。请参阅这篇文章，其中对有效和无效的几何图形进行了更多讨论:https://www.r-spatial.org/r/2017/03/19/invalid.html