Penn MUSA Masterclass: 3D Mapping and Visualization with R and Rayshader
Tyler Morgan-Wall (@tylermorganwall), Institute for Defense Analyses
Personal website: https://www.tylermw.com
Rayshader website: https://www.rayshader.com
Github: https://www.github.com/tylermorganwall
First, let’s install all the required packages for this masterclass. I would recommend R 3.6.0 as a minimum requirement for this class–lower versions of R have a different RDS version, and if you have trouble with some of the raster operations you won’t be able to load the back-up RDS files.
# install.packages(c("ggplot2","raster", "rayrender", "spatstat", "spatstat.utils","suncalc","here", "sp","lubridate","rgdal", "magick", "av","xml2", "dplyr"))
#macOS: xcrun issue? type "xcode-select --install" into terminal
#macOS imager.so or libx11 errors? Install X11: https://www.xquartz.org
# install.packages("remotes")
# remotes::install_github("giswqs/whiteboxR")
# whitebox::wbt_init()
# remotes::install_github("tylermorganwall/rayshader")Now load the packages:
options(rgl.useNULL = FALSE)
library(ggplot2)
library(whitebox)
library(rayshader)
library(rayrender)
library(raster)
library(spatstat)
library(spatstat.utils)
library(suncalc)
library(sp)
library(lubridate)
library(rgdal)
setwd(here::here())Hillshading with Rayshader
We’re going to start by making our own 2D maps with rayshader. Here, we load in a DEM of the River Derwent in Hobart, Tasmania. We’re going to plot it using a few different methods, and then we’re going to go into the third dimension.
Let’s take a raster object and extract a bare R matrix to work with
rayshader. To convert the data from the raster format to a bare matrix,
we will use the rayshader function raster_to_matrix().
loadzip = tempfile()
download.file("https://dl.dropboxusercontent.com/s/8ltz4j599z4njay/dem_01.tif.zip", loadzip)
## Alternate Link
#download.file("https://tylermw.com/data/dem_01.tif.zip", loadzip)
hobart_tif = raster::raster(unzip(loadzip, "dem_01.tif"))
unlink(loadzip)
hobart_mat = raster_to_matrix(hobart_tif)
unlink("dem_01.tif")
hobart_mat[1:10,1:10]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 750 749 748 745 743 743 743 743 744 744
## [2,] 751 751 751 749 749 748 747 748 749 749
## [3,] 750 753 754 753 752 751 751 752 753 754
## [4,] 749 753 756 757 756 755 757 758 758 760
## [5,] 749 754 757 759 760 760 761 763 765 767
## [6,] 755 760 761 765 769 769 767 768 770 773
## [7,] 762 767 768 771 772 773 772 772 775 777
## [8,] 768 771 771 772 774 774 775 776 778 781
## [9,] 769 771 771 771 773 775 776 777 779 782
## [10,] 767 770 770 770 771 773 775 776 779 782
The structure of this data is simply an array of numbers. We will first
visualize the data with the height_shade() function, which performs a
traditional elevation-to-color mapping.
hobart_mat %>%
height_shade() %>%
plot_map()
Here, the default uses what R calls terrain.colors(), although it
doesn’t look like any terrain I’ve ever seen.
Where in nature does color map to elevation? Usually, you get this kind of mapping on large scale topographic features, such as mountains or canyons. Here’s a nice illustration from the 1848 book by Alexander Keith Johnston showing what I mean:
On smaller scale (more human-sized) features like hills or dunes, the landscape is dominated instead by ambient lighting and how the sun hits the surface. Here’s an example from Great Sand Dunes National Park, in Colorado. Note the (sometimes drastic) change in surface color that changes depending on the time of day or direction of the light:
Rather than color by elevation, let’s try coloring the surface by the
direction the slope is facing and the steepness of that slope. This is
implemented in the rayshader function, sphere_shade(). Here’s a video
explanation of how it works:
Taking data and pluging it into the sphere_shade() function, here’s
what we get:
hobart_mat %>%
sphere_shade() %>%
plot_map()
We note right away that it’s much more apparent in this hillshade that
there is a body of water weaving between the two mountains. To add water
to the scene we’ll call two functions in rayshader, detect_water() and
add_water(), that allow you to detect and add bodies of water directly
to the map using only the elevation data. detect_water() works by
looking for large, relatively flat, connected regions in a DEM. It’s
often the case (although not always) that the areas representing water
are far more flat than anything else in the matrix, so we take advantage
of that fact to extract and color in the water areas.
Occasionally these functions can result in false positives, but there
are tuning parameters to help. detect_water() provides the parameters
min_area, max_height, and cutoff in detect_water() so the user
can try to minimize false areas reported as water. min_area specifies
the minimum number of connected flat points that constitute a body of
water. max_height sets the maximum height that a region can be
considered to be water. cutoff determines how vertical a surface as to
be to be defined as flat: cutoff = 1.0 only accepts areas pointing
straight up as water, while cutoff = 0.99 allows for areas that are
slightly less than flat to be classified as water. The default is
cutoff = 0.999.
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat)) %>%
plot_map()
The default min_area is 1/400th the area of the matrix. There’s
nothing special about this number, but it’s a default that tended to
work fairly well for several datasets I tested it on–adjust it to suit
your dataset. We will deal with elevation data later in the class where
the defaults for this function don’t work.
You can also start to see why rayshader was built using the pipe
operator %>%. Here’s what the previous code would look like without
the
pipe:
plot_map(add_water(sphere_shade(hobart_mat),detect_water(hobart_mat)))
The code is far more readable with the pipe, and the step-by-step process of adding layers to our map is much more easy to see.
Both sphere_shade() and add_water() come with several built-in
palettes, and you can create your own with the create_texture()
function. Here’s the built-in textures (top is highlight color, and all
can be specified):
And here’s what a sampling of those textures look like (check the documentation for more):
hobart_mat %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(hobart_mat), color = "desert") %>%
plot_map()
hobart_mat %>%
sphere_shade(texture = "imhof4") %>%
add_water(detect_water(hobart_mat), color = "imhof4") %>%
plot_map()
hobart_mat %>%
sphere_shade(texture = "bw") %>%
add_water(detect_water(hobart_mat), color = "unicorn") %>%
plot_map()
Rayshader’s name comes one of the methods it uses to calculate
hillshades: raytracing, a method that realistically simulates how light
travels across the elevation model. Most traditional methods of
hillshading only use the local angle that the surface makes with the
light, and do not take into account areas that actually cast a shadow.
This basic type of hillshading is sometimes referred to as “lambertian”,
and is implemented in the function lamb_shade(). There’s one problem
with this mapping, though:
hobart_mat %>%
lamb_shade(zscale = 33) %>%
plot_map()
hobart_mat %>%
lamb_shade(zscale = 33, sunangle = 135) %>%
plot_map()
hobart_mat_inverted = -hobart_mat
hobart_mat %>%
lamb_shade(zscale = 33) %>%
plot_map()
# No difference!
hobart_mat_inverted %>%
lamb_shade(zscale = 33, sunangle = 135) %>%
plot_map()To shade surfaces using raytracing, rayshader draws rays originating
from each point towards a light source, specified using the sunangle
and sunaltitude argument. Here are two gifs showing how rayshader
calculates shadows with ray_shade:
Since our data is represented by grid points, rayshader uses a method called “bilinear interpolation” to determine if there are intersections with the ground even between points.
Let’s add a layer of shadows to this map, using the add_shadow() and
ray_shade() functions. We layer our shadows to our lamb_shade() base
layer.
hobart_mat %>%
lamb_shade(zscale = 33) %>%
add_shadow(ray_shade(hobart_mat, zscale = 33,
sunaltitude = 6, lambert = FALSE), 0.3) %>%
plot_map()
# No longer identical--we can tell canyons from mountains, thanks to raytracing
hobart_mat_inverted %>%
lamb_shade(zscale = 33, sunangle = 135) %>%
add_shadow(ray_shade(hobart_mat_inverted, zscale = 33, sunangle = 135,
sunaltitude = 6, lambert = FALSE), 0.3) %>%
plot_map()
We can combine all these together to make our final 2D map:
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat,zscale = 33, sunaltitude = 3,lambert = FALSE), max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat,zscale = 33,sunaltitude = 3), max_darken = 0.5) %>%
plot_map()
We can adjust the highlight/sun direction using the sunangle argument
in both the sphere_shade() function and the ray_shade() function.
We’ll start with the default angle: 315 degrees, or the light from the Northwest. One question you might ask yourself is: Why 315 degrees?
#Default angle: 315 degrees.
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat, zscale = 33, sunaltitude = 5,lambert = FALSE),
max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat,zscale = 33,sunaltitude = 5), max_darken = 0.8) %>%
plot_map()
#225 degrees
hobart_mat %>%
sphere_shade(sunangle = 225) %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat,sunangle = 225, zscale = 33, sunaltitude = 5,lambert = FALSE),
max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat,zscale = 33,sunaltitude = 5), max_darken = 0.8) %>%
plot_map()
You have complete freedom in choosing how your map is lit. And, as we
will show later, you can use the suncalc package to bring in the
actual position of the sun in the sky for a specific times and
locations.
We can also add the effect of ambient occlusion, which is sometimes referred to as the “sky view factor.” When light travels through the atmosphere, it scatters. This scattering turns the entire sky into a light source, so when less of the sky is visible (e.g. in a valley) it’s darker than when the entire sky is visible (e.g. on a mountain ridge).
Let’s calculate the ambient occlusion shadow layer for the Hobart data, and layer it onto the rest of the map.
hobart_mat %>%
ambient_shade() %>%
plot_map()
#No ambient occlusion
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat, zscale = 33, sunaltitude = 5,lambert = FALSE),
max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat,zscale = 33, sunaltitude = 5), max_darken = 0.7) %>%
plot_map()
#With ambient occlusion
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat, zscale = 33, sunaltitude = 5,lambert = FALSE),
max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat,zscale = 33, sunaltitude = 5), max_darken = 0.7) %>%
add_shadow(ambient_shade(hobart_mat), max_darken = 0.1) %>%
plot_map()
We can save these plots straight to disk, without displaying them, using
the save_png() function instead of plot_map().
hobart_mat %>%
sphere_shade() %>%
save_png("filename.png")
unlink("filename.png")3D Mapping with Rayshader
Now that we know how to perform basic hillshading, we can begin the real
fun part: making 3D maps. In rayshader, we do that simply by swapping
out plot_map() with plot_3d(), and adding the heightmap to the
function call. We don’t want to re-compute the ambient_shade() layer
every time we re-plot the landscape, so lets save it to a variable.
ambientshadows = ambient_shade(hobart_mat)
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat, sunaltitude = 3, zscale = 33, lambert = FALSE), max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat, sunaltitude = 3, zscale = 33), max_darken = 0.7) %>%
add_shadow(ambientshadows, max_darken = 0.1) %>%
plot_3d(hobart_mat, zscale = 10,windowsize = c(1000,1000))
render_snapshot()
rgl::rgl.clear()This opens up an rgl window that displays the 3D plot. Draw to
manipulate the plot, and control/ctrl drag to zoom in and out. To close
it, we can either close the window itself, or type in
rgl::rgl.close().
Just visualizing this on your screen is fun when exploring the data, but
we would also like to export our figure to an image file. If you want to
take a snapshot of the current view, rayshader provides the
render_snapshot() function, which is useful for rmarkdown documents
like this. If you use this without a filename, it will write and display
the plot to the current device. With a filename, it will write the image
to a PNG file in the local directory. For variety, let’s also change the
background/shadow color (arguments background and shadowcolor),
depth of rendered ground/shadow (arguments soliddepth and
shadowdepth), and add a title to the plot.
hobart_mat %>%
sphere_shade() %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat, sunaltitude = 3, zscale = 33, lambert = FALSE), max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat, sunaltitude = 3, zscale = 33), max_darken = 0.7) %>%
add_shadow(ambientshadows, max_darken = 0) %>%
plot_3d(hobart_mat, zscale = 10,windowsize = c(1000,1000),
phi = 40, theta = 135, zoom = 0.9,
background = "grey30", shadowcolor = "grey5",
soliddepth = -50, shadowdepth = -100)
render_snapshot(title_text = "River Derwent, Tasmania",
title_font = "Helvetica",
title_size = 50,
title_color = "grey90")
render_snapshot(filename = "derwent.png")
#Delete the file
unlink("derwent.png")If we want to programmatically change the camera, we can do so with the
render_camera() function. We can adjust four parameters: theta,
phi, zoom, and fov (field of view). Changing theta orbits the
camera around the scene, while changing phi changes the angle above
the horizon at which the camera is located. Here is a graphic
demonstrating this relation (and here’s a link to the video from the
presentation ):
zoom magnifies the current view (smaller numbers = larger
magnification), and fov changes the field of view. Higher fov values
correspond to a more pronounced perspective effect, while fov = 0
corresponds to a orthographic camera.
Here’s different views using the camera:
render_camera(theta = 90, phi = 30, zoom = 0.7, fov = 0)
render_snapshot()
render_camera(theta = 90, phi = 30, zoom = 0.7, fov = 90)
render_snapshot()
render_camera(theta = 120, phi = 20, zoom = 0.3, fov = 90)
render_snapshot()
We can also use some post-processing effects to help guide our viewer
through our visualization’s 3D space. The easiest method of directing
the viewer’s attention is directly labeling the areas we’d like them to
look at, using the render_label() function. This function takes the
indices of the matrix coordinate area of interest x and y, and
displays a text label at altitude z. We’ll also use the
title_bar_color argument to add a semi-transparent light bar behind
our title to help it stand out from the background.
#Windows? Turn off freetype.
if(.Platform$OS.type == "windows") {
freetype = FALSE
} else {
freetype = TRUE
}
render_label(hobart_mat, "River Derwent", textcolor ="white", linecolor = "white", freetype = freetype,
x = 450, y = 260, z = 1400, textsize = 2.5, linewidth = 4, zscale = 10)
render_snapshot(title_text = "render_label() demo, part 1",
title_bar_alpha = 0.8,
title_bar_color = "white")
render_label(hobart_mat, "Jordan River (not that one)", textcolor ="white", linecolor = "white", freetype = freetype,
x = 450, y = 140, z = 1400, textsize = 2.5, linewidth = 4, zscale = 10, dashed = TRUE)
render_snapshot(title_text = "render_label() demo, part 2",
title_bar_alpha = 0.8,
title_bar_color = "white")
We can replace all existing text with the clear_previous = TRUE, or
clear everything by calling render_label(clear_previous = TRUE) with
no other arguments.
render_camera(zoom = 0.9, phi = 50, theta = -45,fov = 0)
render_snapshot()
render_label(hobart_mat, "Mount Faulkner", textcolor ="white", linecolor = "white", freetype = freetype,
x = 135, y = 130, z = 2500, textsize = 2, linewidth = 3, zscale = 10, clear_previous = TRUE)
render_snapshot()
render_label(hobart_mat, "Mount Dromedary", textcolor ="white", linecolor = "white", freetype = freetype,
x = 320, y = 390, z = 1000, textsize = 2, linewidth = 3, zscale = 10)
render_snapshot()
render_label(clear_previous = TRUE, freetype = freetype)
render_snapshot()
rgl::rgl.close()This is not the only trick up our sleeve, though. Data visualization is not the first field to tackle the problem of “How do I direct a viewer’s attention in 3D space, projected on a 2D screen?” Cinematographers encounter the same issue when they’re filming (except when using small aperture sizes or lenses with short focal lengths). Below is a shot of Daniel Craig in the movie “Defiance”:
Even though the 3D space this shot is representing is fairly deep (with a distinct foreground, middle, and background), we know our attention is supposed to be focused on the man in the center. The use of a shallow depth of field tells our viewers what isn’t important, and where their attention should be.
We can apply this same technique in rayshader using the render_depth()
function. This applies a depth of field effect, as if our scene was
being photographed by a real camera. Let’s focus in on the river in the
background. The focus parameter should be between 0 and 1, where 1 is
the background and 0 is the foreground. The range of depths will depend
on our camera position–calling the function with preview_focus = TRUE
will print the range in the current view, and you can use the red line
to see exactly where the focal plane will be positioned.
hobart_mat %>%
sphere_shade(sunangle = 60) %>%
add_water(detect_water(hobart_mat), color = "lightblue") %>%
add_shadow(ray_shade(hobart_mat, sunangle = 60, sunaltitude = 3, zscale = 33, lambert = FALSE), max_darken = 0.5) %>%
add_shadow(lamb_shade(hobart_mat, sunangle = 60, sunaltitude = 3, zscale = 33), max_darken = 0.7) %>%
add_shadow(ambientshadows, max_darken = 0.1) %>%
plot_3d(hobart_mat, zscale = 10,windowsize = c(1000,1000),
background = "#edfffc", shadowcolor = "#273633")
render_camera(theta = 120, phi = 20, zoom = 0.3, fov = 90)
render_depth(focus = 0.81, preview_focus = TRUE)## [1] "Focal range: 0.61623-0.977782"
render_depth(focus = 0.9, preview_focus = TRUE)## [1] "Focal range: 0.61623-0.977782"
render_depth(focus = 0.81)
This effect is rather subtle, so let’s increase the focal length of the
camera using argument focallength. This will make for a shallower
depth of field (increase the blurring effect in areas far from the focal
plane). We can also add a title_* arguments like in
render_snapshot(), and we’ll add a lens vignetting effect (which
darkens the edges) by setting vignette = TRUE. This option is also
available in render_snapshot() and
render_movie().
render_depth(focus = 0.81, focallength = 200, title_bar_color = "black", vignette = TRUE,
title_text = "The River Derwent, Tasmania", title_color = "white", title_size = 50)
rgl::rgl.close()Rayshader can be used for more than just topographic data
visualization–we can also easily visualize the intersection between
land and water. Here, we’ll use the built-in montereybay dataset,
which includes both bathymetric and topographic data for Monterey Bay,
California. We include a transparent water layer by setting water = TRUE in plot_3d(), and add lines showing the edges of the water by
setting the waterlinecolor argument.
montereybay %>%
sphere_shade() %>%
plot_3d(montereybay, water = TRUE, waterlinecolor = "white",
theta = -45, zoom = 0.9, windowsize = c(1000,1000),zscale = 50)
render_snapshot(title_text = "Monterey Bay, California",
title_color = "white", title_bar_color = "black")
By default, plot_3d() sets the water level at 0 (sea level), but we
can change this, either by adjusting waterdepth in our call to
plot_3d(), or by calling the render_water() function after the fact.
Here, we make the water slightly less transparent with the wateralpha
argument (wateralpha = 1 is opaque, wateralpha = 0 is transparent).
render_water(montereybay, zscale = 50, waterdepth = -100,
waterlinecolor = "white", wateralpha = 0.7)
render_snapshot(title_text = "Monterey Bay, California (water level: -100 meters)",
title_color = "white", title_bar_color = "black")
render_water(montereybay, zscale = 50, waterdepth = 30,
waterlinecolor = "white", wateralpha = 0.7)
render_snapshot(title_text = "Monterey Bay, California (water level: 30 meters)",
title_color = "white", title_bar_color = "black")
rgl::rgl.close()Although our underlying matrix data is rectangular, we aren’t limited to
rectangular 3D plots. Any entry in the matrix that has a value of NA
will be sliced out of the visualization. Here, we slice out just the
undersea data from the montereybay dataset. We still use the full
matrix to calculate the hillshade and shadows, but pass the sliced
elevation matrix to plot_3d():
mont_bathy = montereybay
mont_bathy[mont_bathy >= 0] = NA
montereybay %>%
sphere_shade() %>%
add_shadow(ray_shade(mont_bathy,zscale = 50, sunaltitude = 15, lambert = FALSE),0.5) %>%
plot_3d(mont_bathy, water = TRUE, waterlinecolor = "white",
theta = -45, zoom = 0.9, windowsize = c(1000,1000))## `montereybay` dataset used with no zscale--setting `zscale=50`. For a realistic depiction, raise `zscale` to 200.
render_snapshot(title_text = "Monterey Bay Canyon",
title_color = "white",
title_bar_color = "black")
rgl::rgl.clear()To slice out just the land, we’ll just swap the comparison:
mont_topo = montereybay
mont_topo[mont_topo < 0] = NA
montereybay %>%
sphere_shade() %>%
add_shadow(ray_shade(mont_topo, zscale = 50, sunaltitude = 15, lambert = FALSE),0.5) %>%
plot_3d(mont_topo, shadowdepth = -50,
theta = 135, zoom = 0.9, windowsize = c(1000,1000))## `montereybay` dataset used with no zscale--setting `zscale=50`. For a realistic depiction, raise `zscale` to 200.
render_snapshot(title_text = "Monterey Bay (sans water)",
title_color = "white",
title_bar_color = "black")
rgl::rgl.clear()Alternatively, if we just want a more interesting base shape and we
don’t mind losing some data, we can have rayshader carve your dataset
into either a hexagon or a circle by specifying the baseshape argument
in plot_3d() (I also throw in some new background colors, shadow
colors, and the vignette effect for variety):
montereybay %>%
sphere_shade() %>%
add_shadow(ray_shade(montereybay,zscale = 50,sunaltitude = 15,lambert = FALSE),0.5) %>%
plot_3d(montereybay, water = TRUE, waterlinecolor = "white", baseshape = "hex",
theta = -45, zoom = 0.7, windowsize = c(1000,1000),
shadowcolor = "#4e3b54", background = "#f7e8fc")## `montereybay` dataset used with no zscale--setting `zscale=50`. For a realistic depiction, raise `zscale` to 200.
render_snapshot(title_text = "Monterey Bay Canyon, Hexagon", vignette = TRUE,
title_color = "white", title_bar_color = "black", clear = TRUE)
montereybay %>%
sphere_shade() %>%
add_shadow(ray_shade(montereybay,zscale = 50,sunaltitude = 15,lambert = FALSE),0.5) %>%
plot_3d(montereybay, water = TRUE, waterlinecolor = "white", baseshape = "circle",
theta = -45, zoom = 0.7, windowsize = c(1000,1000),
shadowcolor = "#4f3f3a", background = "#ffeae3")## `montereybay` dataset used with no zscale--setting `zscale=50`. For a realistic depiction, raise `zscale` to 200.
render_snapshot(title_text = "Monterey Bay Canyon, Circle",
title_color = "white", title_bar_color = "black", clear = TRUE)
Taking snapshots of our map is useful, but a static visualization isn’t
the ideal method to represent our 3D data. For the viewer, the depth
variable has to be inferred from context, which can be ambiguous and
easily misinterpreted. A far better form for our visualization is a
movie/animation–we can guide our viewers on a 3D tour of our data, which
can resolve many visual ambiguities. Rayshader makes this easy with the
render_movie() function, which creates a movie (mp4 file) which can
easily be shared on social media. You can add titles and all of the same
arguments you can to render_snapshot().
By default, render_movie() has two animations built in: a basic orbit
type = "orbit", an oscillating animation type = "oscillate". The
frames argument is the number of frames, while the fps is the frames
per second. The total length of the video is frames/fps. Generally, you
shouldn’t change from either 30 or 60 frames per second if you want to
share your videos on social media. You also shouldn’t make the camera
move too fast or spin too quickly–fast camera movements can cause nausea
and aren’t very interpretable.
montereybay %>%
sphere_shade() %>%
plot_3d(montereybay, water = TRUE, waterlinecolor = "white",
theta = -45, zoom = 0.9, windowsize = c(600,600))## `montereybay` dataset used with no zscale--setting `zscale=50`. For a realistic depiction, raise `zscale` to 200.
#Orbit will start with current setting of phi and theta
render_movie(filename = "montbay.mp4", title_text = 'render_movie(type = "orbit")',
phi = 30 , theta = -45)## [1] "montbay.mp4"
render_movie(filename = "montbayosc.mp4", phi = 30 , theta = -90, type = "oscillate",
title_text = 'render_movie(type = "oscillate")', title_color = "black")## [1] "montbayosc.mp4"
unlink("montbay.mp4")
unlink("montbayosc.mp4")You can also set type = "custom" and pass in a vector of camera values
to the phi, theta, zoom, and fov arguments, and each value will
be treated as a single frame in our animation. Here I’ve also provided
an “easing” function that smoothly transitions between two values, so we
can zoom in and out without jarring
changes.
ease_function = function(beginning, end, steepness = 1, length.out = 180) {
single = (end) + (beginning - end) * 1/(1 + exp(seq(-10, 10, length.out = length.out)/(1/steepness)))
single
}
zoom_values = c(ease_function(1,0.3), ease_function(0.3,1))
#This gives us a zoom that looks like this:
ggplot(data.frame(x = 1:360,y = zoom_values),) +
geom_line(aes(x = x,y = y),color = "red",size = 2) +
ggtitle("Zoom value by frame")
render_movie(filename = "montbaycustom.mp4", type = "custom",
phi = 30 + 15 * sin(1:360 * pi /180),
theta = -45 - 1:360,
zoom = zoom_values)## [1] "montbaycustom.mp4"
rgl::rgl.close()
unlink("montbaycustom.mp4")If you want to do something more advanced, like vary the water depth
between each frame or call render_depth() to create a depth of field
effect, you’ll need to generate the movie manually in a for loop.
Luckily, this is fairly simple using the av package. Let’s vary the
water level from the most recent ice age minimum of -130 meters to the
current sea level. We use the av package to combine them into a movie.
# montereybay %>%
# sphere_shade(texture = "desert") %>%
# plot_3d(montereybay, windowsize = c(600,600),
# shadowcolor = "#222244", background = "lightblue")
#
# render_camera(theta = -90, fov = 70,phi = 30,zoom = 0.8)
#
# for(i in 1:60) {
# render_water(montereybay, zscale = 50, waterdepth = -60 - 60 * cos(i*pi*6/180),
# watercolor = "#3333bb",waterlinecolor = "white", waterlinealpha = 0.5)
# render_snapshot(filename = glue::glue("iceage{i}.png"), title_size = 30, instant_capture = TRUE,
# title_text = glue::glue("Sea level: {round(-60 - 60 *cos(i*pi*6/180),1)} meters"))
# }
#
# av::av_encode_video(glue::glue("iceage{1:60}.png"), output = "custom_movie.mp4",
# framerate = 30)
#
# rgl::rgl.close()
#
# unlink(glue::glue("iceage{1:60}.png"))
# unlink("custom_movie.mp4")Finally, I’m going to show you a really cool feature I’ve just recently
implemented: render_highquality(). This calls a powerful rendering
engine built-in to rayshader to create truly stunning 3D visualizations.
You can control the view entirely through the rgl
window–render_highquality() recreates the scene using the
rayrender package and returns a beautiful, pathtraced rendering with
realistic shadows and lighting. You can adjust the light(s) angle,
direction, color, and intensity, as well as add additional lights and
objects to the scene.
hobart_mat %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(hobart_mat), color = "desert") %>%
plot_3d(hobart_mat, zscale = 10)
render_highquality()
Notice we no longer need any of the shadow generating functions–the
shadows are calculated automatically as part of the rendering process.
You can adjust the light direction, or even add your own lights using
rayrender–it’s an extremely powerful tool to create visualizations.
You can decrease the noise by setting the samples argument to a higher
number in render_highquality() (default is 100).
Adding a custom light:
# `sphere()` and `diffuse()` are rayrender functions
render_highquality(light = FALSE, scene_elements = sphere(y = 150, radius = 30,material = diffuse(lightintensity = 40, implicit_sample = TRUE)))
rgl::rgl.close()Visualizing a 3D model with water:
montereybay %>%
sphere_shade() %>%
plot_3d(montereybay, zscale = 50, water = TRUE)
render_camera(theta = -45, zoom = 0.7, phi = 30,fov = 70)
render_highquality(lightdirection = 100, lightaltitude = 45, lightintensity = 800,
clamp_value = 10, title_text = "Monterey Bay, CA",
title_color = "white", title_bar_color = "black")
rgl::rgl.close()Elevation Data Sources
Here are some sources for downloading elevation/bathymetry data. I have included several how-to PDF guides as well.
Philadelphia Guide
https://github.com/tylermorganwall/MusaMasterclass/raw/master/data_source_guides/philly.pdf Alternate: https://dl.dropboxusercontent.com/s/5ls2gfjssosr2pa/philly.pdf?dl=0
Source: https://maps.psiee.psu.edu/preview/map.ashx?layer=2021
Florida Guide
https://github.com/tylermorganwall/MusaMasterclass/raw/master/data_source_guides/florida.pdf Alternate: https://dl.dropboxusercontent.com/s/adf1dsmytnw9f8c/florida.pdf
Source: http://dpanther2.ad.fiu.edu/Lidar/lidarNew.php
Texas Guide
https://github.com/tylermorganwall/MusaMasterclass/raw/master/data_source_guides/texas.pdf Alternate: https://dl.dropboxusercontent.com/s/89xs8yaet8s0x35/texas.pdf
Source: https://tnris.org/stratmap/elevation-lidar/
SRTM/OpenTopography Global Elevation Data Guide
https://github.com/tylermorganwall/MusaMasterclass/raw/master/data_source_guides/srtm_global.pdf Alternate: https://dl.dropboxusercontent.com/s/zu46o6ilzejy6gx/srtm_global.pdf?dl=0
Source: http://opentopo.sdsc.edu/raster?opentopoID=OTSRTM.082015.4326.1
GEBCO Global Bathymetry Data Guide
https://github.com/tylermorganwall/MusaMasterclass/raw/master/data_source_guides/gebco_global_bathymetry.pdf Alternate: https://dl.dropboxusercontent.com/s/i0regel891xajz3/gebco_global_bathymetry.pdf
Source: https://download.gebco.net
Mapping the Effects of Rising Sea Levels
How are coastal communities impacted by rising sea levels? That’s one of the most pressing problems presented by our changing climate. Actually showing how communities are affected by rising sea levels by visualizing the water level directly is, in my opinion, far more impactful than just discussing the numbers. Here, we’re going to combine lidar data from Miami Beach and predictions about sea level rise to show how that community will be affected.
We’re not just going to take into account sea level rise–climate change results in warmer oceans, which in turn results in more powerful hurricanes. While gradual sea level rise might be countered by levees and seawalls, these hurricanes bring storm surges that can completely overwhelm and inundate coastal cities. And since Miami is directly in the predicted path of many hurricanes, it’s important we take that into account.
Here’s a link to a Google Maps view of the region we’re going to analyze:
We’re going to model the most extreme scenario: 10.7 ft sea level rise by the year 2100, with a 9 ft storm surge. Our data is in the Florida East State Plane coordinate system, so the elevation is in feet. I pulled the following two datasets from this website (follow the guide for Florida to download your own data):
http://dpanther2.ad.fiu.edu/Lidar/lidarNew.php
Now, let’s load some lidar data into R. We’ll do this using the
whitebox package, which has several functions for manipulating and
transforming lidar data, among many other useful features. Many lidar
datasets are far too large to easily visualize, and are often noisy–we
are going to use whitebox to clean the data, and turn it into a simple
elevation model that includes buildings. This can take a minute, so I’ve
commented out the below code–you can run it if you want.
## Download this lidar data
# download.file("https://dl.dropboxusercontent.com/s/hkdxt1zbsjp68jl/LID2007_118755_e.zip",destfile = "LID2007_118755_e.zip")
# download.file("https://dl.dropboxusercontent.com/s/omvb63urfby6ddb/LID2007_118754_e.zip",destfile = "LID2007_118754_e.zip")
##Alternate links
# download.file("https://www.tylermw.com/data/LID2007_118755_e.zip",destfile = "LID2007_118755_e.zip")
# download.file("https://www.tylermw.com/data/LID2007_118754_e.zip",destfile = "LID2007_118754_e.zip")
#
# unzip("LID2007_118755_e.zip")
# unzip("LID2007_118754_e.zip")
#
# whitebox::wbt_lidar_tin_gridding(here::here("LID2007_118755_e.las"),
# output = here::here("miamibeach.tif"),
# resolution = 1, verbose_mode = TRUE,
# exclude_cls = '3,4,5,7,13,14,15,16,18')
#
# whitebox::wbt_lidar_tin_gridding(here::here("LID2007_118754_e.las"),
# output = here::here("miamibeach2.tif"),
# resolution = 1, verbose_mode = TRUE,
# exclude_cls = '3,4,5,7,13,14,15,16,18')While this runs, let’s talk about general limitations you’ll encounter
when dealing with lidar data. If you notice, you’ll see whitebox
actually isn’t loading any data into your environment, and for good
reason. R is memory hungry, and processing a full lidar dataset would be
difficult on most standard machines. The whitebox package operates on
files outside of R, and outputs a processed file when it’s done. Once
the data is transformed from a collection of lidar point measurements to
a regular grid of elevation values, we can load it into R and get to
work.
We’re using whitebox to remove several categories of points in the
data–here’s a table of what these numbers (passed in argument
exclude_cls) are referring to:
| CLS | Classification |
|---|---|
| 0 | Never classified |
| 1 | Unassigned |
| 2 | Ground |
| 3 | Low Vegetation |
| 4 | Medium Vegetation |
| 5 | High Vegetation |
| 6 | Building |
| 7 | Low Point |
| 9 | Water |
| 10 | Rail |
| 11 | Road Surface |
| 13 | Wire - Guard (Shield) |
| 14 | Wire - Conductor (Phase) |
| 15 | Transmission Tower |
| 16 | Wire-Structure Connector (Insulator) |
| 17 | Bridge Deck |
| 18 | High Noise |
Let’s download the pre-processed tifs, if you aren’t generating them from the lidar data:
download.file("https://dl.dropboxusercontent.com/s/rwajxbdwtkcw50c/miamibeach.tif", destfile = "miamibeach.tif")
download.file("https://dl.dropboxusercontent.com/s/39abkh87h05n7rm/miamibeach2.tif", destfile = "miamibeach2.tif")
## Alternate Links
# download.file("https://www.tylermw.com/data/miamibeach.tif", destfile = "miamibeach.tif")
# download.file("https://www.tylermw.com/data/miamibeach2.tif", destfile = "miamibeach2.tif")We’ll load them in with the raster package, and then merge to the two together. Because they came from the same source, they have the same coordinate system, so the merge here happens seamlessly. We then convert the raster to a matrix.
miami1 = raster::raster("miamibeach.tif")
miami2 = raster::raster("miamibeach2.tif")
miami_combined = raster::merge(miami1, miami2)
miami_beach = raster_to_matrix(miami_combined)
dim(miami_beach)## [1] 10000 5000
This is very large! Too large to easily visualize. Let’s downsample the data by a factor of 4.
# 1/4th the size, so 0.25 for the second argument
miami_beach_small = reduce_matrix_size(miami_beach, 0.25)
## Backup:
#download.file("https://www.tylermw.com/data/miami_beach_small.Rds",destfile = "miami_beach_small.Rds")
#miami_beach_small = readRDS("miami_beach_small.Rds")
dim(miami_beach_small)## [1] 2500 1250
miami_beach_small %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(miami_beach_small,zscale = 4)) %>%
add_shadow(ray_shade(miami_beach_small, zscale = 4, multicore = TRUE,
sunaltitude = 10, sunangle = -110),0.3) %>%
plot_map()
Detect water didn’t work great here. This is because the lidar is actually returning the water surface–you can see ripples and waves where the water should be. It’s not perfectly flat, so we need to adjust the parameters to get it to work.
miami_beach_small %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(miami_beach_small, cutoff=0.2, zscale=4,
min_area = length(miami_beach_small)/150,
max_height = 3)) %>%
add_shadow(ray_shade(miami_beach_small, zscale = 4, multicore = TRUE,
sunaltitude = 10, sunangle = -110),0.3) %>%
plot_map()
We should slice this map down to only the region of interest–Miami Beach
proper. We can do this with the crop function in the raster package.
We first create an extent object, which specifies the boundaries of
our map. We then convert the boundaries from lat/long to the coordinate
system of our lidar data, which is the Florida East State Plane. I have
provided a function here, lat_long_to_other, that performs this
transformation. You can use it with other coordinate systems as well–you
will just need to provide the “EPSG” code of the desired coordinate
system.
To figure out the desired latitudes and longitude boundaries for Miami Beach, I went to Google Maps and double clicked on the points I wanted to be the bottom left and top right of our map. Google Maps will display the longitude and latitude at the bottom, which you can then just copy and paste.
lat_long_to_other = function(lat,long, epsg_code) {
data = data.frame(long=long, lat=lat)
coordinates(data) <- ~ long+lat
#Input--lat/long
proj4string(data) <- CRS("+init=epsg:4326")
#Convert to coordinate system specified by EPSG code
xy = data.frame(spTransform(data, CRS(paste0("+init=epsg:", epsg_code))))
colnames(xy) = c("x","y")
return(unlist(xy))
}
# Florida East State Plane EPSG code is 2236
# I grabbed these latitudes and longitudes from double clicking on google maps
bottomleft = lat_long_to_other(25.763675, -80.142499, 2236)
topright = lat_long_to_other(25.775562, -80.127569, 2236)
# Create an extent object:
e = extent(c(bottomleft[1], topright[1], bottomleft[2], topright[2]))
# Use that extent object to crop the original grid to our desired area
crop(miami_combined, e) %>%
raster_to_matrix() %>%
reduce_matrix_size(0.25) ->
miami_cropped
## Backup
#download.file("https://www.tylermw.com/data/miami_cropped.Rds",destfile = "miami_cropped.Rds")
#miami_cropped = readRDS("miami_cropped.Rds")
## R 3.5 and below:
#download.file("https://www.tylermw.com/data/miami_cropped_oldr.Rds",destfile = "miami_cropped_oldr.Rds")
#miami_cropped = readRDS("miami_cropped_oldr.Rds")
miami_cropped %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(miami_cropped, cutoff = 0.3,
min_area = length(miami_cropped)/20,
max_height = 2)) %>%
add_shadow(ray_shade(miami_cropped, zscale = 4, multicore = TRUE,
sunaltitude = 30, sunangle = -110),0.3) %>%
plot_map()
Cropped and ready. Now that we have our data prepped and reduced, let’s visualize it in 3D. The coordinate system here is defined in feet, so that’s what we’ll enter for the values of waterdepth.
miami_cropped %>%
sphere_shade(texture = "desert") %>%
add_shadow(ray_shade(miami_cropped, zscale = 4, multicore = TRUE),0.3) %>%
plot_3d(miami_cropped, zscale = 4, water = TRUE, waterdepth = 0,
zoom=0.85, windowsize = 1000,
background = "grey50", shadowcolor = "grey20")
render_camera(phi = 45,fov = 70,zoom = 0.45,theta = 25)
render_snapshot(title_text = "Modern Day Miami Beach, Sea Level: 0 feet",
title_bar_color = "black",
title_color = "white", vignette = 0.2)
render_camera(phi = 45,fov = 70,zoom = 0.44,theta = 20)
render_water(miami_cropped,zscale = 4,waterdepth = 3)
render_snapshot(title_text = "Miami Beach, Sea Level: 3 feet",
title_bar_color = "black",
title_color = "white", vignette = 0.2)
render_camera(phi = 45,fov = 70,zoom = 0.42,theta = 10)
render_water(miami_cropped,zscale = 4,waterdepth = 7)
render_snapshot(title_text = "Miami Beach, Sea Level: 7 feet (Max pred. 2100 sea level rise)",
title_bar_color = "black",
title_color = "white", vignette = 0.2)
render_camera(phi = 45,fov = 70,zoom = 0.41,theta = 5)
render_water(miami_cropped,zscale = 4,waterdepth = 16)
render_snapshot(title_text = "Miami Beach, Sea Level: 16 feet (Max sea level rise + 9 ft storm surge)",
title_size = 30,
title_bar_color = "black",
title_color = "white", vignette = 0.2)
rgl::rgl.close()Even with only a few feet of sea level rise, what we now consider to be extreme flooding events will become far more common, and many areas by the coast will become unlivable. We focused here on a relatively wealthy part of the nation–visualizing sea level rise in well-known areas brings out-sized attention to the problem–but this type of analysis is important to perform for communities all along the coast. I encourage you to grab elevation/lidar data for smaller coastal communities and see how rising sea levels could affect them.
Mapping the Shadows of Philadelphia
Now let’s use rayshader to perform an analysis using what it does best: calculating shadows! Back in 2015, the New York Times had a fantastic piece by Quoctrung Bui and Jeremy White on mapping the shadows of New York City.
https://www.nytimes.com/interactive/2016/12/21/upshot/Mapping-the-Shadows-of-New-York-City.html
They took three days in the year and mapped out how much each spot in the city spent in shadow: the winter solstice (December 21st), the summer solstice (June 21st), and the autumnal equinox (Sept 22nd). They worked with researchers at the Tandon school of Engineering at NYU to develop a raytracer to perform these calculations. We don’t have access to that tool, but we have lidar data and rayshader–that’s all we’ll need. This type of analysis is particularly important to assess the impact of so-called “pencil towers” that have gone up in Manhattan–extremely tall, skinny skyscrapers lining Central Park and catering to the super wealthy. How do tall buildings affect sunlight, a shared common good and scarce resource in densely packed urban areas?
We are going to perform and visualize a similar analysis for a hypothetical skyscraper in West Philadelphia, using lidar data from Pennsylvania.
Let’s start by loading some lidar data into R. Here’s a link to the PDF file with instructions on how to load data for Philly.
https://github.com/tylermorganwall/MusaMasterclass/raw/master/data_source_guides/philly.pdf
Let’s walk through downloading a specific dataset:
https://maps.psiee.psu.edu/preview/map.ashx?layer=2021
Here, we’re going to load our lidar dataset of Penn Park, right by the Schuylkill. We’re going to again use the whitebox package to prep the data.
#291.3 MB
# download.file("https://dl.dropboxusercontent.com/s/2jge03ptvsley62/26849E233974N.zip", destfile = "26849E233974N.zip")
## Alternate Link
# download.file("https://www.tylermw.com/data/26849E233974N.zip", destfile = "26849E233974N.zip")
unzip("26849E233974N.zip")
whitebox::wbt_lidar_tin_gridding(here::here("26849E233974N.las"),
output = here::here("phillydem.tif"), minz = 0,
resolution = 1, exclude_cls = '3,4,5,7,13,14,15,16,18')## [1] "lidar_tin_gridding - Elapsed Time (including I/O): 42.260s"
Now, let’s load the data into R.
# backups in case whitebox doesn't work for you
download.file("https://dl.dropboxusercontent.com/s/3auywgq93lokurf/phillydem.tif", destfile = "phillydem.tif")
## Alternate Link
# download.file("https://www.tylermw.com/data/phillydem.tif", destfile = "phillydem.tif")
phillyraster = raster::raster("phillydem.tif")
building_mat = raster_to_matrix(phillyraster)A 2640x2640 matrix is easily processed in R, but it’s not ideal for developing a visualization. Large matrices (greater than about 2000x2000, but depends on your machine) take a while to display in 3D. It’s faster to prototype your visualizations using a smaller dataset–once you are happy with the end result, you can re-run the script with the original hi-res data. We’re time limited and our visualization doesn’t require a high resolution dataset, so we won’t use one. Let’s reduce it by a factor of 2.
## Alternate Link
# download.file("https://www.tylermw.com/data/building_mat_small.Rds", destfile = "building_mat_small.Rds")
# building_mat_small = readRDS("building_mat_small.Rds")
building_mat_small = reduce_matrix_size(building_mat, 0.5)
dim(building_mat_small)## [1] 1320 1320
We’ll start our analysis by calculating the shadows from the existing
buildings in this area, so we can compare the relative change with the
addition of our skyscraper. We will be using the suncalc package to
calculate the position of the sun in the sky over Philadelphia.
suncalc::getSunlightPosition() takes a date, time, latitude, and
longitude and returns the angular position of the sun in the sky. You
can also use suncalc::getSunlightTimes() to obtain the times for
sunrise and sunset.
This calculation isn’t that slow (~3 minutes), but it’s slow enough that I’ve precomputed the results:
getSunlightTimes(as.Date("2019-06-21"), lat = 39.9526, lon = -75.1652,tz = "EST")## date lat lon solarNoon nadir
## 1 2019-06-21 39.9526 -75.1652 2019-06-21 12:03:39 2019-06-21 00:03:39
## sunrise sunset sunriseEnd
## 1 2019-06-21 04:33:22 2019-06-21 19:33:57 2019-06-21 04:36:38
## sunsetStart dawn dusk
## 1 2019-06-21 19:30:40 2019-06-21 04:00:31 2019-06-21 20:06:47
## nauticalDawn nauticalDusk nightEnd
## 1 2019-06-21 03:18:49 2019-06-21 20:48:29 2019-06-21 02:30:09
## night goldenHourEnd goldenHour
## 1 2019-06-21 21:37:09 2019-06-21 05:14:06 2019-06-21 18:53:13
#Start and hour after sunrise and end an hour before sunset
philly_time_start = ymd_hms("2019-06-21 05:30:00", tz = "EST")
philly_time_end= ymd_hms("2019-06-21 18:30:00", tz = "EST")
temptime = philly_time_start
philly_existing_shadows = list()
sunanglelist = list()
counter = 1
# while(temptime < philly_time_end) {
# sunangles[[counter]] = suncalc::getSunlightPosition(date = temptime, lat = 39.9526, lon = -75.1652)[4:5]*180/pi
# print(temptime)
# philly_existing_shadows[[counter]] = ray_shade(building_mat_small,
# sunangle = sunangles[[counter]]$azimuth+180,
# sunaltitude = sunangles[[counter]]$altitude,
# lambert = FALSE, zscale = 2,
# multicore = TRUE)
# temptime = temptime + duration("3600s")
# counter = counter + 1
# }
#Downloading the pre-computed data to save time
download.file("https://dl.dropboxusercontent.com/s/ipgcg51ct8esg3w/philly_existing_shadows.Rds",
destfile = "philly_existing_shadows.Rds")
# download.file("https://www.tylermw.com/data/philly_existing_shadows.Rds",
# destfile = "philly_existing_shadows.Rds")
philly_existing_shadows = readRDS("philly_existing_shadows.Rds")
str(philly_existing_shadows)## List of 13
## $ : num [1:1320, 1:1320] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
## $ : num [1:1320, 1:1320] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
## $ : num [1:1320, 1:1320] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
plot_map(philly_existing_shadows[[2]])
# Pithy "data science" way
shadow_coverage = Reduce(`+`, philly_existing_shadows)/length(philly_existing_shadows)
## Verbose "programmer" way
#shadow_coverage = matrix(0, nrow(philly_existing_shadows), ncol(philly_existing_shadows))
#for(i in seq_len(length(philly_existing_shadows)) {
# shadow_coverage = shadow_coverage + philly_existing_shadows[[i]]
#}
#shadow_coverage = shadow_coverage/length(philly_existing_shadows)
shadow_coverage %>%
reshape2::melt(varnames = c("x","y"), value.name = "light") %>%
ggplot() +
geom_raster(aes(x = x,y = y,fill = light)) +
scale_fill_viridis_c("Daily Light\nCoverage") +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
coord_fixed()
Now, let’s build our skyscraper in Penn Park. If we had an actual building we wanted to model, we would create a polygon using the spatial extent of the building and merge that with our raster. However, since this building is entirely imaginary, we’re going to take a shortcut and use a fun trick to draw it in our R graphics device, and merge that into our raster. Let’s load in some functions to help us do this:
owin2Polygons = function(x, id = "1") {
stopifnot(is.owin(x))
x = as.polygonal(x)
closering = function(df) { df[c(seq(nrow(df)), 1), ] }
pieces = lapply(x$bdry,
function(p) {
Polygon(coords = closering(cbind(p$x,p$y)),
hole = is.hole.xypolygon(p)) })
z = Polygons(pieces, id)
return(z)
}
tess2SP = function(x) {
stopifnot(is.tess(x))
y = tiles(x)
nam = names(y)
z = list()
for(i in seq(y))
z[[i]] = owin2Polygons(y[[i]], nam[i])
return(SpatialPolygons(z))
}
owin2SP = function(x) {
stopifnot(is.owin(x))
y = owin2Polygons(x)
z = SpatialPolygons(list(y))
return(z)
}First, let’s plot our original raster. The clickpoly function from
spatstat allows you to click and define polygon vertices directly on
your R graphics device. We are going to create a building with four
vertices, and then use our helper functions to convert those vertices
into a spatial polygon. Then we rasterize that polygon onto the same
grid as our original data with the rasterize function, which also sets
the height of our skyscraper. Finally, using merge, we will combine
them into a single dataset and reduce the
resolution.
plot(phillyraster)
# ployval = clickpoly(nv = 4, add = TRUE)
# saveRDS(ployval,"polygon.Rds")
download.file("https://dl.dropboxusercontent.com/s/76u8ih1njjpxsr0/polygon.Rds", destfile = "polygon.Rds")
## Alternate Link
#download.file("https://www.tylermw.com/data/polygon.Rds", destfile = "polygon.Rds")
# Always check serialized objects for dangerous code before running them:
ployval = readRDS("polygon.Rds")
str(ployval)## List of 5
## $ type : chr "polygonal"
## $ xrange: num [1:2] 2686745 2686941
## $ yrange: num [1:2] 234829 235042
## $ bdry :List of 1
## ..$ :List of 2
## .. ..$ x: num [1:4] 2686941 2686832 2686745 2686848
## .. ..$ y: num [1:4] 234955 235042 234911 234829
## $ units :List of 3
## ..$ singular : chr "unit"
## ..$ plural : chr "units"
## ..$ multiplier: num 1
## ..- attr(*, "class")= chr "unitname"
## - attr(*, "class")= chr "owin"
test2 = rasterize(owin2SP(ployval), phillyraster, 400)
philly_with_building = merge(test2, phillyraster)
philly_with_building_mat = raster_to_matrix(philly_with_building)
philly_with_building_mat_small = reduce_matrix_size(philly_with_building_mat,0.5)Now we’ll run the same analysis on the new lidar + building elevation data. Here’s where the benefit of using a programming language instead of a GUI to do GIS work becomes clear: to rerun this analysis with new data, we just change the variable names but run the same code.
getSunlightTimes(as.Date("2019-06-21"), lat = 39.9526, lon = -75.1652,tz = "EST")## date lat lon solarNoon nadir
## 1 2019-06-21 39.9526 -75.1652 2019-06-21 12:03:39 2019-06-21 00:03:39
## sunrise sunset sunriseEnd
## 1 2019-06-21 04:33:22 2019-06-21 19:33:57 2019-06-21 04:36:38
## sunsetStart dawn dusk
## 1 2019-06-21 19:30:40 2019-06-21 04:00:31 2019-06-21 20:06:47
## nauticalDawn nauticalDusk nightEnd
## 1 2019-06-21 03:18:49 2019-06-21 20:48:29 2019-06-21 02:30:09
## night goldenHourEnd goldenHour
## 1 2019-06-21 21:37:09 2019-06-21 05:14:06 2019-06-21 18:53:13
#Start and hour after sunrise and end an hour before sunset
philly_time_start = ymd_hms("2019-06-21 05:30:00", tz = "EST")
philly_time_end= ymd_hms("2019-06-21 18:30:00", tz = "EST")
# temptime = philly_time_start
# philly_new_shadows = list()
# sunangles = list()
# counter = 1
#
# while(temptime < philly_time_end) {
# sunangles[[counter]] = suncalc::getSunlightPosition(date = temptime, lat = 39.9526, lon = -75.1652)[4:5]*180/pi
# print(temptime)
# philly_new_shadows[[counter]] = ray_shade(philly_with_building_mat_small,
# sunangle = sunangles[[counter]]$azimuth+180,
# sunaltitude = sunangles[[counter]]$altitude,
# lambert = FALSE, zscale = 2,
# multicore = TRUE)
# temptime = temptime + duration("3600s")
# counter = counter + 1
# }
#Downloading the pre-computed data to save time
download.file("https://dl.dropboxusercontent.com/s/rsfyplworlbi4x7/philly_new_shadows.Rds",
destfile = "philly_new_shadows.Rds")
## Alternate Link
# download.file("https://www.tylermw.com/data/philly_new_shadows.Rds",
# destfile = "philly_new_shadows.Rds")
philly_new_shadows = readRDS("philly_new_shadows.Rds")
str(philly_new_shadows)## List of 13
## $ : num [1:1320, 1:1320] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
## $ : num [1:1320, 1:1320] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
## $ : num [1:1320, 1:1320] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ : num [1:1320, 1:1320] 1 1 1 1 1 1 1 1 1 1 ...
new_shadow_coverage = Reduce(`+`, philly_new_shadows)/length(philly_new_shadows)
new_shadow_coverage %>%
reshape2::melt(varnames = c("x","y"), value.name = "light") %>%
ggplot() +
geom_raster(aes(x = x,y = y,fill = light)) +
scale_fill_viridis_c("Daily Light\nCoverage", limits = c(0,1)) +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
coord_fixed()
Let’s plot the decrease in total daily sunlight intensity when we add the skyscraper. To do that, we’ll take the difference between the current shadows and the shadows calculated with the building present. Any negative values here represent an increase in sunlight, and only occur on top of the skyscraper, so we’ll zero those out. We’ll then subtract the resulting difference from 1 to show the relative total daily reduction in sunlight at each point.
shadow_difference = (shadow_coverage - new_shadow_coverage)
shadow_difference[shadow_difference < 0] = 0
(1 - shadow_difference) %>%
reshape2::melt(varnames = c("x","y"), value.name = "light") %>%
ggplot() +
geom_raster(aes(x = x,y = y,fill = light)) +
scale_fill_viridis_c("Relative Daily\nLight Coverage", limits = c(0,1)) +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
coord_fixed() +
ggtitle("Total Daily Sunlight Reduction from \nHypothetical West Philly Skyscraper")
This visualization is okay, but without building outlines there’s a lot of context missing. Luckily, we can use a quick trick to extract building outlines from the elevation data–calculate the normal vector (direction) at each point on the surface, and then only plot those points that are tilted extremely horizontally. If you don’t know what’s going on here, it’s okay–this is just a fun little hack.
# Calculate normal vectors and extract z-component (we also need to reorient
# the matrix--this should be fixed in future versions of rayshader)
abs(calculate_normal(building_mat_small,zscale = 4)$z) %>%
t() %>%
rayshader:::fliplr() ->
building_outlines
plot_map(building_outlines)
# Filter only those points that are facing at least 45 degrees outwards (the z-component
# of the unit-length normal vector is at most 1/2)
building_outlines %>%
reshape2::melt(varnames = c("x","y"), value.name = "normal") %>%
dplyr::filter(normal < 0.5) ->
existing_buildings_df
(1 - shadow_difference) %>%
reshape2::melt(varnames = c("x","y"), value.name = "light") %>%
ggplot() +
geom_raster(aes(x = x,y = y,fill = light)) +
geom_tile(data = existing_buildings_df,aes(x = x,y = y),fill = "black") +
scale_fill_viridis_c("Relative Daily\nLight Coverage", limits = c(0,1)) +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
coord_fixed() +
ggtitle("Total Daily Sunlight Reduction from \nHypothetical West Philly Skyscraper")
We can see this skyscraper doesn’t have much of an effect east of the train tracks due to the stadium and number of buildings already there, but it significantly reduces the amount of sunlight Penn Park recieves during the afternoon. I think we’ll nix this skyscraper.
For fun, we can also plot a 3D visualization of a time slice from the actual shadows from that day:
philly_with_building_mat_small %>%
sphere_shade(colorintensity = 5, sunangle = 93) %>%
add_shadow(philly_new_shadows[[3]],0.3) %>%
add_water(detect_water(philly_with_building_mat_small,zscale = 4,
max_height= 10, cutoff = 0.93, min_area = 1000)) %>%
plot_3d(philly_with_building_mat_small,zscale = 2,
zoom = 0.5,phi = 30,fov = 70,
windowsize = 1000, background = "#d9e7ff", shadowcolor = "#313947")
render_snapshot(title_text = "2019-06-21 8:30 AM, Philadelphia, Hypothetical Skyscraper",
title_bar_color = "white",
title_bar_alpha = 0.7,
vignette = 0.2)
rgl::rgl.close()I’m not going to run the following, but the code below will produce a 3D rotating animation as the shadows move through the day.
#download.file("https://dl.dropboxusercontent.com/s/gqnfs71t0u3pmu7/sunangles.Rds",
# destfile = "sunangles.Rds")
## Alternate Link
# download.file("https://www.tylermw.com/data/sunangles.Rds",
# destfile = "sunangles.Rds")
#
# sunangles = readRDS("sunangles.Rds")
# str(sunangles)
#
# transition_frames = round(seq(1,360,length.out = 14))[-14]
# transition_frames
# counter = 1
# for(i in 1:360) {
# if(i %in% transition_frames) {
# rgl::rgl.clear()
# philly_with_building_mat_small %>%
# sphere_shade(colorintensity = 5, sunangle = sunangles[[counter]]$azimuth+180) %>%
# add_shadow(philly_new_shadows[[counter]],0.3) %>%
# add_water(detect_water(philly_with_building_mat_small,zscale = 4,
# max_height= 10, cutoff = .93, min_area = 1000)) %>%
# plot_3d(philly_with_building_mat_small,zscale = 2,
# zoom = 0.5,phi = 30,fov = 70,
# windowsize = 1000, background = "#d9e7ff", shadowcolor = "#313947")
# counter = counter + 1
# }
# render_camera(theta = i)
# render_snapshot(glue::glue("rotating{i}"))
# }
#
# av::av_encode_video(glue::glue("rotating{1:60}.png"), output = "westphilly_movie.mp4",
# framerate = 30)For the sake of speed, this analysis was performed with hourly slices,
but you could improve the analysis by stepping down to a much smaller
timescale (e.g. temptime = temptime + duration("180s") for three
minute intervals). You could also remove the reduce_matrix_size() step
for increased resolution.
That’s it!
And that’s it for this overview of rayshader! I hope this was enough to get your started on your beautiful journey through the land of cartography and 3D mapping. If you have any questions, ask me in real life or on Twitter, and be sure to post any visualizations you make with the #MusaMasterclass, #rayshader, and #rstats hashtags. I’m always impressed and enjoy seeing what people manage to create!
Good luck, and happy shading!