30 Day Map Challenge, Day 30: Home (Pennsylvania)

pennsylvania_rayshader_example.R

library(rayshader)
library(raster)
library(stringr)
library(sp) 
library(maptools)
library(dplyr)
library(rgdal)
library(sf)
readShapePoly("eastern-states")

#Read in files in local directory without .zip in the name (I unzipped them all already)
list.files() %>%
  (function(x) x[str_detect(x, "zip", negate = TRUE)]) ->
pa_files

#Load all SRTM tiles into R
pa_dems = list()
for(i in 1:length(pa_files)) {
  pa_dems[[i]] = raster::raster(pa_files[[i]])
}

#Merge all SRTM tiles into one raster
pa_all = do.call(raster::merge, pa_dems)

# Download the shapefile of Pennsylvania to crop the SRTM data
# us = shapefile("USA_adm1.shp")
us = getData("GADM", country="USA", level=1)
nestates = "Pennsylvania"
ne = us[match(toupper(nestates), toupper(us$NAME_1)),]

# Crop and mask the raster with the vector shapefile
elevation.sub = raster::crop(pa_all, extent(ne))
elevation.sub.masked = mask(elevation.sub, ne)

# Convert the raster to a matrix and reduce the resolution by a factor of 5
elevation.sub.matrix = raster_to_matrix(elevation.sub.masked)
small_pa = reduce_matrix_size(elevation.sub.matrix, 0.2)

# Now load water shapefile for PA from National Hydrology Database

fgdb = "NHD_H_Pennsylvania_State_GDB/NHD_H_Pennsylvania_State_GDB.gdb"

# List all feature classes in a file geodatabase
fc_list = ogrListLayers(fgdb)
print(fc_list)

# Choose NHDArea for polygons representing water and crop the raster to PA region, using the 
# `elevation.sub.matrix` raster
fc <- readOGR(dsn=fgdb,layer="NHDArea")
water_cropped = crop(fc, elevation.sub.masked)

# Filter out all the small features (I just want the major rivers)
water_cropped_filtered = water_cropped[water_cropped$AreaSqKm > 0.01,]

# Convert spatial dataframe to sf object and rasterize with the fasterize package
water_sf = st_as_sf(water_cropped_filtered)
water_map = fasterize(water_sf, elevation.sub.masked, fun = "any")

# Convert to a matrix and set all NA's to 0 and anything else to 1 for `detect_water()`
water_map_mat = raster_to_matrix(water_map)
water_map_mat[!is.na(water_map_mat)] = 1
water_map_mat[is.na(water_map_mat)] = 0

# Reduce by a factor of 5 to match elevation data
water_map_mat_small = reduce_matrix_size(water_map_mat,0.2)

# Plot map
small_pa %>%
  sphere_shade() %>%
  add_water(rayshader:::fliplr(water_map_mat_small),"blue") %>%
  plot_map()

# Plot map in 3D
small_pa %>%
  sphere_shade(zscale=10) %>%
  add_water(rayshader:::fliplr(water_map_mat_small), "#4287f5") %>%
  add_shadow(ray_shade(small_pa, zscale=10, multicore = TRUE), 0.3) %>%
  add_shadow(ambient_shade(small_pa, zscale = 10, multicore = TRUE), 0.2) %>%
  plot_3d(small_pa, zscale=10 , windowsize = 1000, 
          background = "#d9f6ff", shadowcolor = "#3a5a63", shadowdepth = -150)
          
# Add text label
render_label(small_pa,"Horsham",x=3900,y=1550,z=5000,clear_previous = TRUE, linewidth = 5, textsize = 2,
             textcolor="black", linecolor="black",zscale=10)
       
# Rotate and create animation

for(theta in seq(1, 360, by = 1)) {
  render_camera(zoom = 0.9, fov = 0, phi = 45, theta = 45 - theta)
  render_snapshot(title_text = "My Hometown, Pennsylvania", 
                  title_bar_color = "black", title_color = "white",
                  title_size = 40, title_bar_alpha = 0.7, filename = glue::glue("horsham{theta}"))
}