urbanisation_mountains.R

# Prepare world data

# First up, we need to load the built-up area data that we’re going to be plotting. We download this from the European Commission’s Global Human Settlement Data portal [https://ghsl.jrc.ec.europa.eu/datasets.php] — specifically using the links from this page [http://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_BUILT_LDSMT_GLOBE_R2015B/]. We want the 250m-resolution rasters for 1975 and 2015 (GHS_BUILT_LDS1975_GLOBE_R2016A_54009_250 and GHS_BUILT_LDS2014_GLOBE_R2016A_54009_250).

# Once you’ve downloaded these (they’re BIG, so might take a little while...), we can save ourselves a lot of hassle later on by re-projecting them into the same co-ordinate space as the other data we’re going to be using. Specifically we want to change their units from metres to lat/lon. We do this by:
# 1) Unzipping the archive, and then
# 2) Running the following script on the command-line:

# gdalwarp -t_srs EPSG:4326 -tr 0.01 0.01 path/to/your/built-up-area.tif path/to/your/built-up-area_reprojected.tif

# Make a note of the path to your re-projected .tif file, as we’ll be using this later on.

# Load the country ID raster: this will make it easier to clip to borders and coastlines later on. The GeoTiff was downloaded from here: http://sedac.ciesin.columbia.edu/data/set/gpw-v4-national-identifier-grid-rev10
raster.countries <- raster("~/Downloads/gpw-v4-national-identifier-grid-rev10_30_sec_tif/res.tif")
# Load country boundary polygons from Natural Earth (also useful for cropping later on)
world <- st_read("world/ne_10m_admin_0_countries.shp")

# Load the metadata for the country ID raster (so we know which unique numeric raster value corresponds to which country)
countries_dbf <- foreign::read.dbf("~/Downloads/gpw-v4-national-identifier-grid-rev10_30_sec_tif/gpw_v4_national_identifier_grid_rev10_30_sec.tif.vat.dbf")

# We’re going to be plotting data on urbanisation in China, so let’s get:
# 1) The unique raster value for China from the country ID raster
gridcode <- countries_dbf %>% filter(ISOCODE == "CHN") %>% pull(GRIDCODE)
# 2) China’s borders from the Natural Earth polygons
border <- world %>% dplyr::filter(ADM0_A3 == "CHN")

# Now we can crop the country ID raster to only our country of interest.
# First, we crop the huge country ID raster to the bounding box of China
countries_crop <- crop(raster.countries, border, snap="near")
# Next, we use a little trick to turn all raster cells outside of China into empty cells: we multiply our cropped raster by whether or not any given cell has the numeric value associated with China. Conveniently, this crops out water bodies as well
countries_crop <- countries_crop * (countries_crop == gridcode)
# Finally, we divide this raster by itself, giving us a binary grid of 1 for China and NA for anything else
countries_crop <- countries_crop/countries_crop

# Now we get the coordinates of our specific point of interest: the city of Shenzhen. We’ll use this to get a closer crop of our raster, and for positioning a text label later on. We’re using the ggmap package here, which requires a Google Maps Platform API key. Happily, this is completely free of charge. You can get one by going to this page [https://cloud.google.com/maps-platform/?__utma=102347093.636004825.1540831213.1547140121.1547140121.1&__utmb=102347093.0.10.1547140121&__utmc=102347093&__utmx=-&__utmz=102347093.1547140121.1.1.utmcsr=(direct)%7Cutmccn=(direct)%7Cutmcmd=(none)&__utmv=-&__utmk=123981194#get-started], choosing `places` and continuing from there
register_google(key = your_API_key_here)
coords <- geocode("Shenzhen, China")

