For this project, I was interested in seeing how Rstudio’s functionalities could be used to classify and potentially predict wildfire burn severity. Fire severity is notably defined by USGS research scientist Jon Keeley “the loss or decomposition of organic matter aboveground and belowground” after a fire (Keeley 2009). So, using this definition, high-severity fire is fire that kills a high proportion of life all through the critical zone of a given ecosystem. Due to a confluence of factors surrounding anthropogenic climate change increasing the total amount of fire-weather days a year and a century of fire suppression increasing fuels density, high-severity fire has increased eightfold in California from 1985 to the present (Williams et al. 2018, Parks et al. 2020). Measuring fire severity is important in predicting post-fire vegetative succession, prioritizing areas for fuels treatments, and identifying areas that may be susceptible to post fire landslides, soil erosion, and flooding (Neary et al. 2011).

A common way of stratifying burn severity in remote sensing is through computation of the Normalized Burn Ratio (NBR). This ratio is calculated through taking advantage of the unique reflective properties of vegetation and recently burned areas – healthy vegetation has very high spectral reflectance in the near infrared and very low spectral reflectance in the shortwave infrared, while for burned areas the inverse is true. Because of this the formula (NIR – SWIR) / (NIR + SWIR) allows for clear delineation of vegetation loss post-fire – this formula is the NBR.

I decided to download Landsat 8 data of areas burnt by the Dixie Fire of 2021, a fire that burnt over 380,000 hectares in Plumas and Lassen National Forests. I downloaded images from July 11th and October 25th, before and after the Dixie fire, controlling for images that had less than 10% cloud cover. I first imported this imagery into QGIS and ran a quick NBR classification (something I already knew how to do in GIS software) in order to choose an area with interesting topography and a wide range of potential burn severity classes represented. The study area ended up being a 35,000 hectare parcel that included burnt areas around Bald Eagle Mountain in Plumas. Below you see the process of importing bands 5 and 7 of the images in question, the two Landsat 8 bands needed for Normalized Burn Ratio classification. Additionally, I set the boundary based on the study area I had chosen in QGIS using the “Clip Raster” function in QGIS.

boundary <- ext(634736, 656429, 4413482, 4430448)

Band5Pre <- rast("LC08_L2SP_044032_20210711_20210720_02_T1_SR_B5.TIF")
Band7Pre <- rast("LC08_L2SP_044032_20210711_20210720_02_T1_SR_B7.TIF")
Band5Post <- rast("LC08_L2SP_044032_20211015_20211025_02_T1_SR_B5.TIF")
Band7Post <- rast("LC08_L2SP_044032_20211015_20211025_02_T1_SR_B7.TIF")

Band5Pre <- crop(Band5Pre, boundary)
Band5Post <- crop(Band5Post, boundary)
Band7Pre <- crop(Band7Pre, boundary)
Band7Post <- crop(Band7Post, boundary)

plot(Band5Pre,
     main = "Near-Infrared Imagery of Study Area Pre-Fire")

plot(Band5Post,
     main = "Near-Infrared Imagery of Study Area Post-Fire")

plot(Band7Pre,
     main = "Shortwave Infrared Imagery of Study Area Pre-Fire")

plot(Band7Post,
     main = "Shortwave Infrared Imagery of Study Area Post-Fire")

Next up was the raster algebra potion of my project. In order to get a dNBR to classify fire severity, you must use the equation (NIR - SWIR) / (NIR + SWIR) for each image in question. You can then subtract the NBR after fire from the NBR before fire in order to derive a difference NBR raster. This resultant raster can then have its symbology altered so that burn severity classes are apparent. The NBR ranges you see in my code below come directly from the USGS’s NBR burn severity instructions (USGS). Additionally, some of these methods for classification were inspired by a walkthrough from the UN’s Office for Outer Space Affairs UN-SPIDER Knowledge Portal.

