Fire Risk Analysis using QGIS

NOTE: This is a work in progress, and subject to change.

Purpose

This tutorial will demonstrate how to use multiple tools within QGIS in order to run a fire risk analysis, manipulating vector and raster data, along with the final output of maps. This tutorial uses tools from QGIS' advanced interface, which include QGIS, GRASS GIS, and plugins that you will need internet connection for in order to download them into QGIS. Major tools that will be used in this tutorial include:

Buffer - QGIS, Build Virtual Raster Catalog - QGIS, Clipper - QGIS, Merge Vector Layers - QGIS, Polygonize - QGIS, Raster Calculator - QGIS

Aspect and Slope - GRASS GIS, Reclass - GRASS GIS

Heatmap - Plugin, along with other tools.

Although this tutorial is a walk through, it is highly recommended that you have knowledge of a GIS before attempting this analysis, this is an advanced tutorial.

Introduction

For the purpose of this tutorial, the study area that will be used to demonstrate a Fire Risk Analysis will be approximately 75% of Vancouver Island from just South of Telegraph Cove to Victoria (Figure 1). During this fiscal year (April 1, 2015 – March 31, 2016) there have been 1 843 fires in the entire province of British Columbia which has resulted in approximately 305 000 hectares of land being burn. The area that we will be analyzing for the purpose of this tutorial around Vancouver Island has had just under 200 fires resulting in approximately 25 000 hectares of burnt land according to the Government of B.C. Wildfire Service.

This tutorial will show how to create maps that identify areas that have a high risk of fires based on multiple factors. These maps in turn, can help responders react quicker to help prevent future fires by implementing measures based on the high risk of fires in a certain area on the island.

The way risk will be calculated evenly for the purpose of this analysis will be from 1 to 4. With 1 being the lowest possible risk, and 4 being the highest possible risk. Using reclassification methods, this will be accomplished.

1 = Low Risk, 2 = Moderate-Low Risk, 3 = Moderate-High Risk, 4 = High Risk.

All data and software used in this tutorial are open-sourced and can be found online.

Figure 1: The study area for this tutorial is shown here as a pink grid feature, from just south of Telegraph Cove to Victoria.

Data

Data Identification

In order for the analysis to be run, data sets were identified through based upon research of academic journals and research papers that topics included a fire risk analysis, assessment or use of a GIS with forest fires. Based upon the research, the 5 major data sets that were needed include:

1. Census subdivisions & attribute data (population density)

2. Landsat 8 imagery (Normalized Difference Vegetation Index - NDVI)

3. Digital Elevation Model (aspect and slope - DEM)

4. Past wildfires (spatial pattern/clustering)

5. Mitigation factors (bodies of water)

The data sets were then categorized into three different categories (Table 1), which will create three maps that will be calculated against each other; fire hazards, vulnerability factors & mitigation factors.

Table 1. Data sets that were collected, categorized into three different categories.

Data Set Links & Guide

1. Census Subdivisions Attribute Data Attribute Data *

Download the census subdivision boundary file as a shapefile and download the census subdivision attributes.

2. Landsat 8 Imagery **

Once you have created an account, under Search Criteria select the area of your interest by creating a polygon around your study area. Select an appropriate date range, in our analysis we looked at data sets from July 2014 to August of 2014 when forest fires are most prominent in this region. Next, click Data Sets navigate to Landsat Archive, check the first box under that folder titled L8 OLI/TIRS. Under Additional Criteria you may select the whether the imagery was taken during day or night and the percentage of cloud cover in the imagery. For the purpose of this tutorial, it would be ideal to select imagery during the day and with less than 10% cloud coverage. After all of those processed have completed your data sets should appear and can now select the appropriate imagery to download. Note, you will have to unzip the files after they have downloaded in order to access the bands.

3. Digital Elevation Model

Simply pan the globe to your study area, click the tile and download. The elevation model is at a 90-metre resolution. Note, you will have to unzip the files after they have downloaded in order to access the elevation models.

4. Past Fires

Download the KMZ files from 2010 - 2014 (three separate files). QGIS will not be able to open up a KMZ file, it essentially is a Zipped KML file. You will have to unzip the KMZ file in order to get it into a KML which QGIS will be able to read and import.

5. Water Bodies

Download the zipped folder containing multiple shapefiles of British Columbia, after you have unzipped the folder the shapefile that will be needed for the analysis will be called, waterways.shp.

* Must be a Carleton University alumni, faculty, student or other institution approved by CHASS in order to access this attribute data.

** In order to access the Level 1 data that includes the bands needed to calculate a NDVI, you must create a USGS account.

Note: It is recommended you download 7-Zip to help unzip the Landsat Imagery and KMZ files.

Tutorial

Setting up the Environments

There are three steps in order to set up your environment for this tutorial.

1. The first thing to do after opening the QGIS desktop, navigate to the Project tab and select Project Properties (Figure 2). This is where you can name you project (General Tab) and ensure your coordinate system is correct (CRS Tab). Under the CRS Tab, check 'Enable on the fly CRS transformation', and choose an appropriate coordinate reference system, we chose a WGS 84 (ESPG:4326).

