GEOG 208 Final: Woolsey Fire

This project was created in order to demonstrate the data visualization and interactive mapping capabilities of R. For this project I will be using the November 2018 Woolsey Fire as a case study to demonstrate how R can help in spatial analysis of fire data as well as in displaying that data to a broader audience. I chose this case study since it can benefit from spatial layering of different data formats and multiple levels of analysis can be done to it.

I will be using data aggregatd from multiple sources, all data sources referenced in this project are free to access for anyone interested in duplicating this project or applying its methods towards a similar project. A majority of the documentation will be written in the code chunks to allow for better step-by-step explanation.

This code chunk shows all the packages and libraries that I used for this project.

##################Packages and Libraries##################
install.packages("magick")
library(magick)
library(devtools)
install_github("mtennekes/tmaptools")
install_github("mtennekes/tmap")
#install.packages("tmap")
library(tmap)
library(sp)
library(raster)
library(RColorBrewer)
library(classInt)
install.packages('rasterVis')
library(rasterVis)
install.packages('tigris')
library(tigris)
install.packages('rgdal', type = "source", configure.args=c('--with-proj-include=/Library/Frameworks/PROJ.framework/Headers', '--with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib'))
library(rgdal)
#install.packages("tmaptools")
install.packages("sf")
install.packages("leaflet")
library(tmaptools)
library("sf")
library(leaflet)
library(tigris)
install.packages('devtools')
devtools::install_github('thomasp85/gganimate')
library(ggplot2)
library(gganimate)
library(ggplot2)
library(ggmap)
library(maps)
library(ggthemes)
install.packages("plotly")
library(plotly)

Interactive tmap() map: For this map, I used the tmap package to create an interactive map that layers the Digital Elevation Model, median home value, and fire incidence data on one interactive map. The layer toggle also allows the user to select which layers to display in order to customize the map to only show information that the end user wants to see.

##################Add Fire Data from NASA FIRMS##################
#Create download requests here:https://firms.modaps.eosdis.nasa.gov/download/

#import fire points data
dir("/Users/krgr.df/Desktop/GEOG208/Final Assignment/Fire Data/DL_FIRE_V1_33042(SHP8-24)")
Fire.Points <- readOGR("/Users/krgr.df/Desktop/GEOG208/Final Assignment/Fire Data/DL_FIRE_V1_33042(SHP8-24)", "fire_nrt_V1_33042")

##################Merge and Reclassify Raster##################
#raster data downloaded from from glovis: https://glovis.usgs.gov/app

#import and merge two raster files
DEM_4 <- raster("/Users/krgr.df/Desktop/GEOG208/Final Assignment/DEMs/n34_w120_3arc_v2.tif")
DEM_5 <- raster("/Users/krgr.df/Desktop/GEOG208/Final Assignment/DEMs/n34_w119_3arc_v2.tif")
mosaic(DEM_4, DEM_5, fun = mean, tolerance=0.05, filename="/Users/krgr.df/Desktop/GEOG208/Final Assignment/DEMs/Mosaic_test_1.tif", overwrite=TRUE) 
mosaic_test <- raster("/Users/krgr.df/Desktop/GEOG208/Final Assignment/DEMs/Mosaic_test_1.tif", overwrite=TRUE)

#reclassify raster to set values below 500 to NA, this removes them from display
DEM <- calc(mosaic_test, fun=function(x){ x[x < 500] <- NA; return(x)} )

##################Add Demographic Data##################
#pull ventura & LA county boundary map from library(tigris)
Ven_bound <- tracts(state = "CA", county = "Ventura", cb = TRUE)

LA_bound <- tracts(state = "CA", county = "Los Angeles", cb = TRUE)

#import csv file with median home value data (Ventura)
HW4 <- file.path ("/Users/krgr.df/Desktop/GEOG208/Assignment 4/ventura_median_val.csv")
median.val <- read.csv(HW4, stringsAsFactors = FALSE)
median.val$median.value <- as.numeric(median.val$median.value)
median.val$id2 <- as.character(median.val$id2)