NBRpre <- (Band5Pre - Band7Pre) / (Band5Pre + Band7Pre)
NBRpost <- (Band5Post - Band7Post) / (Band5Post + Band7Post)
dNBR <- (NBRpre) - (NBRpost)

dNBR_scaled <- 1000*dNBR

class.matrix <- matrix(c(-251, -101, 2, -101, 99, 3, 99, 269, 4, 269, 439, 5, 439, 659, 6, 659, 1300, 7), ncol=3, byrow = TRUE) 

dNBR_reclass <- classify(dNBR_scaled, rcl = class.matrix, right=NA)

leg  <- c("Enhanced Regrowth, Low", "Unburned", "Low Severity", "Moderate-low Severity", "Moderate-high Severity", "High Severity") 

levels(dNBR_reclass) <- leg

my_col=c("darkolivegreen4","limegreen", "yellow2", "orange2", "red", "purple")

plot(dNBR_reclass, type="classes", plg=list(legend=c("Enhanced Regrowth, Low", "Unburned", "Low Severity", "Moderate-low Severity", "Moderate-high Severity", "High Severity") , x="bottomright"), main = "dNBR Stratified Burn Severity \nDixie Fire, Plumas National Forest", col = my_col)

hist(dNBR_scaled,
     main = "Histogram of dNBR Values, Scaled by 1000",
     xlab = "dNBR Value",
     col = "dark red")

I wanted to add a hillshade to my dNBR map in order to get a sense of the influence of topography on fire severity in my AOI. Below you see the import and treatment of DEM, the creation of a hillshade, and my final dNBR map of study area present with tmapper in interactive viewing mode. Zooming in on this classified image you can see that there are indeed a variety of severity classes represented, with high severity to the west of highway 70 and the highest severity appearing at the easternmost crests of a couple of medium-sized ridges.

DEM <- rast("AOIDEM.tif")
crs(DEM) <- "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"

AOIDEM <- crop(DEM, boundary)

slope <- terrain(AOIDEM, v="slope")
aspect <- terrain(AOIDEM, v="aspect")
hilsh <- Pshade <- shade(slope/180*pi, aspect/180*pi, angle=40, direction=330)

plot(AOIDEM)

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(hilsh) + tm_raster(palette = "-Greys",legend.show = F) +
tm_shape(dNBR_reclass) +
  tm_raster(style = "cat", labels = leg, palette = my_col, alpha = 0.65, legend.show = T) +
  tm_layout(title = "dNBR Stratified Burn Severity, Dixie Fire, Plumas National Forest")
## stars object downsampled to 1312 by 762 cells. See tm_shape manual (argument raster.downsample)

After creating this map, I was interested in seeing if I could find continuous rasters on fuel loads in the interest of creating a rudimentary predictive linear model of fire severity. I found three rasters through LandFire, presented by the Forest Service, on Canopy Base Height, Canopy Bulk Density, and Fuel Vegetation Height. The linear model seen here is more of a proof of concept with the data I wrangled than anything else, given how complex the issue of fire spread and fire severity is. The R-squared value for this glm is very low, around 0.02, however it is interesting to note that with each added raster the R-squared input does increase somewhat. To see the true statistical significance, sampling should be used to reduce the size of this data. I would love attempt this kind of sampling, as well as integrate data on surface fuel loads, if I continue analyzing this particular area.

CanopyBaseHeight <- rast("CBH.tif")
crs(CanopyBaseHeight) <- "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"

CanopyBulkDensity <- rast("CBD.tif")
crs(CanopyBulkDensity) <- "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"

FuelVegetationHeight <- rast("FVHEIGHT.tif")
crs(FuelVegetationHeight) <- "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"

CBH <- crop(CanopyBaseHeight, boundary)
CBD <- crop(CanopyBulkDensity, boundary)
FVH <- crop(FuelVegetationHeight, boundary)



