First of all, where are we now?

mpi <- opencage::opencage_forward("Am Obstberg 1 78315 Radolfzell", limit = 1)$results
class(mpi)
## [1] "tbl_df"     "tbl"        "data.frame"
head(names(mpi))
## [1] "annotations.DMS.lat"    "annotations.DMS.lng"   
## [3] "annotations.MGRS"       "annotations.Maidenhead"
## [5] "annotations.Mercator.x" "annotations.Mercator.y"

Birding in a bird hide?

Where to find a bird hide?

bbox <- bbox::lonlat2bbox(mpi$geometry.lng,
                          mpi$geometry.lat,
                          dist = 10, method = "lawn")
library("osmdata")
## 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 : [ Tue 3 Sep 2018 10:29:46 ]
##            $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
results$osm_points
##              leisure            geometry
## 5004940425 bird_hide 8.920901, 47.741569

Visualizing our location and the bird hide

library("osmplotr")
## 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")

points_map <- add_osm_objects(map0, results$osm_points, 
                              col = 'salmon',
                              size = 5)
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)
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")

Birding where birds should be?

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, ])
t1 <- Sys.time()
d <- geodist::geodist (xy_W, xy_N)
# so fast!!!
Sys.time() - t1
## Time difference of 21.38027 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), ]
# 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")

Conclusion

Open geographical data in R

Other R packages for spatial analyses

More birding soon!