In August 2021 Greece encountered severe wildfires, some of them resulting to an unbelievably enormous damaged area. I was watching the news and I was feeling so sad about the wildlife and the forests and the villages and the people’s properties all burnt down to ashes.
I was not only watching the news. I was daily watching the satellites and I was stuck in front of my PC waiting for the most updated release to come, to see the impact of the wildfires through multi-spectral imagery. My most frequent process was to download imagery, captured before a fire event and immediately after it, for a particular area, and then process it in ArcGIS Pro to monitor the damage and to calculate the burnt area.
In this blog post I present my approach and a step-by-step process I follow to calculate the burnt area from the Vilia wildfire. Not any rocket science here, this is a simple process, but very useful and effective to detect the affected area and, most important, it is scientifically correct, since it is based on the official Sentinel-2 products.
Before reading this post, I would recommend you firstly read the How to Make Outstanding Maps with Sentinel-2 and ArcGIS Pro – Part 1: Band Combinations, which I wrote a couple of months ago. It is also worth to mention that here I applied a new trick I recently learned in the Imagery in Action MOOC by Esri, with which I calculate the burnt area directly from a raster layer, without having to convert it to a vector polygon layer!
So, let’s dive into the process.
Preparing the Data
I create a new map in ArcGIS Pro, I specify its Coordinate System to Greek Grid (epsg:2100) and I navigate to my area of interest, which is the burnt area of Vilia at the Western Attica Regional Unit.
If you want to follow along with the same area, you can set the map extent to Top 38.1896114°N, Left 23.2267564°E, Right 23.4778065°E, Bottom 38.0501505°N and then zoom to it.
Next, I open the Catalog pane and in the Living Atlas tab I search for the Sentinel-2 Level-2A multispectral and multitemporal atmospherically corrected imagery, which I add to the current map. The imagery layer is added on the Contents pane and on the map, where I can clearly see the burn scar of Vilia (Picture 1).
I click on the Sentinel-2 imagery layer and I go to the Image Display Order group at the Data contextual tab. I click on the Sort button and from the dropdown list I select the last option to display the Most Recent image.
Then, at the Processing group I click on the Processing Templates button and from the dropdown list I select the Short-wave Infrared with DRA option to change its band combination. I finally zoom-out to see a wider area around Vilia and I add on the map the Sentinel-2 Tiling Grid (Picture 2).
My area of study falls in between the 34SFH and the 34SGH tiles of the Sentinel-2 Grid. So, I go to the Copernicus Open Access Hub I navigate to that area and I search for the most recent imagery captured in August 28, 2021. The platform returns the two products, the “after” the fire tiles, which I download and save in my Project’s folder (Picture 3).
I then perform a second query which includes the same two tiles, but for a former date, before the fire, which is the August 16, 2021. I also download and save the “before” tiles in my Project’s folder.
To make your life easier, in case you want to follow along with this case study, here are the download links of the products I use (login required):
Before the fire:
After the fire:
Create the Band Composites
Now, I have to create two band composites. One made out of the Sentinel-2 imagery captured before the fire and one after the fire. And for each one of them I have to mosaic the various bands, since the area of study falls in between two different Sentinel-2 tiles.
To do this, I click on the Function Editor button, at the Analysis group of the Imagery tab, to open a new Raster Function Template (see Picture 4).
I will firstly create the band composite for the after-the-fire image, captured in August 28. In the Raster Function Template I add the Blue Band (B02) from the 34SFH tile, as well as the Blue Band (B02) from the 34SGH tile and I chain them with a Mosaic Rasters function, which I insert from the Raster Functions pane (see Picture 4).
At the Mosaic Rasters Properties I rename the new mosaic to B02 and for Output Pixel Type I select 16 Bit Unsigned from the dropdown list (Picture 4).
I repeat the same process for the Green (B03), Red (B04), NIR (B08), Narrow NIR (B8A), SWIR-1 (B11) and SWIR-1 (B12) bands from the two Sentinel tiles and I rename respectively the output mosaics (see Picture 5).
I then drop a Composite Bands function, from the Raster Functions pane, and I chain all seven output mosaics with it, in the Function Editor. The seven mosaics are added on the Rasters list of the Composite Bands Properties, where I sort them with the order B04, B03, B02, B08, B8A, B11, B12 (see Picture 5).
It is also important to select Min Of from the dropdown list at the Cellsize Type option, to ensure that the output band composite will have the smaller possible cell size, hence the higher resolution, which in Sentinel-2 is 10 meters (Picture 5) .
I now click on the Save As button at the menu of the Raster Function Template to save is a as a new function in my Project. In the Save As dialog that appears, I write Band_Composite_2021_08_28 at the Name field, I select Project at the Category dropdown list, for Description I write Band Composite Aug 28, 2021 and I click OK (Picture 6).
If I open the Raster Functions pane and go to the Project tab, I will see my newly created and saved Raster Function, named Band_Composite_2021_08_28 (Picture 7).
I click on the Band_Composite_2021_08_28 function to open it and I click on the Create new layer button (Picture 8).
The raster function collects all seven bands from the two Sentinel-2 tiles, stored in my Project’s folder. It stitches the bands from the two tiles in seven new mosaic outputs and then combines the outputs in a new band composite, which is being added on the Contents pane and on the map (Picture 9).
Half of the work is done. I now have to repeat the exact same process but for the bands collected in Aug 16, before the fire. So, following the same steps, I create a new Raster Function, which I name to Band_Composite_2021_08_16.
The Band_Composite_2021_08_16 Raster Function is added on the Project tab of the Raster Functions pane. (Picture 10).
Running the Band_Composite_2021_08_16 Raster Function results with the Band_Composite_2021_08_16 layer added on the Contents pane and on the map (Picture 11).
I just created two new multiband raster layers. One that shows my study area before the fire (Band_Composite_2021_08_16) and on that shows the same area after the fire (Band_Composite_2021_08_28). The damage is clearly visible when I compare these two layers.
Calculating the Normalized Burn Ratio Index (NBRI)
The Normalized Burn Ratio Index (NBRI) uses the NIR and SWIR bands to emphasize burned areas, while reducing illumination and atmospheric effects. Its formula is simple and goes like
NBR = (NIR – SWIR) / (NIR + SWIR)
NIR = The pixel values from the Near Infrared band, which in Sentinel-2 is band 8 (B08) and
SWIR = The pixel values from the short-wave infrared band, which in Sentinel-2 is band 11 (B11) or band 12 (B12).
I will calculate the NBR for both raster layers, I created at the previous steps, because I want to have the NBR before and after the fire.
For this, I will firstly clip them to the current map extent to ensure faster processing. So, I open the Clip function where for Raster I select the Band_Composite_2021_08_28, for Clipping Type I select Outside and for Clipping Geometry / Raster I select the Active-Map-View-Extent option from the dropdown list (Picture 12).
I then click on the Create new layer button and the clipped version of the Band_Composite_2021_08_28 is added on the Contents pane and on the map (Picture 13).
To calculate the NBR for the clipped after-the-fire Band Composite I select it on the Contents pane and I click on the Indices Gallery button at the Tools group of the Imagery tab on the ribbon. At the list of spectral indices I scroll down to the Landscape group from where I select the NBR (Picture 14).
The NBR dialog appears in the middle of my screen, where I have to define which band of my Band Composite will represent the NIR band at the index calculation and which one will represent the SWIR band. I select Band_4 for the NIR band and Band_6 for the SWIR band (Picture 15).
Note that the Band_4 and Band_6 I selected have been defined at the step where I created the Band Composite in the Function Editor.
I click OK and after a couple of seconds a new raster layer is added on the Contents pane and on the map. This is the NBR for the after-the-fire Band Composite and by default it is rendered in grayscale.
To better visualize it, I select it on the Contents pane and I open the Symbology. For Primary symbology I select Classify, for Method I select Natural Breaks (Jenks) and for Classes I select 5 (Picture 16).
Natural Breaks (Jenks) do a very good job in separating the classes in a meaningful way, so the burn scar is clearly distinguished. I also use the Swipe tool to compare the NBR with the Natural Color Band Composite (Picture 17).
I repeat the exact same steps for the before-the-fire Band Composite, namely Band_Composite_2021_08_16. Briefly, I clip it in the same map extent and I calculate the NBR from the Indices Gallery (Picture 18).
Calculating the Burnt Area
In order to calculate the burnt area I have to create a change detection layer, derived from substracting the before-the fire NBR layer from the after-the-fire NBR layer.
From the Raster Functions pane I open the Calculator function. At the Raster Variables I create two variables, one named After and one named Before. I then assign the after-the-fire NBR layer to the After variable and the before-the-fire NBR layer to the Before variable (see Picture 19).
At the Expression field I type the simple calculation After – Before, which will subtract the Before NBR from the After NBR (Picture 19).
I click on the Create new layer button and a new difference layer is added on the Contents pane and on the map. I Symbolize it with Natural Breaks (Jenks) in two classes to distinguish the pixels representing the burnt area and all the other pixels (Picture 20).
The difference layer is a float raster, but in order to calculate the area I must have a unique values raster. For this I will use the Remap function to reclassify the pixel values of the difference layer.
Basically, I want to create a binary raster with only two pixel values, burnt and not burnt. I need a threshold for this, which has been already defined with the two classes Natural Breaks (Jenks). Of course, I can define a more precise threshold value with trial and error, but normally the NBR indicates burnt areas for pixel values below -0.25, so in my case I am OK.
I open the Remap function and for Raster I select the difference layer for Remap Definition Type I select List and at the list I enter the Minimum and Maximum values of the two classes of the difference layer, as defined by symbology (see Picture 21).
At the list, at the Output field, I enter the value 1 for the class representing the burnt area and the value 0 for the class representing the non burnt area, for which I also click on the NoData checkbox (Picture 21).
I click on the Create new layer button and the new reclassified raster is added on the Contents pane and on the map (see Picture 22).
Since I had selected NoData for the pixel values representing non burnt areas, the new raster has only one pixel value, the value 1 and all the rest has no data (Picture 22).
I now have to export the reclassified raster layer in order to create an attribute table for it. I select it on the Contents pane and from the Data contextual tab I click on the Export Raster button.
The Export Raster dialog appears where I write Burnt_Area.tif for the Output Raster Dataset, I select 16 bit Unsigned for Pixel Type and I write the value 0 at the NoData value field (Picture 23).
I click on the Export button and the reclassified Burnt_Area.tif raster is saved in my Project’s folder as a new raster file. If the saved raster is added one the Contents pane I will remove it in order to disable any lock that ArcGIS Pro might do, to ensure the next step will run properly.
I now open the Build Raster Attribute Table Geoprocessing tool and for Input Raster I navigate to my local folder and find the Burnt_Area.tif (see Picture 24).
The Build Raster Attribute Table reason will not run if the pixel type is floating point or double precision and this is the reason I selected 16 bit Unsigned for Pixel Type when exported the Burnt_Area.tif raster file (see Picture 23).
The Burnt_Area.tif is added on the Content pane and on the map and if I right-click on it I can see the Attribute Table option activated. So, I open its Attribute Table, where I can read the count of pixels with the value 1, which represent the burnt area. In my case the count is 825,547 pixels (Picture 25).
Now, I know that every pixel on my raster has a cell size of 10 meters. This means that every pixel has an area of 10 × 10 = 100 square meters.
So, the total of the pixels with value 1 have an area of 825,547 × 100 = 82,554,700 square meters, which equals to 82.5547 square kilometers!
Finally, I can use the Burnt_Area.tif with appropriate symbology to highlight the burnt area on the map (Picture 26).
Another way (the classic approach) to calculate the burnt area is to convert the Burnt_Area.tif raster file to a vector polygon file and then derive the Shape Area from that polygon. This can be easily done with the Raster to Polygon Geoprocessing Tool.
I performed this process to all the burnt areas in Greece, that were damaged in August 2021 and I calculated the burn scar for all of them. I then created a simplified vector polygon layer, with which I created the Wildfires in Greece 2021 Google Map. This was my little contribution to understanding the magnitude and severity of this year’s wildfires and to give common people the opportunity to realize the spatial distribution and shape area of the burnt areas.
I hope this blog post will help – or at least inspire – other geospatial professionals around the World to monitor fire events or other physical destruction events (such as floods etc) that may happen where they live. I would also be very happy to receive any feedback or suggestions for improvement.
Kindest regards from Crete, Greece
Your support is welcome!
If you really enjoy my blog and have benefited from the knowledge and experience I share here, you are more than welcome to make a donation or become a Patron, to help keep this content open, free and up-to-date!