# https://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html
# the github version of ggmap, which recently pulled in a small fix I had for a bug 
# devtools::install_github("dkahle/ggmap")
 
library(ggmap)
library(maps)
library(mapdata)
library(gganimate)
load("/cloud/project/Data/CalFire.RData")
# glimpse(CalFire.df)
kable(head(CalFire.df))
OBJECTID Damage Street.Number Street.Name Street.Type Street.Suffix City State Zip.Code CAL.FIRE.Unit County Community Battalion Incident.Name Incident.Number Incident.Start.Date Hazard.Type Where.did.fire.start. What.started.fire. Structure.Defense.Actions.Taken Structure.Type Structure.Category Units.in.Structure..if.multi.unit. of.Damaged.Outbuildings…120.SQFT of.Non.Damaged.Outbuildings…120.SQFT Roof.Construction Eaves Vent.Screen Exterior.Siding Window.Pane Deck.Porch.On.Grade Deck.Porch.Elevated Patio.Cover.Carport.Attached.to.Structure Fence.Attached.to.Structure Distance…Propane.Tank.to.Structure Distance…Residence.to.Utility.Misc.Structure..gt..120.SQFT Fire.Name..Secondary. APN..parcel. Assessed.Improved.Value..parcel. Year.Built..parcel. Site.Address..parcel. GLOBALID Latitude Longitude
1 No Damage 8376 Quail Canyon Road Winters CA NA LNU Solano 8 Quail CALNU 008419 2020-06-06 0:00:00 Fire Single Family Residence Multi Story Single Residence 1 NA NA Asphalt Unenclosed Mesh Screen <= 1/8” Wood Single Pane Wood Wood No Patio Cover/Carport No Fence Quail 0101090290 510000 1997 8376 QUAIL CANYON RD VACAVILLE CA 95688 {193828F2-B0EB-4404-B80F-FD2B56549C78} 38 -122
2 Affected (1-9%) 8402 Quail Canyon Road Winters CA NA LNU Solano Quail CALNU 008419 2020-06-06 0:00:00 Fire Deck on Grade Unknown Hand Crew Fuel Break Single Family Residence Single Story Single Residence NA NA NA Asphalt Unenclosed Mesh Screen <= 1/8” Wood Multi Pane Masonry/Concrete No Deck/Porch No Patio Cover/Carport Combustible N/A Quail 0101090270 573052 1980 8402 QUAIL CANYON RD VACAVILLE CA 95688 {15E713E9-D33F-4013-A405-1787444FACBE} 38 -122
3 No Damage 8430 Quail Canyon Road Winters CA NA LNU Solano Quail CALNU 008419 2020-06-06 0:00:00 Fire Single Family Residence Single Story Single Residence NA NA NA Asphalt Enclosed Mesh Screen > 1/8” Wood Single Pane No Deck/Porch No Deck/Porch No Patio Cover/Carport No Fence Quail 0101090310 350151 2004 8430 QUAIL CANYON RD VACAVILLE CA 95688 {A9DC6A21-2CC4-493D-84DC-A6F7553A54FA} 38 -122
4 No Damage 3838 Putah Creek Road Winters CA NA LNU Solano Quail CALNU 008419 2020-06-06 0:00:00 Fire Single Family Residence Single Story Single Residence NA NA NA Asphalt Unenclosed Mesh Screen > 1/8” Wood Single Pane No Deck/Porch No Deck/Porch Combustible No Fence Quail 0103010240 134880 1981 3838 PUTAH CREEK RD WINTERS CA 95694 {41C17AAC-D4CE-423B-B3AA-EA37F0F8B406} 38 -122
5 No Damage 3830 Putah Creek Road Winters CA NA LNU Solano Quail CALNU 008419 2020-06-06 0:00:00 Fire Single Family Residence Single Story Single Residence NA NA NA Tile Enclosed Mesh Screen > 1/8” Wood Multi Pane Wood Wood Combustible No Fence Quail 0103010220 346648 1980 3830 PUTAH CREEK RD WINTERS CA 95694 {307B34F1-5447-4340-B9DD-D47CC03EED50} 38 -122
6 No Damage 3834 Putah Creek Road Winters CA NA LNU Solano Quail CALNU 008419 2020-06-06 0:00:00 Fire Single Family Residence Single Story Single Residence NA NA NA Asphalt Enclosed Mesh Screen > 1/8” Wood Single Pane Wood Wood Combustible No Fence Quail 0103010230 159577 1978 3834 PUTAH CREEK RD WINTERS CA 95694 {4BB35E2D-2A20-40F2-9A08-8217AA202887} 38 -122

