1

I am plotting points on a map using Leaflet, and would like to collect any data I can concerning the water bodies that encompass these points. I have been making use of polygons from HydroSHEDS (https://hydrosheds.org) to identify waterbody outlines, but this is not a complete list of waterbodies (The default base map used in leaflet seems to have polygon info on many more lakes).

As an example, if I make a map like this:

library(leaflet)

Lat<-c(43.847863,43.847199,43.846288,43.852605)
Lon<-c(-74.886135,-74.886111,-74.897731,-74.899169)
LakeID<-c(1,2,3,4)


leaflet() %>%
  addCircleMarkers(lat = Lat, lng = Lon,
                   color = "blue",
                   radius = 2)%>%
  addTiles() 

For points that fall in a water body, I'd like to get info on the polygon used for that lake, and the name of the lake as well. I can see that the base map has the larger and smaller waterbodies colored blue, and has the name for the larger one "Twitchell Lake", so I assume that data can be queried in some way.

Is there a way I can access this data from the basemap so I could find out:

Polygon used for individual waterbodies Name of waterbody

2

If you look at the credit line at the bottom of the map you'll see that the basemap is data from OpenStreetMap (OSM). There is a network API you can call to query the underlying vector data that is used to draw the map.

The R package osmar has some functionality to read OSM data but I can't get it to work at the moment. So I wrote a quick 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"))
}

This function takes a single latitude and longitude, together with a tag string you can use to select elements. You probably want to select elements with the "natural=water" tag. So for your first point:

w1 = getosm(Lat[1],Lon[1],"natural=water")

The returned object is a complex list object. For example the name of the first element (since it can possibly find several features at a location) is:

> w1$elements[[1]]$tags$name
[1] "Twitchell Lake"
> 

The geometry defining the lake is stored in w1$elements[[1]]$members object. It has five elements:

> length(w1$elements[[1]]$members)
[1] 5

which (probably) define the main shoreline and maybe some islands inside the lake. There might be code in osmar to construct spatial features from this data, or some other conversion code somewhere.

If there's no water at the queried point then the returned elements will be length zero:

> w3 = getosm(Lat[3],Lon[3],"natural=water")
> length(w3$elements)
[1] 0
> 

R could benefit from a package to do queries using this API, especially since osmar looks broken. I might investigate the osmar problems further.

1

Another solution using osmdata package. I decided to use tmap for plotting just because I'm more familiar with its syntax but there is not real difference between interactive version of tmap and leaflet.

# packages
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(tmap)

tmap_mode("view")
#> tmap mode set to interactive viewing

# point data
Lat<-c(43.847863, 43.847199, 43.846288, 43.852605)
Lon<-c(-74.886135, -74.886111, -74.897731, -74.899169)

I need to define a bounding box to create a query for OSM (read the help page of ?osmdata::opq and ?osmdata::getbb for more details) so I manually create a bbox based on those points. You should manually expand that bbox if you want to include more features.

# define a bbox for those points
my_points <- st_sfc(st_multipoint(cbind(Lon, Lat)), crs = 4326)
my_bbox <- st_bbox(my_points)

Now I have to create an Overpass Query specifying the appropriate keys and values. I used the same query as @Spacedman but you can read mode details here to modify it.

# define tags
my_lakes <- opq(my_bbox) %>% 
  add_osm_feature(key = "natural", value = "water") %>% 
  osmdata_sf()

This is the result. You should read the vignettes of osmdata package for a more detailed explanation of that output.

my_lakes
#> Object of class 'osmdata' with:
#>                  $bbox : 43.846288,-74.899169,43.852605,-74.886111
#>         $overpass_call : The call submitted to the overpass API
#>                  $meta : metadata including timestamp and version numbers
#>            $osm_points : 'sf' Simple Features Collection with 128 points
#>             $osm_lines : 'sf' Simple Features Collection with 4 linestrings
#>          $osm_polygons : 'sf' Simple Features Collection with 2 polygons
#>        $osm_multilines : NULL
#>     $osm_multipolygons : 'sf' Simple Features Collection with 1 multipolygons

Let's check now the polygons and multipolygons objects.

lakes_polygon <- my_lakes$osm_polygons
lakes_polygon
#> Simple feature collection with 2 features and 3 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -74.90088 ymin: 43.84064 xmax: -74.87768 ymax: 43.85952
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>              osm_id        name natural                       geometry
#> 198403657 198403657        <NA>    <NA> POLYGON ((-74.88185 43.8501...
#> 222234302 222234302 Oswego Pond   water POLYGON ((-74.90039 43.853,...

lakes_multipolygon <- my_lakes$osm_multipolygons
lakes_multipolygon
#> Simple feature collection with 1 feature and 4 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -74.90088 ymin: 43.84064 xmax: -74.87768 ymax: 43.85952
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>          osm_id           name natural         type
#> 4056964 4056964 Twitchell Lake   water multipolygon
#>                               geometry
#> 4056964 MULTIPOLYGON (((-74.88166 4...

The polygon object consists of 2 POLYGONS while the multipolygon sf object consists of just 1 MULTIPOLYGON. We can see that both outputs include a variable called name that specify the name of the corresponding object according to OpenStreetMap data. Now we can also plot these results. There are a few problems with leaflet and osmdata (i.e here and here) but we can easily solve them.

# and plot the result
names(lakes_multipolygon$geometry) <- NULL
names(lakes_multipolygon$geometry[[1]][[1]]) <- NULL

tm_shape(lakes_multipolygon) + 
  tm_polygons(col = "darkred", alpha = 0.5) + 
tm_shape(lakes_polygon) + 
  tm_polygons(col = "darkred") + 
tm_shape(my_points) + 
  tm_dots(col = "darkblue") + 
tm_basemap("OpenStreetMap")

Created on 2019-12-28 by the reprex package (v0.3.0)

The MULTIPOLYGON object represent the Twitchell lake, the first polygon represent a small island inside the Twitchell Lake while the second polygon represent the other lake.

  • 1
    I looked at doing a bounding-box query but the questioner asked for polygons containing a point. If you want to do it by bounding box you first have to decide how big a box you need, and that needs adapting since a bounding box query needs to be large enough to catch some of the surrounding polygon, and also do a point-in-polygon test to filter out things that don't have the point in them (using st_intersects). – Spacedman Dec 28 '19 at 12:03
  • Oh ok, sorry, I didn't understand the OP question. – agila Dec 28 '19 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.