This post is the 1st post of a series showcasing various rOpenSci
packages as if Maëlle were a birder trying to make the most of R in
general and rOpenSci in particular. Although the series use cases will
mostly feature birds, it’ll be the occasion to highlight rOpenSci’s
packages that are more widely applicable, so read on no matter what your
field is! Moreoever, each post should stand on its own.
A logical first step in this birder’s guide to rOpenSci is to use R to
find where to observe birds! In this blog post, we shall use rOpenSci
packages accessing open geographical data to locate and visualize
where to observe birds near a given location.
First of all, where are we now?
The Max Planck Institute for Ornithology (henceforth shortened to MPI),
where Maëlle will give a talk is located at Am Obstberg 1 78315
Radolfzell. Let’s geolocate it using rOpenSci’s
package that interfaces the
OpenCage Geocoder, a commercial service
based on open data. When choosing to
get only one result via
limit = 1, one gets what the API considers to
be the best one.
mpi <- opencage::opencage_forward("Am Obstberg 1 78315 Radolfzell", limit = 1)$results class(mpi)
##  "tbl_df" "tbl" "data.frame"
##  "annotations.DMS.lat" "annotations.DMS.lng" ##  "annotations.MGRS" "annotations.Maidenhead" ##  "annotations.Mercator.x" "annotations.Mercator.y"
This gets us Am Obstberg 1, 78315 Radolfzell, Germany (
which is in ?? (
mpi$annotations.flag gets us a flag!).
Birding in a bird hide?
Where to find a bird hide?
You can most certainly go birding anywhere, but if you can find a bird
hide (or bird blind depending on the English you speak), it might be a
very appropriate observation point. Now that we know where the MPI is,
we can look for bird hide(s) in the vicinity. For that, we shall use
osmdata package by
Mark Padgham and collaborators! Note that incidentally, Mark did his
This package is an interface to OpenStreetMap’s Overpass
API. OpenStreetMap is
a collective map of the world. It contains information about towns’
limits, roads, placenames… but also tags of everything, from bars as
seen in this blog post to
trees. You can
browse existing features via OpenStreetMap’s
wiki. Some parts of the
world are better mapped than others depending on the local OpenStreetMap
community. Actually, OpenCage’s blog features an interesting series of
To look for a bird hide, we first create a bounding box of 10km around
the MPI, using rOpenSci’s
bbox <- bbox::lonlat2bbox(mpi$geometry.lng, mpi$geometry.lat, dist = 10, method = "lawn")
We then use the key and value associated with bird hides in
respectively leisure and bird_hide.
## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
library("magrittr") (results <- opq(bbox = bbox) %>% add_osm_feature(key = 'leisure', value = 'bird_hide') %>% osmdata_sf ())
## Object of class 'osmdata' with: ## $bbox : 47.6865,8.8753,47.8494,9.1177 ## $overpass_call : The call submitted to the overpass API ## $timestamp : [ Thu 5 Jul 2018 08:06:35 ] ## $osm_points : 'sf' Simple Features Collection with 1 points ## $osm_lines : 'sf' Simple Features Collection with 0 linestrings ## $osm_polygons : 'sf' Simple Features Collection with 0 polygons ## $osm_multilines : 'sf' Simple Features Collection with 0 multilinestrings ## $osm_multipolygons : 'sf' Simple Features Collection with 0 multipolygons
## leisure geometry ## 5004940425 bird_hide 8.920901, 47.741569
Yay, we now know where to find a bird hide not too far from the MPI!
Visualizing our location and the bird hide
So one could enter the coordinates of that bird hide in one’s favourite
mapping software or app but to show you where the bird hide is we can
actually step back and use
osmplotr, another package
contributed to rOpenSci by Mark Padgham!
osmplotr works is letting you create a basemap, on which
you’ll add different layers extracted from OpenStreetMap using
osmdata functions directly. Its
strengths are therefore the use of open data, and the customization of
what you’re using as background!
Let’s create a basemap for our bounding box, and then add roads and
buildings to it.
## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
bbox <- get_bbox(bbox) dat_B <- extract_osm_objects (key = 'building', bbox = bbox) dat_H <- extract_osm_objects (key = 'highway', bbox = bbox) map0 <- osm_basemap(bbox = bbox, bg = 'gray20') %>% add_osm_objects (dat_B, col = 'gray40') %>% add_osm_objects (dat_H, col = 'gray80') map0 %>% add_axes() %>% print_osm_map (filename = 'map_a1.png', width = 600, units = 'px', dpi = 72) library("magrittr") magick::image_read('map_a1.png') %>% magick::image_annotate("Map data © OpenStreetMap contributors", color = "white", boxcolor = "black", size = 15, gravity = "south")
Quite pretty! The lakes can be seen because of the absence of roads and
buildings on them.
Now, let’s plot the bird hide and the MPI on them. Since we used
osmdata::osmdata_sf, we had gotten the data in a receivable class,
points_map <- add_osm_objects(map0, results$osm_points, col = 'salmon', size = 5)
For plotting the MPI, we’ll convert
opencage output to an
with the same coordinate reference system as the OpenStreetMap data
coords <- data.frame(lon = mpi$geometry.lng, lat = mpi$geometry.lat) crs <- sf::st_crs(results$osm_points) mpi_sf <- sf::st_as_sf(coords, coords = c("lon", "lat"), crs = crs) points_map <- add_osm_objects(points_map, mpi_sf, col = 'white', size = 5)
We can now visualize both points on the map, the MPI in white and the
bird hide in salmon, South-West from the MPI.
points_map %>% add_axes() %>% print_osm_map (filename = 'map_a2.png', width = 600, units = 'px', dpi = 72) magick::image_read('map_a2.png') %>% magick::image_annotate("Map data © OpenStreetMap contributors", color = "white", boxcolor = "black", size = 15, gravity = "south")
Aha, now we see where the bird hide is, fantastic! But as Mark noted,
birds can actually be observed from other places.
Birding where birds should be?
Birds are most likely to be found where water lies close to natural
areas, and we can translate this to R code! We shall get all water and
(separately) all non-watery natural areas and find the shortest
distances between them before plotting the results using
First, let’s get all water and (separately) all non-watery natural
dat <- opq(bbox = bbox) %>% add_osm_feature(key = 'natural') %>% osmdata_sf (quiet = FALSE)
## Issuing query to Overpass API ... ## Rate limit: 2 ## Query complete! ## converting OSM data to sf format
indx_W <- which (dat$osm_polygons$natural %in% c ("water", "wetland")) indx_N <- which (!dat$osm_polygons$natural %in% c ("water", "wetland")) xy_W <- sf::st_coordinates (dat$osm_polygons [indx_W, ]) xy_N <- sf::st_coordinates (dat$osm_polygons [indx_N, ])
Then we use Mark’s
package to get pairwise distances
between all points, find minima, and make a data.frame to submit to
add_osm_surface(). We have 7068 watery points and 10065 non-watery
points so the speed of
geodist is crucial here!
t1 <- Sys.time() d <- geodist::geodist (xy_W, xy_N) # so fast!!! Sys.time() - t1
## Time difference of 13.57983 secs
d1 <- apply (d, 1, min) d2 <- apply (d, 2, min) xy <- cbind (rbind (xy_W, xy_N), c (d1, d2)) xy <- xy [, c (1, 2, 5)] colnames (xy) <- c ("x", "y", "z") xy <- xy [!duplicated (xy), ]
Finally we plot the results on the roads we had gotten earlier, although
we do not recommend staying on the middle of a road to observe birds! We
re-add the points corresponding to the MPI and bird hide after the
surface coloring. With
osmplotr, order matters because layers are
added on top of each other. Note that plotting the distances takes a
# colorblind-friendly palette! cols <- viridis::viridis_pal (direction = -1) (30) add_osm_surface (map0, dat_H, dat = xy, col = cols) %>% add_axes() %>% add_colourbar(cols = cols, zlim = range(xy[,"z"])) %>% add_osm_objects(mpi_sf, col = 'white', size = 5) %>% add_osm_objects(results$osm_points, col = 'white', size = 5)%>% print_osm_map (filename = 'map_a3.png', width = 600, units = 'px', dpi = 72) magick::image_read('map_a3.png') %>% magick::image_annotate("Map data © OpenStreetMap contributors", color = "white", boxcolor = "black", size = 15, gravity = "south")
On the map, the yellower/lighter a road is, the better it is to observe
birds according to Mark’s assumption that birds are most likely to be
found where water lies close to natural areas. With this metric, the MPI
itself is not too badly located, which after all is not too surprising
for an MPI of ornithology. Maybe one should just walk to the one of
the nearest lakes to meet some birds.
Open geographical data in R
In this post we showcased three rOpenSci packages helping you use open
geographical data in R:
osmplotr, therefore mostly
making use of the awesome OpenStreetMap data (The OpenCage Geocoder uses
a bit more, but it only warrants attributing
OpenStreetMap). We therefore added
the annotation “Map data © OpenStreetMap contributors” to all maps of
this post using rOpenSci’s
magick package. You might also have noticed
in the code earlier that both
osmplotr have start-up
messages indicating the data origin and licence.
Other R packages for spatial analyses
We also used the
package for ultra lightweight,
ultra fast calculation of geo distances. This package is developed in
the hypertidy GitHub organization which
is a good place to find useful R geospatial packages. Other good
resources for geospatial analyses with R include the r-spatial.org
website and this great book by Robin
Lovelace, Jakub Nowosad and Jannes
Muenchow, and more links
presented in this blog post of Steph
More birding soon!
Stay tuned for the next post in this series, about getting bird
observation data in R! In the meantime, happy birding!