This is an analysis of the CalFire data set listing “residential” fire incidents from 2013 to 2020.

Data Preparation

In the original data set, break up the Time field into its components: year, day, etc. Also, for convenience in working with the data, let’s restrict the columns and eliminate those that don’t seem useful for analytics.

CalFire.df <- CalFire.df %>% 
#  filter(Damage != "No Damage") %>% # keep only with damage
#  filter(Damage == "Destroyed (>50%)") %>% # only Destroyed
  rename(When = Incident.Start.Date) %>% # shorter name
  select(c(1:4, 7:11, 15:17, 37, 43:44)) %>% # fewer columns
  mutate(WDay = wday(When, label=TRUE, abbr=TRUE), Year = year(When), Week = week(When), Month = month(When, label=TRUE), Day = day(When))

kable(head(CalFire.df))
OBJECTID Damage Street.Number Street.Name City State Zip.Code CAL.FIRE.Unit County Incident.Number When Hazard.Type Fire.Name..Secondary. Latitude Longitude WDay Year Week Month Day
1 No Damage 8376 Quail Canyon Winters CA NA LNU Solano CALNU 008419 2020-06-06 0:00:00 Fire Quail 38 -122 Sat 2020 23 Jun 6
2 Affected (1-9%) 8402 Quail Canyon Winters CA NA LNU Solano CALNU 008419 2020-06-06 0:00:00 Fire Quail 38 -122 Sat 2020 23 Jun 6
3 No Damage 8430 Quail Canyon Winters CA NA LNU Solano CALNU 008419 2020-06-06 0:00:00 Fire Quail 38 -122 Sat 2020 23 Jun 6
4 No Damage 3838 Putah Creek Winters CA NA LNU Solano CALNU 008419 2020-06-06 0:00:00 Fire Quail 38 -122 Sat 2020 23 Jun 6
5 No Damage 3830 Putah Creek Winters CA NA LNU Solano CALNU 008419 2020-06-06 0:00:00 Fire Quail 38 -122 Sat 2020 23 Jun 6
6 No Damage 3834 Putah Creek Winters CA NA LNU Solano CALNU 008419 2020-06-06 0:00:00 Fire Quail 38 -122 Sat 2020 23 Jun 6

Create summaries

Each fire incident is listed separately as a row in this data set. It is useful to group and count the number of incidents by County and by week/month/year etc.

# convert county names to standard lower case to match names from map data below; also complete the data to include regions with 0 incidents (so that they show up in map plots)

county.Monthly <- rename(count(CalFire.df, Year, Month, County), Freq = n) %>% ungroup() %>% complete(Year, Month, County, fill = list(Freq = 0)) %>% mutate(County = apply(as.data.frame(County),1,tolower)) 

county.Yearly <- county.Monthly %>% group_by(County, Year) %>% summarise(Freq = sum(Freq))

# county.Yearly <- rename(count(CalFire.df, Year,County), Freq = n) %>% ungroup() %>% complete(Year,County, fill = list(Freq = 0)) %>% mutate(County = apply(as.data.frame(County),1,tolower)) 

county.Weekly <- rename(count(CalFire.df, Year, Week,County), Freq = n) %>% ungroup() %>% complete(Year, Week, County, fill = list(Freq = 0)) %>% mutate(County = apply(as.data.frame(County),1,tolower))

Build Map Data

Now we’ll get the California county map data (i.e., lat-long locations).

# map("county", regions="California") # draws it (not data)

# cm <- get_map(location= "California", source="stamen", maptype="toner", crop=FALSE) 

county.map <- map_data("county", interior=TRUE, regions="california") %>% 
  rename(County=subregion) # %>% 
#  filter(lat > 39) # when focusing on Northern Counties

and then extend each county with the corresponding Counts of the number of fire incidents.

county.Yearly <- inner_join(county.map, county.Yearly, by = "County") # [, c("County", "Freq")]

