3d-map

A short description of the post.

The R package rayshader uses different hill shading techniques by calculating a colour for each point on the surface using:

• A direct elevation-to-colour mapping: height_shade()
• The surface normals and hemispherical UV mapping: sphere_shade()
• Lambert Shading Mapping by calculating the dot product between light direction and the surface normal vector: lamb_shade()
• shadows calculated by Leland Brown’s “texture shading” method: texture_shade(). This better defines ridges and drainage networks
• Ambient occlusion layer darkens the valleys to account for less scattered atmospheric light reaching the valley floor: ambient_shade()

load the reuired packages

library(dplyr)
library(ggplot2)
library(osmdata)
library(raster)
library(rayshader)
library(sf)
library(tmap)

1. load the elevation data for area of interest and convert to a matrix

dem <- raster("D:/spatial-data/dem/n30_e076_1arc_v3.tif") |>
  merge(raster("D:/spatial-data/dem/n30_e077_1arc_v3.tif")) |> 
  crop(extent(76.8, 77.1, 30.6, 30.9))

dem_mat <- raster_to_matrix(dem)

2. Build a base map with rayshader

dem_mat |> 
  height_shade() |> 
  plot_map()

3. add a spherical aspect shading colour

dem_mat |> 
  height_shade() |> 
  add_overlay(
    sphere_shade(dem_mat, texture = "desert", zscale=4, colorintensity = 5), 
    alphalayer=0.5) |>
  plot_map()

4. overlay a standard hillshade or Lambertian (cosine) shading

dem_mat |> 
  height_shade() |> 
  add_overlay(
    sphere_shade(dem_mat, texture = "desert", zscale = 4, colorintensity = 5), 
    alphalayer = 0.5) |>
  add_shadow(lamb_shade(dem_mat, zscale = 6), 0) |>
  plot_map()

5. add shadows calculated by Leland Brown’s “texture shading” method

dem_mat |> 
  height_shade() |> 
  add_overlay(
    sphere_shade(dem_mat, texture = "desert", zscale = 4, colorintensity = 5), 
    alphalayer = 0.5) |>
  add_shadow(lamb_shade(dem_mat, zscale = 6), 0) |>
  add_shadow(texture_shade(dem_mat, detail = 8/10, contrast = 9, brightness = 11), 
             0.1) |>
  plot_map()

6. add an ambient occlusion layer

base_map <- dem_mat |> 
  height_shade() |> 
  add_overlay(
    sphere_shade(dem_mat, texture = "desert", zscale = 4, colorintensity = 5), 
    alphalayer = 0.5) |>
  add_shadow(lamb_shade(dem_mat, zscale = 6), 0) |>
  add_shadow(ambient_shade(dem_mat), 0) %>%
  add_shadow(texture_shade(dem_mat, detail = 8/10, contrast = 9, brightness = 11), 
             0.1)
base_map |> plot_map()

7. extract major roads from osm and add them

# define a bounding box for area of interest
osm_bbox <- c(76.8, 30.6, 77.1, 30.9) # xmin, ymin, xmax, ymax

# extract highway from osm
#dem_highway <- opq(osm_bbox) |> 
  #add_osm_feature(key = "highway") |> 
  #osmdata_sf()

# filter major roads
#trunk_roads <- dem_highway$osm_lines |> filter(highway == "trunk")
#secondary_roads <- dem_highway$osm_lines |> filter(highway == "secondary")
#st_write(trunk_roads, "data/trunk_road.gpkg")
#st_write(secondary_roads, "data/secondary_road.gpkg")

trunk_roads <- st_read("data/trunk_road.gpkg")

Reading layer `trunk_road' from data source 
  `D:\R\abhikumar86.github.io\_posts\2021-07-21-3d-map\data\trunk_road.gpkg' 
  using driver `GPKG'
Simple feature collection with 188 features and 32 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 76.78954 ymin: 30.58918 xmax: 77.1175 ymax: 30.91891
Geodetic CRS:  WGS 84

secondary_roads <- st_read("data/secondary_road.gpkg")

Reading layer `secondary_road' from data source 
  `D:\R\abhikumar86.github.io\_posts\2021-07-21-3d-map\data\secondary_road.gpkg' 
  using driver `GPKG'
Simple feature collection with 126 features and 32 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 76.75612 ymin: 30.59039 xmax: 77.12317 ymax: 30.9053
Geodetic CRS:  WGS 84

# remove dem_highway
rm(dem_highway)

# prepare the road layer for overlaying
roads_layer <- generate_line_overlay(trunk_roads, extent = extent(76.8, 77.1, 30.6, 30.9),
                          heightmap = dem_mat,
                          linewidth = 5, color = "white") |>
  add_overlay(
    generate_line_overlay(secondary_roads, extent = extent(76.8, 77.1, 30.6, 30.9),
                          heightmap = dem_mat,
                          linewidth = 2.5, color = "white"))
# add to base map
base_map |>
  add_overlay(roads_layer) |>
  plot_map()

9. similarly extract waterways from osm and add them to base-map

# define a bounding box for area of interest
osm_bbox <- c(76.8, 30.6, 77.1, 30.9) # xmin, ymin, xmax, ymax

# extract waterway (rivers) from osm
#dem_water <- opq(osm_bbox) |> 
  #add_osm_feature(key = "waterway") |> 
  #osmdata_sf()

# filter major river and streams
#rivers <- dem_water$osm_lines |> filter(waterway == "river")
#streams <- dem_water$osm_lines |> filter(waterway != "river")
#st_write(rivers, "data/rivers.gpkg")
#st_write(streams, "data/streams.gpkg")

rivers <- st_read("data/rivers.gpkg")

Reading layer `rivers' from data source 
  `D:\R\abhikumar86.github.io\_posts\2021-07-21-3d-map\data\rivers.gpkg' 
  using driver `GPKG'
Simple feature collection with 30 features and 16 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 74.56863 ymin: 29.4721 xmax: 77.11868 ymax: 30.91332
Geodetic CRS:  WGS 84

streams <- st_read("data/streams.gpkg")

Reading layer `streams' from data source 
  `D:\R\abhikumar86.github.io\_posts\2021-07-21-3d-map\data\streams.gpkg' 
  using driver `GPKG'
Simple feature collection with 41 features and 16 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 73.5626 ymin: 29.22947 xmax: 77.07197 ymax: 30.8885
Geodetic CRS:  WGS 84

# remove dem_water
#rm(dem_water)

# prepare the river and stream layer for overlaying
stream_layer <- generate_line_overlay(rivers, extent = extent(76.8, 77.1, 30.6, 30.9),
                          heightmap = dem_mat, linewidth = 8, color = "skyblue2") |>
  add_overlay(
    generate_line_overlay(streams, extent = extent(76.8, 77.1, 30.6, 30.9),
                          heightmap = dem_mat, linewidth = 4, color = "skyblue"))
# add to base map
base_map |>
  add_overlay(stream_layer, alphalayer = 0.8) |>
  add_overlay(roads_layer) |>
  plot_map()