Lesson 6. Calculate and plot difference normalized burn ratio (dNBR) from Landsat remote sensing data in R

Learning Objectives

After completing this tutorial, you will be able to:

  • Calculate dNBR in R
  • Describe how the dNBR index is used to quantify fire severity.

What you need

You will need a computer with internet access to complete this lesson and the data for week 6 of the course.

Download Week 6/7 Data (~500 MB)

As mentioned in the previous lesson, we can use NBR to map the extent and severity of a fire. Let’s explore creating NBR using Landsat data.

Calculate dNBR using Landsat data

First, let’s setup our spatial packages.

# load spatial packages
library(raster)
library(rgdal)
library(rgeos)
library(RColorBrewer)
# turn off factors
options(stringsAsFactors = F)

Next, we open up our landsat data and create a spatial raster stack.

# create stack
all_landsat_bands_pre <- list.files("data/week_06/Landsat/LC80340322016189-SC20170128091153/crop",
           pattern=glob2rx("*band*.tif$"),
           full.names = T) # use the dollar sign at the end to get all files that END WITH
all_landsat_bands_pre
## [1] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band1_crop.tif"
## [2] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band2_crop.tif"
## [3] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band3_crop.tif"
## [4] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band4_crop.tif"
## [5] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band5_crop.tif"
## [6] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band6_crop.tif"
## [7] "data/week_06/Landsat/LC80340322016189-SC20170128091153/crop/LC80340322016189LGN00_sr_band7_crop.tif"

# stack the data
landsat_stack_pre <- stack(all_landsat_bands_pre)

Next we calculate dNBR using the following steps:

  1. Open up pre-fire data and calculate NBR
  2. Open up the post-fire data and calculate NBR
  3. Calculate dNBR (difference NBR) by subtracting post-fire NBR from pre-fire NBR.
  4. Classify the dNBR raster using the classification table below.

. Note the code to do this is hidden. You will need to figure out what bands are required to calculate NBR using Landsat.

landsat derived NDVI plot

You can export the NBR raster if you want using writeRaster().

writeRaster(x = landsat_nbr_pre,
              filename="data/week_06/outputs/landsat_nbr",
              format = "GTiff", # save as a tif
              datatype='INT2S', # save as a INTEGER rather than a float
              overwrite = T)

Next, we can open the post-fire landsat data to calculate post-fire NBR.

all_landsat_bands_post <- list.files("data/week_06/Landsat/LC80340322016205-SC20170127160728/crop",
           pattern=glob2rx("*band*.tif$"),
           full.names = T) # use the dollar sign at the end to get all files that END WITH
all_landsat_bands_post
## [1] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band1_crop.tif"
## [2] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band2_crop.tif"
## [3] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band3_crop.tif"
## [4] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band4_crop.tif"
## [5] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band5_crop.tif"
## [6] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band6_crop.tif"
## [7] "data/week_06/Landsat/LC80340322016205-SC20170127160728/crop/LC80340322016205LGN00_sr_band7_crop.tif"

# stack the data
landsat_stack_post <- stack(all_landsat_bands_post)

Then we calculate NBR on the post data - note the code here is purposefully hidden. You need to figure out what bands to use to perform the math!

landsat derived NBR post fire

Finally, calculate the DIFFERENCE between the pre and post NBR!!

# calculate difference
landsat_nbr_diff <- landsat_nbr_pre - landsat_nbr_post
plot(landsat_nbr_diff,
     main = "Difference NBR map \n Pre minus post Cold Springs fire",
     axes=F, box=F)

Difference NBR map

When you have calculated dNBR or the difference in NBR pre minus post fire, classify the output raster using the classify() function and the classes below.

SEVERITY LEVEL dNBR RANGE
Enhanced Regrowth > -.1
Unburned -.1 to +.1
Low Severity +.1 to +.27
Moderate Severity +.27 to +.66
High Severity > +1.3

NOTE: your min an max values for NBR may be slightly different from the table shown above! If you have a smaller min value (< -700) then adjust your first class to that smallest number. If you have a largest max value (>1300) then adjust your last class to that largest value in your data.

Alternatively, you can use the Inf to specify the smallest -Inf and largest Inf values.

Your classified map should look something like:

classified NBR output

Compare to fire boundary

As an example to see how our fire boundary relates to the boundary that we’ve identified using MODIS data, we can create a map with both layers. I’m using the shapefile in the folder:

data/week_06/vector_layers/fire-boundary-geomac/co_cold_springs_20160711_2200_dd83.shp

Add fire boundary to map.

classified NBR output

Make it look a bit nicer using a colorbrewer palette. I used the RdYlGn palette:

brewer.pal(5, 'RdYlGn')

classified NBR output

Note that you will have to figure out what date these data are for! I purposefully didn’t include it in the title of this map.

barplot(nbr_classified,
        main = "Distribution of Classified NBR Values",
        col=the_colors)

plot barplot of fire severity values

Add labels to your barplot!

barplot(nbr_classified,
        main = "Distribution of Classified NBR Values",
        col=the_colors,
        names.arg = c("Enhanced \nRegrowth", "Unburned", "Low \n Severity", "Moderate \n Severity", "High \nSeverity"))

plot barplot of fire severity values with labels

Optional challenge - NBR using MODIS

The table below shows the band ranges for the MODIS sensor. We know that the NBR index will work with any multispectral sensor with a NIR band between 760 - 900 nm and a SWIR band between 2080 - 2350 nm. What bands should we use to calculate NBR using MODIS?

BandWavelength range (nm)Spatial Resolution (m)Spectral Width (nm)
Band 1 - red620 - 6702502.0
Band 2 - near infrared841 - 8762506.0
Band 3 - blue/green459 - 4795006.0
Band 4 - green545 - 5655003.0
Band 5 - near infrared1230 – 12505008.0
Band 6 - mid-infrared1628 – 165250018
Band 7 - mid-infrared2105 - 215550018

Leave a Comment