Figure 2. Opening the properties tab to name your project and define an appropriate coordinate system for your study area, by enabling 'on the fly projection'. Here a WGS 84 coordinate system was used.

2. The next step will be to add the plugin you will need for an analysis. Navigate to the Plugin tab (ensure you have internet connection) and select Manage and Install Plugins. The plugin interface will now appear, in the search box located at the top type in 'heatmap' and check the box beside the heatmap plugin (Figure 3). This tool creates a heatmap (raster) based on the proximity of vector points to one another, detecting spatial pattern/clustering which will be used for the past fires data.

Figure 3. Setting up the heatmap plugin needed in order to create a heatmap for the past fires points.

3. The final step in setting up your environment is to turn on the adavnced interface for your toolbox in order to access the GRASS GIS tools. Navigate to the Processing tab and select Toolbox and the toolbox interface should appear, and at the bottom of the toolbox Simple Interface should be on, click that box and turn Advanced Interface on (Figure 4).

Figure 4. Setting up the advanced toolbox interface.

Importing the Data

Since we have multiple types of data to be imported, we will have to use the Add Vector Layer and Add Raster Layer , located along the left hand-side of the QGIS user interface.

The data to be imported using the Add Vector Layer includes the following; Census subdivision cartographic boundary shapefile, three past fires kml files, and the waterways shapefile.

The data to be imported using the Add Raster Layer includes the following; Three Digital Elevation Models tiff files, and the three Landsat 8 Imagery tiff files. In order to calculate the NDVI, we will need the Red (4) and Near Infrared (5) bands, they will both need to be selected per imagery (Figure 5).

Figure 5. Importing the Red and Near Infrared bands per imagery, in order to calculate NDVI.

Organizing the Data

In order to keep your data organized, and to avoid confusion it is a good idea to group your data based on the five data types that were imported. After importing the data it will look quite messy and may not make much sense, so it is recommended that you change the name of the data immediately after importing to a name that makes sense or relates to set data layer. In the layer on the left hand-side right click your data, click properties and rename your data for each - it is recommended you do this for each data set imported. (Figure 6)

Figure 6. The figures showing the messiness of data, along with renaming and the ideal layer renaming scheme.

The next step is to group the data based on the original 5 categories that the data were imported in, as depicted in Figure 7. Hold the control button on PC or command on Mac and select the appropriate data layers to group, right click and select group and rename accordingly.

Figure 7. The figures show how to group the data along with the ideal grouping scheme.

Working with the Data

Now we can begin to use tools to manipulate the data. We first need to work with the Census subdivision (CSD) boundary shapefile, and need to create a CSD of Vancouver Island. However, this will not be our area of interest. In order to create a new shapefile of just Vancouver Island you must turn the attribute toolbar on, it should be automatically turned on - if it is not right click just under the properties and check the attribute toolbar on. After you have done so, move to select features (it can be done by polygon, freehand etc. - whatever you feel most comfortable with) and select all census subdivisions that reside on Vancouver Island (Figure 8).

Figure 8. The figures show how to turn the attribute toolbar and select features in order to create the new shapefile.

After you select the method you wish to select features with, in the case shown with Figure 9 hit control or command while selecting the polygons. After all features are selected, right click the CSD layer and save the selected features as a new shapefile, give it an appropriate name.

Figure 9. The figures show Vancouver Island being selected and how to save the selected features as a new shapefile.

Creating an Area of Interest & Building Virtual Raster Catalogs

In order to create our area of interest we will need to build virtual rasters based on the Landsat imagery and clip it to the newly created Vancouver Island CSD shapefile. The first step is to turn all of your Landsat 8 bands on, then navigate to the Raster tab, click Miscellaneous and finally click Build a Virtual Raster (VRT). For all of the rasters imported you will need to create a virtual catalog for each group. Turn the group on that you wish to convert into one raster; you will need to do it for the Landsat 8 Imagery: Band 4, Landsat 8 Imagery: Band 5 and the three Digital Elevation Models. While creating the virtual raster (VRT), ensure that no data is set to 0 and the resolution is set to highest, these rules can be applied to all of the above (Figure 10).

Figure 10. These figures show how to create a Virtual Raster (VRT), and the rules can be applied to all of the raster data sets.

After the Virtual Rasters are created, you will want to clip it to just Vancouver Island. For this we will used the Band 4 VRT, navigate to the Raster tab, click Extraction and use the Clipper tool, input Band 4 as the raster to be clipped and use the Vancouver Island as the mask layer, ensure that no data is set to 0 (Figure 11).

Figure 11. The figures showing how to use the clipper tool, using band 4 as the raster, Vancouver Island as the mask layer and setting 0 as no data and the newly created band 4.

After using the clipper tool and creating the new band 4 raster, you will want to reclassify all the values to 1 in order to create the area of interest. Navigate to Raster and click Raster Calculator and all values that are equal and/or greater than 0 shall equal 1. After this has been created you will want to create a polygon based on this value which will be your area of interest/study area, by using the polygonize function that is located under Raster and Conversion (Figure 12).