county.Monthly <- inner_join(county.map, county.Monthly, by = "County") # [, c("County", "Freq")]

county.Total <- county.Yearly %>% group_by(County, order, group, long, lat) %>% summarise(Total=sum(Freq))

# county.map <- merge(county.map, county.count[, c("County", "Freq")], by="County")

# county.map$Freq <- county.count$Freq[match(county.map$subregion, county.count$County)]

Visuals

I want to make a geo-map of the number of fire incidents by county, using a lat-long map data (from map_data) to create the boundaries of each county (using ggplot and geom_polygon) and color to denote the metric.

There are a few ways to generalize this plot: 1) the time layer e.g., monthly or yearly totals (i.e., which of the dataframes defined above will be used), 2) which counties are shown, from all to those just north of a given latitude (e.g., 39), and 3) the metric itself, e.g., all fire incidents, or just those where damage was tagged in a certain way (another filter on the data frame).

I would like to write up a function that produces the plot with the above 3 generalizations in mind. But, the first one is not easy to pull off because the Yearly plot needs a facet_wrap while a Monthly needs a facet_grid addition. The third one is problematic because it requires a condition on Damage, but this column is not captured in the data frame used.

Yearly

# https://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html

# county.totals.plot(county.Monthly, 39, damage != "No Damage")
county.Yearly.plot <- function(latmin) {
  df <- county.Yearly %>% filter(lat > latmin)
  
  ggplot()  + geom_polygon(data = df, aes(x=long, y = lat, group=group), fill="white", color = "brown") + geom_polygon(data = df, aes(x=long, y = lat, group=group, fill=log(1+Freq)), color = "brown", size=0.5)  + scale_fill_gradient(low = "white", high = "blue")  + facet_wrap(~Year) + coord_fixed(1.3) + theme_classic() #  + theme(legend.position = "none")
}
county.Yearly.plot(20)

p1 <- ggplot()  + geom_polygon(data = county.Yearly, aes(x=long, y = lat, group=group), fill="white", color = "brown") + geom_polygon(data = county.Yearly, aes(x=long, y = lat, group=group, fill=log(1+Freq)), color = "brown", size=0.5)  + scale_fill_gradient(low = "white", high = "blue")  + facet_wrap(~Year) + coord_fixed(1.3) + theme_classic() #  + theme(legend.position = "none")

# + geom_polygon(data = county.map, aes(x=long, y = lat, group=group), color = "brown")

grid.arrange(p1,  ncol=1)

Animate by Year.

p1 <- ggplot()  + geom_polygon(data = county.Yearly, aes(x=long, y = lat, group=group), fill="white", color = "brown") + geom_polygon(data = county.Yearly, aes(x=long, y = lat, group=group, fill=log(1+Freq)), color = "brown", size=0.5)  + scale_fill_gradient(low = "white", high = "blue")  + coord_fixed(1.3) + theme_classic() #  + theme(legend.position = "none")

p1 + transition_time(Year)

Monthly

# county.totals.plot(county.Monthly, 39, damage != "No Damage")
county.Monthly.plot <- function(latmin) {
  df <- county.Monthly %>% filter(lat > latmin)
  
  ggplot()  + geom_polygon(data = df, aes(x=long, y = lat, group=group), fill="white", color = "brown") + geom_polygon(data = df, aes(x=long, y = lat, group=group, fill=log(1+Freq)), color = "brown", size=0.5)  + scale_fill_gradient(low = "white", high = "blue")  + facet_grid(vars(Year), vars(Month)) + coord_fixed(1.3) + theme_classic() #  + theme(legend.position = "none")
}
county.Monthly.plot(39)

# https://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html

p1 <- ggplot()  + geom_polygon(data = county.Monthly, aes(x=long, y = lat, group=group), fill=NA, color = "brown") + geom_polygon(data = county.Monthly, aes(x=long, y = lat, group=group, fill=log(1+Freq)), color = "brown", size=0.5)  + scale_fill_gradient(low = "white", high = "blue")  + facet_grid(vars(Year), vars(Month)) + coord_fixed(1.3) + theme_classic() #  + theme(legend.position = "none")

# + geom_polygon(data = county.map, aes(x=long, y = lat, group=group), color = "brown")

grid.arrange(p1,  ncol=1)