To calculate difference Normalized Burn Ratio (dNBR), you first need to calculate NBR for the pre and post fire data. This, of course, presumes that you have data before and after the area was burned from the same remote sensing sensor. Ideally, this data also does not have clouds covering the fire area.

Open up and stack the Landsat post-fire data.

# Import and stack post fire Landsat dataall_landsat_bands=glob("data/cold-springs-fire/landsat_collect/LC080340322016072301T1-SC20180214145802/crop/*band*.tif")all_landsat_bands.sort()landsat_post_fire_path="data/cold-springs-fire/outputs/landsat_post_fire.tif"es.stack_raster_tifs(all_landsat_bands,landsat_post_fire_path)# You are not cropping the data in this lesson so you can just use .read()withrio.open(landsat_post_fire_path)assrc:landsat_post_fire=src.read(masked=True)landsat_post_meta=src.profilelandsat_post_bounds=src.boundslandsat_extent=plotting_extent(src)# Open fire boundary layer and reproject it to match the Landsat datafire_boundary_path="data/cold-springs-fire/vector_layers/fire-boundary-geomac/co_cold_springs_20160711_2200_dd83.shp"fire_boundary=gpd.read_file(fire_boundary_path)# If the CRS' are not the same be sure to reprojectfire_bound_utmz13=fire_boundary.to_crs(landsat_post_meta['crs'])

Next, you can calculate NBR on the post fire data. Remember that NBR uses different bands than NDVI but the calculation formula is the same. For landsat 8 data you will be using bands 7 and 5. And remember because python starts counting at 0 (0-based indexing), that will be bands 6 and 4 when you access them in your numpy array.

Below the es.normalized_diff() function is used to calculate NBR. You can also calculate NBR like you did NDVI with raster math :

Next, calculate NBR for the pre-fire data. Note that you will have to download the data that is being used below from Earth Explorer following the lessons in this tutorial series.

Also note that you will need to clip or crop the data so that you can subtract the post fire data from the pre fire data. The code to do this is hidden but you did this last week so you should know what to do!

True

Next you can:

calculate NBR on the pre-fire data and then

calculate difference dNBR by subtracting the post fire data FROM the pre (pre-post)

Finally you can classify the data. Remember that dNBR has a set of classification bins and associated categories that are commonly used. When you have calculated NBR - classify the output raster using the np.digitize() function. Use the dNBR 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

> +.66

NOTE: your min and max values for NBR may be slightly different from the table shown above. In the code example above, np.inf is used to indicate any values larger than .66.

The code to classify below is hidden! You learned how to classify rasters in week 2 of this course.

Finally you are ready to plot your data. Note that legends in python can be tricky so you are provided with two legend options below. In the first plot, you create a legend with individual boxes. Creating a legend this way forces you to use matplotlib patches. Effectively you are drawing each box to create a legend from colors you’ve used in your image.

a ListedColormap: this is a discrete colormap. Thus it will only contain n colors rather than a gradient of colors.

A list of titles for each category in your classification scheme: this is what will be “written” on your legend

You are ready to plot. Below you add a custom legend with unique boxes for each class. This is one way to create a legend. It requires a bit more knowledge of matplotlib under the hood to make this work!

Classified difference normalized burn ratio (dNBR) calculated for the Cold Springs fire images from Landsat, with legend created using earthpy.

Create a Colorbar Legend

Alternatively, you can create a discrete colorbar with labels. This method might be a bit less technical to follow. You can decide what type of legend you prefer for your homework.

# Grab raster unique values (classes)values=np.unique(dnbr_landsat_class).tolist()# Add another index value because for n categories you need n+1 values to create binsvalues=[0]+values# Make a color map of fixed colorsnbr_colors=["g","yellowgreen","peachpuff","coral","maroon"]nbr_cmap=ListedColormap(nbr_colors)# But the goal is the identify the MIDDLE point of each bin to create a centered tickbounds=[((a+b)/2)fora,binzip(values[:-1],values[1::1])]+[5.5]# Define normalizationnorm=colors.BoundaryNorm(bounds,nbr_cmap.N)

Then, plot the data but this time add a colorbar rather than a custom legend.

Calculate Total Area of Burned Area

Once you have classified your data, you can calculate the total burn area.

For this lesson, you are comparing the area calculated using modis that is “high severity” to that of Landsat.

You can calculate this using a loop, a function or manually - like so:

# To calculate area, multiply the number of pixels in each bin by image resolution# Result will be in total square meterslandsat_pixel=landsat_pre.res[0]*landsat_pre.res[0]burned_landsat=(dnbr_landsat_class[dnbr_landsat_class==5]).sizeburned_landsat=np.multiply(burned_landsat,landsat_pixel)print("Landsat Severe Burn Area:",burned_landsat)