Figure 12. These figures show how to use the raster calculator to create a simple reclassification (all values that are equal and/or greater than 0 equal 1), which was then converted to a polygon for the area of interest.

After successfully create your area of interest/study area. You can again apply the clipper function to the other two VRT created using the area of interest shapefile as your mask (Figure 13).

Figure 13. The figures above show Band 5 (top) and the Digital Elevation Model (bottom) clipped to the study area.

Calculating Normalized Difference Vegetation Index and Reclassification

For this part of the tutorial, we will be using the raster calculator which was touched upon in the previous section in order to create a Normalized Difference Vegetation Index (NDVI), along with reclassifying where we first see the GRASS GIS tools being applied. The data that will be used in this section will be the newly clipped Band 4 (Red) and newly clipped Band 5 (NIR). A NDVI calculation requires the following: NDVI = NIR - RED / NIR + RED.

Using the raster calculator we will be able to do this calculation, the output values could/should range from -1.0 to 1.0. The reclassification tool from GRASS GIS has a difficult time reclassifying decimals and usually rounds up. So in order to reclassify more accurately, you will have to multiple the NDVI created by an appropriate value and then reclassify. We chose to multiply it by 100 (Figure 14).

Figure 14. The original NDVI calculation, the NDVI calculation multiplied by 100 in order for reclassification and the output values & visualization map.

Based on the values outputted by the last raster calculation, we now can reclassify based upon those values. Justification for the reclassification is as follows;

please refer to introduction for definition of risk numbers #Introduction

35 to 90 = 4

25 to 35 = 3

10 to 25 = 2

0 to 10 = 1

Now navigate to the advanced toolbox interface that we set up in section 5.1, type in r.reclass and open the tool. Choose the appropriate raster for reclassification, in this case the raster that was multiplied by 100. You will notice that it requires a text rule, simply open any word processor - Notepad is always quite easy to work with. Open the word processor and type in the justification above instead of to use thru, save it as a .txt file, input it into the tool and run the process and you should have reclassified your NDVI which will await further calculations (Figure 15).

Figure 15. The above figures show where to find the reclassification tool, how to reclassify using a text file along with using the tool and finally the output from the reclassification.

Aspect, Slope and Reclassification

In this section we will be using the DEM that we clipped to our area of interest to preform aspect and slope analyses. Aspect is the compass direction slope is facing, and output generated usually represents degrees from 0 to 360 and will help us determine which way the wind traverses on land, where Slope is incline/steepness of the surface and the output generated usually represents degrees from 0 - Flat to 90 - Straight up, The steepness of the slope affects both the rate and direction of the fire spread. Fires usually move faster uphill than downhill and the steeper the slope, the faster the fire will move. (1)

To create the Aspect and Slope rasters we will use the grass tool called r.slope.aspect which can be found in the toolbox (advanced interface), after you have opened the tool input your DEM that was clipped to your area of interest, and give appropriate names to your aspect and slope models. If you scroll further down on the tool, you will notice it outputs many other types of terrain analyses - however they will not be needed for our overall analyses, so turn all the boxes off except for aspect and slope, run the tool (Figure 16).

Figure 16. Type in r.slope.aspect in the advanced interface toolbox, and input all of the necessary components in order to run the analyses.

The next step will be to reclassify the aspect and slope rasters just created. Since Westerly winds prevail over Vancouver island (2) aspect that faces east have the highest risk for this analysis, the justification for slope is as follows:

0 to 23 = 2 (North)

23 to 68 = 3 (Northeast)

68 to 113 = 4 (East)

113 to 158 = 3 (Southeast)

158 to 203 = 2 (South)

203 to 248 = 1 (Southwest)

248 to 293 = 1 (West)

293 to 338 = 1 (Northwest)

338 to 360 = 2 (North)

Now use the same methods for slope, slope in degrees with greater than or equal to 35 has high risk, the justification for slope reclassification is as follows;

0 to 10 = 1

10 to 25 = 2

25 to 35 = 3

35 to 90 = 4

After both have been inputted into a .txt file, exactly like the processes in section 5.6 #Calculating_Normalized_Difference_Vegetation_Index_and_Reclassification, open r.reclass input your aspect raster and reclassification rule and run the tool, do the same for slope. The output maps should look similar to the figures in Figure 18.

Figure 18. Slope and Aspect Maps after reclassification and change of symbology.

In order to create the output that will be used in the hazard calculation, we will need to combine aspect and slope to identify areas that face east and have high slope values associated with it, we will use the Raster Calculator in order to complete the process. Multiply both rasters by each other, and since both reclassified rasters have maximum values of 4 the maximum output value of the new raster map (we called it wind) will be 16. This can then be reclassified by using r.reclass and you will have your aspect and slope output (Figure 19), the reclassification scheme follows;

1 to 4 = 1

4 to 8 = 2

8 to 12 = 3

12 to 16 = 4

Figure 19. Map for reclassified slope and aspect calculation, then the reclassification scheme for the newly created raster and the final output of the wind raster which will be used in the final calculation.