# Now we have all of our rwaw materials, and we can focus on the steps required to actually render our 3D plot. First, let’s create a function that will take us from our raw built-up area raster to a cropped and prepared raster that can be read as an elevation map.
prepare_raster <- function(filepath){
  # First we load the global built-up areas raster
  raster_master <- raster(filepath)
  # Then we crop it to the Chinese border
  raster_crop <- crop(raster_master, extent(border), snap="near")
  # Then we multiply it by our binary 1/NA country ID raster. By now we’ve got rid of every cell on the earth’s surface that isn’t on land in China
  combined <- countries_crop * raster_crop
  # Not we use the wonderful `crop_raster_square` function from Neil Charles’ geoviz package to crop a nice 100km x 100km square around our point of interest — the city of Shenzhen
  crop <- crop_raster_square(combined, coords$lat, coords$lon, square_km=100)
  # This next stage is not essential, but I find it helpful. I’m multiplying the built-up area values by 1000, which gets them into the kind of value space we might usually associate with land heights.
  crop <- (crop * 1000)
  # This next step is a little trick that makes sure `rayshader` doesn’t think "flat" land (in our case — land that is not built-up) is treated as water. We’re telling R to take any raster cell that is on land, but has a built-up value of less than 1 (less than 0.001 in our original units) and assign to it a random value between 0 and 1. This is essentially creating a layer of very fine elevation "noise" over our map, so large rural areas are no longer "flat" in our pseudo elevation map, and thus will not be treated as water
  crop[crop < 1] <- runif(sum(crop@data@values<1, na.rm=T))
  # Finally, we make all NA values (cells that are not Chinese land) negative one, i.e below sea level. This will allow `rayshader` to easily render water later on. NB: if our area of interest covered more than one country, we might not want to treat land on secondary countries as land below sea-level like this, but since Shenzhen is on the coast, and no other countries are inside our crop, we’re safe to do this
  crop[is.na(crop)] <- -1
  crop
}

# We can now use the function we’ve just written to create pseudo-elevation rasters of built-up areas around Shenzhen as they were in 1975 ... (NB: the filepath you pass to this function will be the path to wherever you saved the reprojected version of your Built-Up Grid GeoTiff)
crop_1975 <- prepare_raster("~/Documents/geoData/DEMS/global/GHSL/built-up/GHS_BUILT_LDS1975_GLOBE_R2016A_54009_250_v1_0/GHS_BUILT_LDS1975_GLOBE_R2016A_54009_250_v1_0_reprojected.tif")
# ... and 2015
crop_2015 <- prepare_raster("~/Documents/geoData/DEMS/global/GHSL/built-up/GHS_BUILT_LDS2014_GLOBE_R2016A_54009_250_v1_0/GHS_BUILT_LDS2014_GLOBE_R2016A_54009_250_v1_0_reprojected.tif")
# Let’s plot those as flat "elevation" maps, just to get a rough idea of what we’re working with
plot(crop_1975)
plot(crop_2015)

# Next, we’re going to start the fun bit: converting our flat heightmap into 3D shapes!

# First up, we’re going to want to map our height values onto a colour scale. For a striking "high is hot" effect, I’m using a colour ramp that works with the FT’s colour palette, and ramps up from our famous pale pink to a dark red.

# This is where we come to another of our little tricks: we need to make sure our elevation palettes are consistent across our two graphics (1975 and 2015), i.e if central Shenzhen is 100% built-up today, but was only 75% built-up in 1975, we want to make sure that the peak today is a darker red than the peak in 1975 — 75% as dark, to be precise. To do this, we need to build temporary copies of the rasters we just plotted, where we manually ensure that they each have the same maximum value in order to render their colour-maps (don’t worry — we’ll go back to the real data when we render the graphics).

# First up, we store the min and max values across both rasters:
all_min <- min(c(crop_1975@data@values, crop_2015@data@values))
all_max <- max(c(crop_1975@data@values, crop_2015@data@values))

# Then we create temporary copies of each first raster, where we manually assign the maximum value (the most built-up raster cell) to the top-left cell of the raster
crop_1975_temp <- crop_1975
crop_2015_temp <- crop_2015
crop_1975_temp[1] <- all_max
crop_2015_temp[1] <- all_max

# That’s done — we’ll come back to those shortly.

