#30DayMapChallenge 2021 - Day 23 - GHSL

ghsl_nz.R

# Create a 3d render of NZ built up area density
# For #30DayMapChallenge 2021 - Day 23 - Global Human Settlement Layer
#  -- David Friggens, 2 December 2021

# This is a slight alteration of David Schoch's code:
# https://github.com/schochastics/30DayMapChallenge/blob/main/23_ghsl/23_ghsl.R

library(raster)
library(rayshader)

CRS_MERC <- "+proj=longlat +datum=WGS84 +no_defs"
CRS_MOLL <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# A world raster in Mollweide. Great for Europe, not so much for Aotearoa
world <- 
  raster("Day_23/GHS_BUILT_LDS2014_GLOBE_R2018A_54009_1K_V2_0.tif")

# We want to crop NZ in this box, which sounds good...
nz_bb <-
  as(extent(166.42, 178.58,
            -47.29, -34.39),
     'SpatialPolygons')
crs(nz_bb) <- CRS_MERC

# But we have to crop in Mollweide, so we end up picking up some of Australia
nz_bb_moll <-
  projectExtent(nz_bb,
                CRS_MOLL)
nz_moll <-
  crop(world, nz_bb_moll)

# So lets convert back to Mercator and crop again
nz <- 
  nz_moll %>%
  projectRaster(crs = CRS_MERC) %>% 
  crop(nz_bb)

nz_matrix <-
  raster_to_matrix(nz)

nz_matrix[is.na(nz_matrix)] <- 0


nz_matrix %>%
  sphere_shade(texture = create_texture("#000000", "#000000", "#000000", "#000000", "#ffffff")) %>%
  plot_3d(nz_matrix,
          zscale = 5,
          fov = 0, 
          theta = -80,
          zoom = 0.65,
          phi = 30,
          soliddepth = -10, solidcolor = "white", shadow = TRUE, shadowdepth = -12,
          shadowcolor = "black", background = "white",
          windowsize = c(1600, 1200)
  )

render_snapshot("Day_23/Day_23_ghsl_nz.png",
                software_render = FALSE,
                vignette = FALSE)