#import csv file with median home value data (Los Angeles)
HW4.2 <- file.path ("/Users/krgr.df/Desktop/GEOG208/Assignment 4/LA_median_val.csv")
LA.median.val <- read.csv(HW4.2, stringsAsFactors = FALSE)
LA.median.val$median.val <- as.numeric(LA.median.val$median.val)
LA.median.val$id2 <- as.character(LA.median.val$id2)

#merge census tract boundaries with median home values (Ventura)
Ventura.Stats <- merge(Ven_bound, median.val, by.x = "AFFGEOID", by.y = "id")

#merge census tract boundaries with median home values (Los Angeles)
LA.Stats <- merge(LA_bound, LA.median.val, by.x = "AFFGEOID", by.y = "id")

##################Create Interactive tmap()##################
#set tmap to interactive viewing mode and set theme to 'cobalt'
tmap_mode("view")
tmap_style("cobalt")

#allows toggling between static aka 'plot' and interactive aka 'view' modes
ttm() 

#the following lines of code adds layers into tmap and edits the appearance of the map
tmap_map <- tm_shape(Ventura.Stats)+
  tm_fill(col = "median.value", palette = "PuBu", title = "Ventura County Median Home Value", 
          style = "fixed",
          breaks = c(0, 250000, 500000, 750000, 1000000, Inf),
          textNA = "No Data", 
          colorNA = "gray",
          alpha = .7) +
  tm_borders(col = "white") +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.text.color = "gray20",
            legend.bg.color = "white",
            legend.width = -1,
            legend.bg.alpha = 1)+
  tm_view(view.legend.position = c("right", "bottom")) +
  tm_shape(LA.Stats)+
  tm_fill(col = "median.val", palette = "YlGn", title = "LA County Median Home Value", 
          style = "fixed",
          breaks = c(0, 250000, 500000, 750000, 1000000, Inf),
          textNA = "No Data", 
          colorNA = "gray",
          alpha = .7) +
  tm_borders(col = "white") +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.text.color = "gray20",
            legend.bg.color = "white",
            legend.width = -1,
            legend.bg.alpha = 1) +
  tm_view(view.legend.position = c("left", "bottom")) +
  tm_shape(Fire.Points) +
  tm_dots("BRIGHT_TI4", 
          alpha = 1,
          border.col = NA,
          title = "Fire Intensity",
          popup.vars = "ACQ_DATE",
          palette = "YlOrRd") +
  tm_view(view.legend.position = c("right", "top")) +
  tm_shape(DEM, unit = "mi") +
  tm_raster("layer", 
            midpoint = NA, 
            title = "DEM",
            legend.show = FALSE) +
  tmap_options(max.raster = c(plot = 2883601, view = 2883601)) +
  tm_scale_bar(size = 1, 
               text.color = NA,
               color.dark = "black", color.light = "white", 
               lwd = 1, 
               position = c("bottom", "left"))

tmap_map

I created a histogram of fire frequency over time in order to pinpoint which days had the highest number of fire incidents

##################Plot Histogram Data##################

#import NASA fire points .csv data
fire_csv <- read.csv("/Users/krgr.df/Desktop/GEOG208/Final Assignment/Fire Data/NASA VIIRs CSV/fire_nrt_V1_32986.csv")

#change acquired date column from numeric to date
fire_csv$acq_date <- as.Date(fire_csv$acq_date)

#plot a histogram of fire frequency over time
hist(fire_csv$acq_date, "days", 
     freq = TRUE,
     breaks = 15,
     xlab = "Date Acquired",
     col = "firebrick3",
     border = "orange"
     )

I also created an animated map that displays the distribution of fire instances over time during the extent of the Woolsey Fire.

##################Animated Map##################
#created using gganimate package from: https://github.com/thomasp85/gganimate

#change format of date column from numeric to date
fire_csv$acq_date <- as.Date(fire_csv$acq_date)
str(fire_csv)

#create animated map using gganimate
animated_map <- ggplot(fire_csv, (aes(longitude, latitude, color = bright_ti4))) +
  borders("county", ".*ventura",colour = "gray85", fill = "gray80", ) +
  theme_map() +
  borders("county", ".*los angeles",colour = "gray85", fill = "gray80", ) +
  theme_map() +
  geom_point(alpha = .5) +
  scale_color_gradient(low="yellow", high="red", name = "Brightness Value") + 
  labs(title = 'Date: {frame_time}', size = 3) +
  transition_time(acq_date) +
  enter_fade() + 
  ease_aes("linear")