# We now have everything we need to prepare our 3D scenes, so let’s write a function that will take care of that pipeline. We’re going to pass this function two rasters — the true rasters that we created in our first function, and our temporary manually-adjusted rasters that ensure a common colour ramp.
prepare_scene <- function(prepped_raster, scaled_elevation){
  # First up, we want to create a Digital Elevation Model from our true raster data
  DEM = matrix(raster::extract(prepped_raster,raster::extent(prepped_raster)),
                    nrow=ncol(prepped_raster),ncol=nrow(prepped_raster))
  
  # Next, we create an elevation overlay — this is our elevation (or in our case urbanisation) -based colour ramp.
  # We use our manually-adjusted raster data here. First we run it through our pink -> red colour palette. The top-left cell of this will be the darkest colour across our two palettes, ensuring that all other cell colours share the same linear colour ramp ...
  elevation_overlay <- elevation_shade(scaled_elevation, elevation_palette = c("#fff1e5", "orange", "darkorange", "firebrick"))
  # Then we get rid of any extraneous data bands (our final elevation overlay should have four bands — red, green, blue and alpha) ...
  elevation_overlay <- elevation_overlay[,,-5]
  # An then we assign a uniform value to the alpha channel, determining how opaque this will be in our final image
  elevation_overlay[,,4] <- 0.8
  
  # Now it’s time to use our colour ramp trick:
  # First, we create a second elevation layer, using our true raster. The top-left cell of this will be the true colour that it should have been before we made it the maximum value
  elevation_overlay_temp <- elevation_shade(prepped_raster, elevation_palette = c("#fff1e5", "orange", "darkorange", "firebrick"))
  
  # Next, we overwrite our false top-left cell value with the true one from our second elevation layer. This ensures the best of both worlds: all of our colours are correctly scaled with a common minimum and maximum, and all also show their true colour
  elevation_overlay[,,1][1] <- elevation_overlay_temp[,,1][1]
  elevation_overlay[,,2][1] <- elevation_overlay_temp[,,2][1]
  elevation_overlay[,,3][1] <- elevation_overlay_temp[,,3][1]
  
  # Now we can use the power of `rayshader` to calculate global (ray_shade) and ambient (ambient_shade) shadow maps. The `zscalez` and `sunangle` arguments set the vertical exaggeration of our "peaks", and the direction from which the sun is shining as it creates our shadows
  raymat = ray_shade(DEM, zscale=20, sunangle=135)
  ambmat = ambient_shade(DEM, zscale=20)
  
  # Finally, we take all of the above and prepare our 3D render scene. We’re adding an initial `sphere_shade` to render our 3D hillshade surface (using a custom colour texture here)
  scene <- DEM %>%
    sphere_shade(texture = create_texture("#fff1e5","#fff1e5","#fff1e5","#fff1e5","#fff1e5"), sunangle=135) %>%
    add_overlay(elevation_overlay) %>%
    # add_water(detect_water(DEM, zscale=1, cutoff=0.999), color="#F3FAFF") %>%
    add_shadow(raymat) %>%
    add_shadow(ambmat)
  
  return(list(scene=scene, DEM=DEM))
}

# Now we can use our new function to prepare scenes for our 1975 and 2015 graphics:
scene_1975 <- prepare_scene(crop_1975, crop_1975_temp)
scene_2015 <- prepare_scene(crop_2015, crop_2015_temp)

# We grab the co-ordinates for a centrally-positioned label over Shenzhen:
coords_cell <- as.numeric(rowColFromCell(crop_1975, cellFromXY(crop_1975, coords)))

# And, we plot our beautiful 3D urbanisation mountains!

# First up, 1975.
# Step one: we render the plot to a 3D graphics device
rayshader::plot_3d(
  scene_1975$scene,
  scene_1975$DEM,
  zscale = 30,
  solid = FALSE,
  shadow = TRUE,
  shadowdepth = -10,
  phi = 30,
  water = TRUE,
  waterdepth = 0,
  zoom = 0.7,
  windowsize = c(1000,800)
)
# Then we add our label:
render_label(scene_1975$DEM,x=coords_cell[1],y=coords_cell[2], z=1000,zscale=10, text = "Shenzhen, 1975",textsize = 2,linewidth = 0, alpha = 0)
# Then we render the current snapshot of our 3D plot to the R plotting device (this will let us export a flat image)
render_snapshot()
# Then we clear our 3D rendering device ...
rgl::rgl.clear()
# And close it...
rgl::rgl.close()
# And the two lines below will export a flat PNG image of the 3D plot
dev.copy(png, "shenzhen_1975.png", width=1000, height=800, res=144)
dev.off()

# And now the same for 2015.
rayshader::plot_3d(
  scene_2015$scene,
  scene_2015$DEM,
  zscale = 30,
  solid = FALSE,
  shadow = TRUE,
  shadowdepth = -10,
  phi = 30,
  water = TRUE,
  waterdepth = 0,
  zoom = 0.7,
  windowsize = c(1000,800)
)
render_label(scene_1975$DEM,x=coords_cell[1],y=coords_cell[2], z=1000,zscale=10, text = "Shenzhen, 2015",textsize = 2,linewidth = 0, alpha = 0)
render_snapshot()
rgl::rgl.clear()
rgl::rgl.close()
dev.copy(png, "shenzhen_2015.png", width=1000, height=800, res=144)
dev.off()

# Voila!