Почему пространственный фрейм данных, полученный из OSM / Nominatim, недействителен? - PullRequest
Новая
       121

Почему пространственный фрейм данных, полученный из OSM / Nominatim, недействителен?

2 голосов
/ 20 февраля 2020

Я хочу получить графические пространственные (граничные / полигональные) данные из OpenStreetMap, используя пакет nominatim R .

Запрос nominatim::osm_search_spatial работает успешно, но полученные данные Несмотря на то, что он обнаружен как «SpatialPolygonsDataFrame», он отображается как недопустимый, и я не могу получить его для построения с plot или tmap. Кажется, он также отсутствует CRS.

Я не уверен, пропустил ли я один или два шага или действительно что-то не так с данными, полученными из OSM (кажется маловероятным, но возможным) .

Редактировать: Кажется, есть проблема с широтой и долготой в полученных данных, поэтому этот вопрос был отредактирован с новым представлением. Похоже, что на самом деле, возможно, его следует задавать на ГИС-форуме, а не здесь на R.

Удалено предыдущее представление , показывающее лишние шаги и результаты (перенесено в суть здесь ) - это то, на что отвечал Юджин Чонг в своем первом ответе.

Вот мой новый, более сфокусированный, репрезентативный показ проблемы (ключ OSM API скрыт). Вы можете видеть, что многоугольник @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...

Создано в 2020-02-20 пакетом представитель (v0.3.0)

Следующий шаг: я иду попробовать прямую загрузку с nominatim, а не через пакет R. и посмотреть, что я получу.

1 Ответ

1 голос
/ 20 февраля 2020

Использование st_make_valid() из пакета lwgeom приводит к предупреждению go прочь:

Добавьте 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"

Что касается самой геометрии, то это все же всего лишь линия, проходящая от России до западного побережья Африки, а не полигон Эшфилда, как можно было бы ожидать.

И насколько почему Во-первых, я не уверен, что геометрия была неверной, но ошибка самопересечения предполагает, что многоугольник пересек сам себя в некоторой точке (как было бы изгибу ie). Смотрите этот пост, в котором более подробно обсуждаются действительные и недействительные геометрии: https://www.r-spatial.org/r/2017/03/19/invalid.html

...