animated_map

ANIMATED MAP HERE (open link in new tab): https://photos.google.com/share/AF1QipNZN2aIUorKUYKMkM5VbLYoWQUNbYvYlYubbeNgeh5SEIdfSBAkNcCuuEAu__RiHg/photo/AF1QipM8bCC5syWDq5n3_9v2GmFVQO7jzpSIIFAcZCIM?key=NDFycU16eExNSl81dUpRLXAwVEtWZmtFMkhuX3V3

A heatmap is useful in this case because the name fits perfectly with my case study (get it, “heat” “map” ok im done) and because its useful for seeing which areas had the highest concentration of fire incidents.

##################Heatmap##################
#created by using techniques from:
#http://trucvietle.me/r/tutorial/2017/01/18/spatial-heat-map-plotting-using-r.html

#import fire data in .csv format
fire_csv <- read.csv("/Users/krgr.df/Desktop/GEOG208/Final Assignment/Fire Data/NASA VIIRs CSV/fire_nrt_V1_32986.csv")

#set map extent by taking mean lat long of fire points (genius method created by Truc Viet)
mean.longitude <- mean(fire_csv$longitude)
mean.latitude <- mean(fire_csv$latitude)

#get Google basemap
if(!requireNamespace("devtools")) install.packages("devtools")
devtools::install_github("dkahle/ggmap", ref = "tidyup")
register_google(key = "AIzaSyDpOleZ5sm6jctsh7fHZ8vRKamcCW33nNM")
affected_area <- get_map(location = c(mean.longitude, mean.latitude), zoom = 9, scale = 2)
affected_area <- ggmap(affected_area, extent="device", legend="none")

#create heatmap, in-depth documentation from website linked above
heatmap1 <- affected_area + 
            stat_density2d(data=fire_csv,
                 aes(x=longitude, y=latitude, fill=..level.., alpha=..level..), 
                 geom="polygon") +
            scale_fill_gradientn(colours=(brewer.pal(7, "YlOrRd"))) +
            guides(size=FALSE, alpha = FALSE)
heatmap1

#transform static ggplot() heatmap to interactive map using ggplotly()
heatmaply <- ggplotly(heatmap1)
heatmaply

This final map is a GIF showing satellite imagery of the area before and after the fire, the fire points from the NASA dataset are also added to highlight which areas had a higher fire incident frequency.

#created by using magick package see:
#https://cran.r-project.org/web/packages/magick/vignettes/intro.html
#satellite imagery (Sentinel2) downloaded from Earth Explorer:
#https://earthexplorer.usgs.gov

#created by cycling through screenshotted images from QGIS
before_fire <- image_scale(image_read("/Users/krgr.df/Desktop/Screen Shot 2018-12-04 at 2.20.06 AM.png"), "x450")
before_points <- image_scale(image_read("/Users/krgr.df/Desktop/Screen Shot 2018-12-04 at 2.07.35 AM.png"), "x450")
with_points <- image_scale(image_read("/Users/krgr.df/Desktop/Screen Shot 2018-12-04 at 2.07.52 AM.png"), "x450")
frames <- image_morph(c( before_fire, with_points, before_points, before_points, before_points, before_fire), frames = 25)
image_animate(frames)

ANIMATED MAP HERE (open link in new tab): https://photos.google.com/share/AF1QipNZN2aIUorKUYKMkM5VbLYoWQUNbYvYlYubbeNgeh5SEIdfSBAkNcCuuEAu__RiHg/photo/AF1QipNkvlR-hjBYCTTb4qZYSz-cO5hbNM0GNOsES2rP?key=NDFycU16eExNSl81dUpRLXAwVEtWZmtFMkhuX3V3

Next Steps: For the next steps I plan on implementing other factors into the first interactive map such as wind strength and direction (especially the Santa Ana winds), infrastructure types, and vegetation/land use data.