lm1 <- lm(values(dNBR)[,1] ~ values(CBD)[,1] + values(CBH)[,1] + values(FVH)[,1])
summary(lm1)
## 
## Call:
## lm(formula = values(dNBR)[, 1] ~ values(CBD)[, 1] + values(CBH)[, 
##     1] + values(FVH)[, 1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33298 -0.13135 -0.02144  0.11400  0.53655 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.039e-01  1.002e-03  103.64   <2e-16 ***
## values(CBD)[, 1] -6.148e-04  7.738e-06  -79.45   <2e-16 ***
## values(CBH)[, 1] -5.612e-04  1.683e-05  -33.35   <2e-16 ***
## values(FVH)[, 1]  1.809e-04  2.101e-06   86.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1532 on 408491 degrees of freedom
## Multiple R-squared:  0.02184,    Adjusted R-squared:  0.02183 
## F-statistic:  3040 on 3 and 408491 DF,  p-value: < 2.2e-16

Thank you for your interest in my project!

References:

“Dixie Fire (CA) Information - InciWeb the Incident Information System.” Accessed December 10, 2021. https://inciweb.nwcg.gov/incident/7690/.

Keeley, Jon E., and Jon E. Keeley. “Fire Intensity, Fire Severity and Burn Severity: A Brief Review and Suggested Usage.” International Journal of Wildland Fire 18, no. 1 (February 17, 2009): 116–26. https://doi.org/10.1071/WF07049.

“LANDFIRE Program: Data Products - Fuels - Forest Canopy Bulk Density.” Accessed December 10, 2021. https://landfire.gov/cbd.php.

“Landsat Normalized Burn Ratio | U.S. Geological Survey.” Accessed December 13, 2021. https://www.usgs.gov/landsat-missions/landsat-normalized-burn-ratio.

Lydersen, Jamie M., Brandon M. Collins, Matthew L. Brooks, John R. Matchett, Kristen L. Shive, Nicholas A. Povak, Van R. Kane, and Douglas F. Smith. “Evidence of Fuels Management and Fire Weather Influencing Fire Severity in an Extreme Fire Event.” Ecological Applications 27, no. 7 (2017): 2013–30.

Meneses, Bruno M. “Vegetation Recovery Patterns in Burned Areas Assessed with Landsat 8 OLI Imagery and Environmental Biophysical Data.” Fire 4, no. 4 (December 2021): 76. https://doi.org/10.3390/fire4040076.

Neary, Daniel G, Karen A Koestner, and Anne Youberg. “Hydrologic Impacts of High Severity Wildfire: Learning from the Past and Preparing for the Future,” n.d., 8.

Odion, Dennis C., and Chad T. Hanson. “Fire Severity in Conifer Forests of the Sierra Nevada, California.” Ecosystems 9, no. 7 (November 1, 2006): 1177–89. https://doi.org/10.1007/s10021-003-0134-z.

Parks, S. A., and J. T. Abatzoglou. “Warmer and Drier Fire Seasons Contribute to Increases in Area Burned at High Severity in Western US Forests From 1985 to 2017.” Geophysical Research Letters 47, no. 22 (2020): e2020GL089858. https://doi.org/10.1029/2020GL089858.

“Step by Step: Burn Severity with R-Studio | UN-SPIDER Knowledge Portal.” Accessed December 14, 2021. https://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/Step-By-Step/RStudio.

Tran, Bang Nguyen, Mihai A. Tanase, Lauren T. Bennett, and Cristina Aponte. “Evaluation of Spectral Indices for Assessing Fire Severity in Australian Temperate Forests.” Remote Sensing 10, no. 11 (November 2018): 1680. https://doi.org/10.3390/rs10111680.

Williams, A. Park, John T. Abatzoglou, Alexander Gershunov, Janin Guzman-Morales, Daniel A. Bishop, Jennifer K. Balch, and Dennis P. Lettenmaier. “Observed Impacts of Anthropogenic Climate Change on Wildfire in California.” Earth’s Future 7, no. 8 (August 1, 2019): 892–910. https://doi.org/10.1029/2019EF001210.