3-layers Open Street Map Data & Rayshader Map

rm(list = ls())

1 Loading library

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.8
## ✓ tidyr   1.2.0     ✓ stringr 1.4.0
## ✓ readr   2.1.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUE
library(rayshader) # To install rayshader, must install XQuartz first: https://github.com/tylermorganwall/rayshader/issues/86
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
library(elevatr)

2 Create Bounding Box

med_bbox <- st_bbox(c(xmin = -75.676575, xmax = -75.462715, 
                      ymin = 6.113246, ymax = 6.372624),
                    crs = 4326)

med_bbox_df <- data.frame(x = c(-75.676575, -75.462715),
                       y = c(6.113246, 6.372624))


extent_zoomed <- raster::extent(med_bbox)

3 Get elevation data

prj_dd <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

elev_med <- get_elev_raster(med_bbox_df, prj =  prj_dd, z = 12, clip = "bbox")
## Mosaicing & Projecting
## Clipping DEM to bbox
## Note: Elevation units are in meters.
elev_med_mat <- raster_to_matrix(elev_med)

4 Get road data

med_roads <- med_bbox %>%
  opq() %>%
  add_osm_feature("highway",
                  c("motorway",
                    "trunk",
                    "primary", 
                    "secondary")) %>%
  osmdata_sf()
## Error in curl::curl_fetch_memory(url, handle = handle): HTTP/2 stream 0 was not closed cleanly: PROTOCOL_ERROR (err 1)
## Request failed [ERROR]. Retrying in 1 seconds...
med_road_lines <- med_roads$osm_lines

5 Get river data

med_rivers <- med_bbox %>%
  opq() %>%
  add_osm_feature("waterway",
                  c("river")) %>%
  osmdata_sf()

med_river_lines <- med_rivers$osm_lines

6 Get building data

med_buildings <- med_bbox %>%
  opq() %>%
  add_osm_feature("building") %>%
  osmdata_sf()
## Error in curl::curl_fetch_memory(url, handle = handle): HTTP/2 stream 0 was not closed cleanly: PROTOCOL_ERROR (err 1)
## Request failed [ERROR]. Retrying in 1.4 seconds...
## Error in curl::curl_fetch_memory(url, handle = handle): HTTP/2 stream 0 was not closed cleanly: PROTOCOL_ERROR (err 1)
## Request failed [ERROR]. Retrying in 3.5 seconds...
med_building_polygons <- med_buildings$osm_polygons

7 Create basemap

base_map <- elev_med_mat %>% 
  height_shade() %>% 
  add_overlay(
    sphere_shade(elev_med_mat,
                 texture = rayshader::create_texture(
                   lightcolor = "#b8ff78",
                   shadowcolor = "#193600",
                   leftcolor = "#80d453",
                   rightcolor = "#80d453",
                   centercolor = "#568a27"),
                 sunangle = 0,
                 colorintensity = 5)
  )

8 Add roads, rivers, and building

final_map <- base_map %>%
  add_overlay(
    generate_line_overlay(
      med_road_lines, extent = extent_zoomed,
      linewidth = 3, color = "white",
      heightmap = elev_med_mat
    )) %>%
  add_overlay(
    generate_line_overlay(
      med_river_lines, extent = extent_zoomed,
      linewidth = 10, color = "lightblue",
      heightmap = elev_med_mat
    )) %>%
  add_overlay(
    generate_polygon_overlay(
      med_building_polygons, extent = extent_zoomed,
      linewidth = 0, palette = "red",
      heightmap = elev_med_mat
    )) %>%
  plot_3d(elev_med_mat, zscale = 10, fov = 0, 
          theta = -45, zoom = .5, phi = 30, 
          windowsize = c(1000,800))


# render_snapshot(
#   filename = '14-3-dimensional/plot.png',
#   title_text = "Medellín, Colombia\nFuente: OSM | Elaborado por Camilo Martínez (@camartinezbu)",
#   title_offset = c(20, 20),
#   title_color = "black",
#   title_size = 25,
#   title_font = "Assistant",
#   title_position = "north"
# )