querying features that enclose a point

[barryrowlingson]

The OpenStreetMap web page has two query outputs - features near the query point and features enclosing the query point. The "enclosing" features list means you can click somewhere and get what's under that point. This can return multiple features, like a lake, the county the lake is in, the country the county is in and soon.

I reverse-engineered what the web page was querying to come up with the following function:

getosm <- function(lat, long, tag){
    s = paste0('https://overpass-api.de/api/interpreter?data=%5Bout%3Ajson%5D%3Bis_in(',lat,'%2C',long,')-%3E.a%3Bway(pivot.a)%3Bout+tags+bb+geom%3Bout+ids+geom%3Brelation(pivot.a)%5B',tag,'%5D%3Bout+geom%3B')
    rjson::fromJSON(paste(readLines(url(s)),collapse="\n"))
}

which I used to answer a GIS StackExchange question here: https://gis.stackexchange.com/questions/345994/get-waterbody-info-from-base-map-in-leaflet

I don't pretend to really understand the OverpassAPI syntax thoroughly but here it is in plain text:

data=[out:json];
is_in(lat,long)->.a;
way(pivot.a);
out tags bb geom;
out ids geom;
relation(pivot.a)[tag=value];
out geom;

where lat,long is the point coordinate and tag=value is something like natural=water if the user wants to search for lakes, for example. I'm not sure if the multiple out lines are necessary. I don't know if the way lines needs a tag filter. I don't understand the magic behind the pivot.a thing. All I understand (see the GIS SE question) is that this R function returns the lakes under a point.

Currently my code returns the raw JSON, so its a bit useless in R without some post-processing, but there must be code in osmdata to do all that. Could this query of finding features enclosing a point be included in osmdata?

[mpadge]

Thanks @barryrowlingson, that is a really good idea that i hadn't thought of, the good news for which is, yes, that would be really easy to incorporate within osmdata. Similar to the previous issue on extracting by OSM ID values, this functionality would not really rest neatly within any current functional structures, so would likely best have its own, separate functionality. Something as simple as osm_nearby() should do, with an output directly able to be piped to the extraction functions:

osm_nearby(x, y) %>%
   osmdata_sf/sp/sc()

but also, in line with your example, accept additional defining features,

osm_nearby(x, y) %>%
    add_osm_feature(tag = "natural", value = "water") %>%
   osmdata_sf/sp/sc()

That function could then also immediately accept an additional osm_id parameter, rather than coordinates:

osm_nearby(osm_id = <long long int>) %>%
    osmdata_sf/sp/sc()

The only real potential incongruity there would be that add_osm_feature is not really adding as such, rather restricting the result to the specified key-value pair, but that's a minor semantic issue that could be addressed once functionality is in place. Any comments or further suggestions welcome in the meantime - I'll try to find some time early in the new year to implement that. Thanks for the suggestion, and as always, for using the package!

[mpadge]

BTW @barryrowlingson i don't have enough rep points to comment on gis.se, but osmar has long been completely and irrevocably broken. The API it used to reply on no longer even exists, and the package does nothing. It also contains no tests whatsoever, and so this failure just slides on through the CRAN servers ... alas. The package should be removed and forgotten about.

[barryrowlingson]

I don't think "nearby" is the right word for this, since that implies a location just outside a polygon would also return the polygon. The query as written finds things that enclose the point, and hence I think it can't return point features or line features (nodes/ways).

Its more like osm_enclosing(x,y) or maybe osm_whats_at(x,y)

You could probably also do a separate osm_nearby - the query could be obtained by looking at what openstreetmap.com generates for a click for that list - here's that:

data | [timeout:10][out:json];(node(around:15,53.94542,-2.52017);way(around:15,53.94542,-2.52017););out+tags+geom(53.94429,-2.52499,53.94605,-2.51902);relation(around:15,53.94542,-2.52017);out+geom(53.94429,-2.52499,53.94605,-2.51902);

which returns things within 15 metres of that location. Wow this is nice...

[mpadge]

@barryrowlingson am now finally getting on to this. The above commits implement the example given in the new opq_enclosing function which does this:

library (osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

lat <- 54.33601
lon <- -3.07677
key <- "natural"
value <- "water"
q <- opq_enclosing (lon, lat, key, value) %>%
    opq_string ()
q
#> [1] "[out:xml][timeout:25];\n(\nis_in(54.33601,-3.07677)->.a;relation(pivot.a) [\"natural\"=\"water\"];);\n(._;>;);\nout;"
x <- osmdata_sf (q)
x
#> Object of class 'osmdata' with:
#>                  $bbox : 
#>         $overpass_call : The call submitted to the overpass API
#>                  $meta : metadata including timestamp and version numbers
#>            $osm_points : 'sf' Simple Features Collection with 1241 points
#>             $osm_lines : NULL
#>          $osm_polygons : 'sf' Simple Features Collection with 2 polygons
#>        $osm_multilines : NULL
#>     $osm_multipolygons : 'sf' Simple Features Collection with 1 multipolygons

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

way-type enclosures can be obtained like this:

library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

lat <- 54.04827
lon <- -2.82288
key <- value <- NULL
q <- opq_enclosing (lon, lat, key, value, enclosing = "way") %>%
    opq_string ()
q
#> [1] "[out:xml][timeout:25];\n(\nis_in(54.04827,-2.82288)->.a;way(pivot.a););\n(._;>;);\nout;"
x <- osmdata_sf (q)
x
#> Object of class 'osmdata' with:
#>                  $bbox : 
#>         $overpass_call : The call submitted to the overpass API
#>                  $meta : metadata including timestamp and version numbers
#>            $osm_points : 'sf' Simple Features Collection with 51 points
#>             $osm_lines : NULL
#>          $osm_polygons : 'sf' Simple Features Collection with 1 polygons
#>        $osm_multilines : NULL
#>     $osm_multipolygons : NULL

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

As expected, the default enclosing = "relation" generally returns multipolygon objects, while enclosing = "way" will always return polygon objects only. Two questions:

1. What do you think of that implementation? Any comments?
2. I think it's best to separate the relation and way queries, and to leave the default as relation only, because inadvertent application of enclosing = "way" can return huge amounts of data, and is probably not wise. That said, it would of course be possible to combine the two and remove this parameter, but that would risk submitting unnecessarily large requests to the overpass server. Thoughts?

---

TODO

• implement osm_nearby() as you describe