osm_search_spatial returns invalid and inaccurate lon/lat data for a UK admin area

[francisbarton]

I asked a question about this on Stack Overflow but haven't got to the bottom of it.

I want to use nominatim to return a valid spatial object that I can plot as a polygon.

Here's a reprex of my problem with osm_search_spatial :

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 = "KtacQGOTApeFbfDhKaq5cGk8T16V4ioP")
class(ashfield)
#> [1] "list"

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

tmap::qtm(ashfield)
#> Warning: The shape ashfield is invalid. See sf::st_is_valid
#> Warning: Currect projection of shape ashfield unknown. Long-lat (WGS84) is
#> assumed.

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-171-248.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

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

# set CRS
st_crs(ashfield_sf) <- 4326
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...
tmap::qtm(ashfield_sf)
#> Warning: The shape ashfield_sf is invalid. See sf::st_is_valid

^{Created on 2020-02-26 by the reprex package (v0.3.0)}

Notice how the bbox section of the geometry is incorrect (below): the bbox_left number (a longitude) is using the latitude number for my area, and the bbox_bottom number (a latitude) is using a longitude number for my area. Hence the polygon plotted stretches from west Africa to Russia!

#> $ bbox_left    <fct> 53.0080617
#> $ bbox_top     <fct> 53.1714343
#> $ bbox_right   <fct> -1.3445928
#> $ bbox_bottom  <fct> -1.1642542

Here's a second reprex showing the result of the equivalent "direct" query to the API for comparison: I have not managed to get the result of this to be a plottable/mappable spatial object though (would appreciate some advice on this by the way!)

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(jsonlite)
library(purrr)
#> 
#> Attaching package: 'purrr'
#> The following object is masked from 'package:jsonlite':
#> 
#>     flatten

# try with direct query
ashfield_fromjson <- jsonlite::fromJSON("https://open.mapquestapi.com/nominatim/v1/search.php?key=KtacQGOTApeFbfDhKaq5cGk8T16V4ioP&format=json&q=ashfield&limit=1&polygon_geojson=1")
class(ashfield_fromjson)
#> [1] "data.frame"
# ^^^ I don't know how to turn this data frame into a plottable/tmappable spatial object.
# I have tried various things to try to use the geometry data in this result
# as a spatial object but I have not yet found the right way to do this

ashfield_coords <- ashfield_fromjson %>%
  pluck("geojson", "coordinates", 1) %>%
  matrix(., ncol = 2, dimnames = list(NULL, c("lon", "lat")))
# compare to the result from osm_search_spatial above
head(ashfield_coords)
#>            lon      lat
#> [1,] -1.344593 53.06326
#> [2,] -1.344409 53.06299
#> [3,] -1.344221 53.06276
#> [4,] -1.344161 53.06242
#> [5,] -1.344115 53.06231
#> [6,] -1.344068 53.06225

^{Created on 2020-02-26 by the reprex package (v0.3.0)}

[francisbarton]

This is what the district should look like (from Nominatim search): https://nominatim.openstreetmap.org/details.php?osmtype=R&osmid=154043&class=boundary i.e. the source data should be valid here!?