A GIS Plugin for Susceptibility Modeling: Case Study of Wildfires in Vila Nova de Foz Côa

Abstract

Risk mapping is a crucial part of spatial planning, as it optimizes the allocation of resources in its management. It is, therefore, of great interest to build tools that enhance its production. This work focuses on the implementation of a susceptibility model for different types of spatially distributed risk in a geographic information systems (GIS) Python plugin. As an example, the susceptibility model was applied to study the occurrence of wildfires in the municipality of Vila Nova de Foz Côa, Portugal. The plugin was developed to simplify the production and evaluation of susceptibility maps regarding the available geographical information. Regarding our case study, the data used corresponds to three training areas, ten years of burned areas and nine environmental variables. The model is applied to different combinations of these factors. The validation, performed with receiver operating characteristic (ROC) curves, resulted in an area under the curve (AUC) of 74% for a fire susceptibility model, calculated with the same environmental factors used in official Portuguese cartography (land use and slope) and with the optimal training area, years of information on burned area and level of land use classification. After experimenting with four variable combinations, a maximum AUC of 77% was achieved. This study confirms the suitability of the variables chosen for the production of official fire susceptibility models but leaves out the comparison between the official methodology and the methodology proposed in this work.

1. Introduction

Wildfires constitute a phenomenon of crucial relevance in many countries. Due to a predominantly Mediterranean climate, a large number of forest fires occurs in Portugal mainland during the summer season, with a strong tendency to ravage shrubland [1,2,3,4]. With increasing rural abandonment, large-scale forestation programs and other climatic factors, the frequency and intensity of these wildfires is rising [5,6], with Portugal being one of the countries with the highest density of ignition and burned area [2,7]. The occurrence of fire events is not, as it is frequently assumed, a strictly negative phenomenon, as it can be essential for the regeneration of ecosystems classified as sensitive or even dependent on fire [8]. However, the inadequacy of spatial planning and anthropic bad practices may lead to wildfires that are prejudicial to both natural ecosystems and human civilization. In fact, most of the known causes of wildfire in northern Portugal involve intentional tort or negligence [9], and these events may lead to catastrophes, such as the one experienced by Portugal in 2017, in which Portugal was affected by two major forest fire events that occurred outside the typical forest fire season [10].

Catastrophic wildfires demand risk mapping. Here, it is important to distinguish “risk cartography” (a generalist term) from the terms “susceptibility map” (representing the propensity of an area to be affected by a given event in function of its characteristics), “hazard map” (representing the product between susceptibility and probability of occurrence, as conditioned by previous events) and “risk map” (representing specifically the product between hazard and potential damage) [11,12].

Geographical information systems (GIS) are an essential tool to generate fire risk cartography [13,14,15,16], allowing the quick and accurate analysis and combination of data from multiple sources, the manipulation of the resulting geoinformation and the generation of new data. However, the generation of complex spatial models implies the use of GIS tools that require considerable time (when using the tools one by one) and knowledge about the respective software. It would be, therefore, of great interest to have applications that facilitate the production of such maps, while still using GIS software.

In that context, some applications have been developed to facilitate the creation of fire risk cartography. For instance, Baranovskiy and Yankovich (2018) [17] created an embedded GIS software tool (under ArcGIS software) for forecasting, monitoring, and evaluating forest fire occurrence probability in Iran, using Python language. Mahmud et al. (2009) [13] developed an extensive Avenue programming script to deliver the fire vulnerability mapping in Malaysia, while allowing authorized users to edit, add or modify parameters whenever necessary, supporting fire hazard mapping using ArcView software. Bonazountas et al. (2007) [18] developed an integrated computer system based on semi-automatic satellite image processing (for fuel maps creation), socio-economic risk modelling and probabilistic models for forest fire prevention, planning and management in the island of Evoia, Greece. Gulluce and Celik (2020) [19] proposed a new fire detection method and monitoring software, FireAnalyst, for an early warning fire detection system aimed at valuable forested areas in Turkey, using the libraries of Google Maps’ application programming interface (API) in a cloud. Volokitina et al. (2021) [20] developed a fire simulation software to identify inventory plots ready to burn as well as to spread the rate for fire parts dependent upon weather conditions, predict fire intensity and fire development and calculate the required manpower and resources for fire suppression in Kazakhstan Altai.

In Portugal, a GIS open-source application [15] was already developed for the generation of risk cartography according to the specifications of the Portuguese authorities [21], which define susceptibility maps as the product between slope maps and the Corine land cover (CLC), according to the tabled values defined as a priori. Nevertheless, one can conceive of a methodology that can complement those environmental variables with others as well as using fire favorability scores that are adapted to local conditions. Based on the literature consulted, there is no other GIS application that provides intuitive tools to manipulate and generate fire risk cartography, which underlines the novelty of the proposed methodology for that purpose. Although there are other tools that are able to generate susceptibility models, to our knowledge, none of them have the flexibility to experiment with different training areas, years of occurrence and environmental variables that our tool displays.

The main objective of this work is to implement a susceptibility model considering environmental variables, that can be applied to different types of natural phenomena in a GIS open-source application under QGIS software, using Python language, which is able to: (i) be applied to different types of spatially distributed risk (such as fire occurrences, landslides or other); (ii) generate a susceptibility model, for a given study area, by calculating the susceptibility scores associated with multiple environmental variables; (iii) evaluate the training and prevision models generated, calculating the area under the curve (AUC) associated to the respective regressive operational characteristic (ROC) curves [22] and (iv) optimize the model generated by selecting different training areas, years of occurrence information and environmental variables used in the model.

2. Materials and Methods

The methodology proposed in this work comprises the development of a GIS plugin in an open-source environment. Figure 1 presents the workflow of the developed GIS plugin. The processes are further discussed in this chapter, as well as the preceding data acquisition and treatment.

2.1. Study Area

The study area corresponds to the municipality of Vila Nova de Foz Côa, located on the southern bank of the Douro river (Figure 2), being an example of a territory that should have a natural relationship to wildfires, due to its strongly Mediterranean climate, yet suffered great loss in 2017.

2.2. Dataset

Geographical information relative to the study and training areas (Figure 2) were extracted from the national geographic information system (sistema nacional de informação geográfica; SNIG) search engine [23].

Data on burned areas, considering 9 years (2009–2018) of occurrences in Portugal, was obtained from institute for the conservation of nature and forests (Instituto da Conservação da Natureza e das Florestas; ICNF) geocatalog [24]. That information was converted from vectorial to raster format, attributing the value 1 to pixels burned in this period and the value 0 to the remaining pixels. Although acknowledging that not all wildfires are equal, as larger fires are more associated with extreme weather conditions and denser fuel [25,26], this work does not make a distinction between burned areas.

This work sought to include the environmental variables, which were more plausible of contributing to the prediction of a territory’s fire susceptibility, represented in Figure 3. Susceptibility models created with the official Portuguese specifications [21] uses land use and terrain slope, as do most works on fire modeling [11,27,28,29,30]. Terrain aspect is also used often [11,26,27]. Verde and Zêzere [11,26] also explore mean average precipitation (which in Portugal seems to have a positive correlation with fire when considering fuel limited fire regimes, favoring the growth of fuel) and the number of days with minimum temperatures equal or superior to 20 °C (assuming that night time suppression efforts are more likely to succeed, thanks to lower temperatures and higher air humidity). Adab (2013) [27] also suggests moisture indexes and distance to roads and urbanizations. The normalized difference vegetation index (NDVI) is mostly used for burned area detection but can also be used to access susceptibility (as it expresses fuel availability) recurrently used to study fire risk [28,31,32,33,34] and is explored in this work considering both winter and summer values. As this is a work intended to be of wider application that can be used in several and different contexts, it does not include variables, such as elevation, topographic roughness index, annual average and maximum temperature, demographic variables, fuel connectivity, fire recurrence patterns and daily climate data, mentioned in other works [2,25,29].

Land use data were obtained, in vectorial format, from Portugal’s 2018 land use map (Carta de Ocupação do Solo; COS2018) [26]. These data were converted to raster assigning different integer values to each class, for each of the 4 hierarchical levels of classification available. This has resulted in four land use images, from the most general classification to the most specific.

A map of distances to artificialized territories, representative of roads and urbanizations, was obtained from the COS2018. That map was reclassified according to

Table 1. Classification of distance to artificialized territories.

Distance Interval (m)	Reclassified Value
Between 0 and 100	1
Between 100 and 200	2
Between 200 and 400	3
Between 400 and 800	4
Between 800 and 1600	5
Over 1600	6

.

A digital elevation model (DEM) with a 30 m of spatial resolution was obtained from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) satellite data, available on the United States Geological Survey (USGS) [35]. The slope and aspect maps were generated from the DEM (using slope native algorithm from QGIS software that considers the Horn’s method by default [36]), and both were reclassified according to

Table 2. Terrain slope classification.

Slope Interval (°)	Reclassified Value
Between 0 and 4	1
Between 4 and 8	2
Between 8 and 16	3
Between 16 and 32	4
Between 32 and 90	5

and

Table 3. Terrain aspect classification.

Aspect Interval (° North)	Reclassified Value
Plain (999)	0
Between 315 and 360	1
Between 0 and 45	1
Between 45 and 135	2
Between 135 and 225	3
Between 225 and 315	4

, respectively.

Satellite images with a 10 m of spatial resolution, obtained from the Multispectral Instrument (MSI) sensor onboard the Sentinel-2 were downloaded in the Copernicus European Space Agency portal [37]. Bands 8 and 4 of these images, respectively, correspondent to near infrared (NIR) and RED bands in the electromagnetic spectrum, were used to estimate NDVI for June and December 2018 (climatic extremes) according to Equation (1). That information was reclassified according to

Table 4. NDVI and GVMI classification.

NDVI/GVMI Interval	Reclassified Value
Between −1.0 and −0.3	1
Between −0.3 and −0.1	2
Between −0.1 and 0.1	3
Between 0.1 and 0.3	4
Between 0.3 and 1.0	5

. Bands 8 and 11, respectively, correspondent to NIR and short-wave infrared (SWIR), were used to calculate the Global Vegetation Moisture Index (GVMI) for June 2018, according to Equation (2) [38]. That information was also reclassified according to Table 4.

$$NDVI=\frac{NIR-RED}{NIR+RED}$$

(1)

$$GVMI=\frac{\left(NIR+0.1\right)-\left(SWIR+0.02\right)}{\left(NIR+0.1\right)+\left(SWIR+0.02\right)}$$

(2)

Data on total precipitation (mean annual values, in mm, for the period 1931–1960) was obtained from the national environmental information system (Sistema Nacional de Informação de Ambiente, SNIAmb) [39]. That information was converted to raster with the correspondence presented in

Table 5. Precipitation classification.

Precipitation Interval (mm)	Reclassified Value
Between 0 and 400	1
Between 400 and 500	2
Between 500 and 600	3
Between 600 and 700	4
Between 700 and 800	5
Between 800 and 1000	6
Between 1000 and 1200	7
Between 1200 and 1400	8

.

Hourly data on the temperature registered by the meteorological station (with relevant data) nearest to the study area, about 29 km northwest (in Folgares), was collected from the national hydrological resources information system (Sistema Nacional de Informação de Recursos Hídricos, SNIRH) website for 2015–2021 [40].

A Python script was developed to automatically generate a .csv file with the daily minimum temperatures registered in 2018. These data were associated, in vectorial format, to the point corresponding to the station’s location. Afterward, with the DEM previously used and Montgomery K.’s rule of thumb for the variation of temperature with altitude (−1 °C for every 200 m increase) [41], 365 raster files with minimum temperatures were generated for each day of 2018 in Vila Nova de Foz Côa (Equation (3)).

$$T_{raster}=T_{measured}-\Delta T_{sea}+\Delta T_{alt}$$

(3)

where T_raster represents the value obtained in each pixel (for each of the 365 raster files), T_measured represents the minimum temperature measured by the station on each day, ∆T_sea represent the difference in temperature expected between the average sea level and the station (which is at a 744 m altitude) and ∆T_alt represents a difference in temperature expected between the average sea level and any specific point, whose altitude was extracted from the DEM.

From these images, 365 other raster files were generated, reclassified with the value 0 in pixels where T <= 20 °C and with the value 1 in pixels where T > 20 °C. The sum of the values in those images represents the number of days with a minimum temperature equal or superior to 20 °C in each pixel of the study area. The maps of number of days with minimum temperatures equal or superior to 20 °C contained only four values (3, 4, 5, 6 and 7 days). Therefore, it was not necessary to define intervals (as it would be for a continuous scale, for example).

2.3. Implementation of the GIS Plugin

The GIS open-source application was developed under QGIS software, Version 3.16 [42]. Python programming language and several APIs and libraries were used. The “Plugin Builder” plugin was used to generate a plugin template, while the Qt Designer framework was used to stylize its graphical user interface (GUI), as seen in Figure 4. The developed GIS application is composed by a main scroll area where the user can browse for the input files—shapefiles corresponding to training and modelling areas, and raster files corresponding to a null variable (a raster image where all pixels in the extent of the study and training areas have the value “0”), at least 2 environmental variables (up to 4) and at least 2 occurrence years (up to 10)—and an area where the user can browse for the directory where the output files will be created. The “Run” button assigns those paths to string variables and passes them into the scripts described in Section 2.3.1 and Section 2.3.2, which create a susceptibility model and provide data for its evaluation.

The plugin can be downloaded at the GitHub repository linked as “Supplementary Materials” as well as a folder with some geographic information as examples. To execute the plugin on a QGIS software, the instructions in the PDF included in the folder should be followed.

2.3.1. Main Script

This work centers on the development of a Python script, named main.py, that uses the QGIS interface and processing tools to automatically generate susceptibility models. Besides importing the relevant libraries (“os” [43] and “pandas” [44]), all of the script’s content involves the definition of the “main” function, which was adapted to incorporate a plugin for QGIS (see Section 2.3). The function takes, as parameters, the instance of the plugin’s class (self), the “ROC” function (see Section 2.3.2.), the input file paths and the output folder directory.

The script starts by creating the folder where the resulting information will be saved if it does not yet exist. In this phase, it also defines the “clip” function (used recurrently in the script) and accommodates the variables that correspond to undefined occurrences and environmental variables (depending on whether or not the fields marked as optional are filled), adapting the expressions used in following calculations to those vacancies.

The script adds the geographic information corresponding to occurrences to the GIS project and calculates the sum of occurrences for each pixel (according to the expression initially defined). Afterward, it calculates a raster where the pixels with occurrences assume the value 1, and the pixels without assume the value 0 (boolean.tif) as well as a raster with the inverse logic (antiboolean.tif).

The information corresponding to the environmental variables and the occurrences (boolean and antiboolean calculated in the previous phase) is added to the project. That information is clipped by the defined training area, and the information relative to the area with occurrences is analyzed.

To calculate the susceptibility of a given area to occurrences, the informative value method was adopted. That method is usually used for analyzing the susceptibility of a territory to landslides [45,46] but also to fire [47] and other types of risk, including non-environmental [48]. The raster layer zonal statistics tool, used in the analysis mentioned in the previous paragraph, creates a table that relates the different classes of the reclassified variable (returned as “zone”) to statistics regarding the area with occurrences within such class, namely “count” (number of pixels belonging to the class) and “sum” (sum of the pixel values in the class). The function “calcular_scores”, which applies the informative value method, calculates favorability scores for each class based on the obtained statistics (Equation (4)).

$$ln\left(\frac{burned pixels in the zone}{total pixels in the zone}/\frac{burned pixels in the study area}{total pixels in the study area}\right)$$

(4)

After the calcular_scores function is applied to every environmental variable, the susceptibility map is generated according to an expression that sums the score values assigned to the different pixels of the specified variables, defined initially. A reclassified map is also generated, with the values “0” changed to “0.1” for them not to interfere with the validation described in Section 2.3.2.

From that second susceptibility map, the script extracts the pixels corresponding to susceptibility in the area with occurrences (true_part.tif) and in the area without (false_part.tif), which were analyzed with the rasterlayeruniquevaluesreport tool. The resulting tables are used in the ROC function, defined in Section 2.3, to generate the first validation graphics [22].

Lastly, the script uses the calculated favorability scores to generate reclassification tables for each variable, corresponding a given score to each class. The raster files corresponding to environmental variables are then clipped by the study area and reclassified according to the scores presented in the reclassification tables. Finally, the susceptibility map corresponding to the predictive model and the associated validation graphics are generated similarly to what was described for the training model.

2.3.2. Validation Script

As discussed in Section 2.3.1., the main script generates two tables, each describing the susceptibility of each pixel to a given type of risk in the areas with occurrences (positive) or without occurrences (negative) in the study area, respectively. The information in those tables is essential for the validation of the model with the “ROC” function, imported from the validacao.py scrip.

The “ROC” function uses the “csv” [49], “matplotlib” [50] and “numpy” [51] libraries. It starts by analyzing the information contained in the tables referred to in the previous paragraph, ignoring the 0 values (in the table with the negative pixel values, the “0” corresponded to positive values; in the table with the positive pixel values, the “0” corresponded to positive values) and afterward converting the value 0.1 into 0, correcting the reclassification made in the main script. With these data, it also returns statistics correspondent to the total number of positive pixels, the total number of negative pixels and the maximum and minimum scores for the negative pixels.

The first graphic simply plots two lines, corresponding to positive and negative occurrence pixels, which represent the frequency of each susceptibility value, as seen in Figure 5. As for the second graphic, the function defines 100 thresholds—imaginary lines marking the value that separates the pixels not susceptible to the risk (on the left) from the pixels susceptible to the risk (on the right)—for the interval between the maximum and minimum negative scores. For each, it determines the false positive ratio (1-specificity) and the true positive ratio (sensibility), according to Equations (5) and (6).

$$false positive ratio=\frac{negative pixels considered as positive}{total negative pixels}$$

(5)

$$true positive ratio=\frac{positive pixels considered as positive}{total negative pixels}$$

(6)

This generates 100 points of ordinates (y) defined by the true positive ratio and abscissae (x) defined by the false positive ratio, where the first point corresponds to y = 0 (no pixel is considered to be positive) and the last point corresponds to y≈1 (nearly all pixels are considered positive). The points are used to generate the ROC curve and, finally, the AUC is calculated using functions from NumPy library and the corresponding graphic is generated.

As an AUC around 0.5 represents a useless classifier (success indistinguishable from that of a random classifier), we can assume that an AUC over 0.6 is acceptable and over 0.7 is good [22]. Naturally, superior quality patterns may be used in other areas (in epidemiology, for example, a classifier with an AUC between 0.7 and 0.9 is considered to be only moderately good [52]). Nevertheless, and keeping the high complexity inherent to the theme of this study in mind, this work considers the first classification framework.

3. Results

3.1. Determination of the Ideal Training Area

To test the performance of the GIS application developed, the study area of Vila Nova de Foz Côa was used as an example. In a first phase, with the preliminary use of 6 years of burned areas and level 2 COS2018 land use classification, the developed application was applied for different training areas to evaluate which one returns the best wildfire susceptibility map for Vila Nova de Foz Côa in 2018, evaluating the AUC of the ROC curve associated to the predictive model. The municipalities west of the municipality in the same river margin (training area 1—TA1), the two municipalities in the northern margin of the river (training area 2—TA2), the four municipalities adjacent at northwest (training area 3—TA3) and the six municipalities adjacent with margins on the river (training area 4—TA4) were considered (Figure 1).

The training areas 2, 3 and 4 have produced a similar prediction AUC (

Table 6. Results of the validation of susceptibility models generated with 6 years of burned areas, level 2 COS2018 land use classification and different training areas.

Training Area	AUC (Training)	AUC (Prediction)
TA1	0.80	0.70
TA2	0.69	0.72
TA3	0.72	0.72
TA4	0.71	0.72

). Training area 3 was chosen as it generates the best AUC for the training model.

3.2. Determination of the Most Adequate Wildfire Occurrence Period

The methodology was applied, according to the premise of using TA3 and level 2 of the COS2018, for a different number of years of information on the burned area, in order to evaluate which of them allowed the generation of a better fire susceptibility map in Vila Nova de Foz Côa in 2018. Five options were considered: 2, 4, 6, 8 and 10 years of information (

Table 7. Results of the validation of susceptibility models generated with training area 3, level 2 of COS2018 land use classification and different years of information on wildfire occurrences.

Years	AUC (Training)	AUC (Prediction)
2	0.70	0.67
4	0.69	0.67
6	0.72	0.72
8	0.72	0.72
10	0.72	0.72

).

The models generated with 2 or 4 years had a lower performance than the others—albeit they were acceptable—perhaps because there was an extensive burned area in 2017, and because some quality is derived from the simple fact that artificialized territories and water bodies are largely excluded from the burned areas. In the remaining years, the quality of the model has stabilized, indicating that the additional information becomes redundant. This may be related to the changes in the land cover, which make the information about burned areas less relevant. We chose to select the 6 years of information, respecting Ocam’s razor principle.

3.3. Determination of the Environmental Variables to Use

To determine the environmental variables that generate the most adequate susceptibility map, the developed script was applied using the combination of TA3 and 6 years of information on burned areas.

3.3.1. Determination of the Ideal COS Level

To validate the susceptibility models with the different training areas, levels 1, 2, 3 and 4 of the COS2018 were considered (

Table 8. Results of the validation of susceptibility models generated with training area 3, 6 years of occurrences and different COS levels.

COS Level	AUC (Training)	AUC (Prediction)
1	0.70	0.70
2	0.72	0.72
3	0.73	0.73
4	0.73	0.5

).

Level 3 exhibited the best AUC for both the predictive model and the training model. Level 4 confirmed the tendency of an excessive class specification to sacrifice the quality of the predictive model for that of the training model (overfitting phenomenon), generating a prediction model equivalent to a perfectly random classification.

With that, level 3 was considered as the COS2018 analysis level that allows the generation of the best wildfire susceptibility model for Vila Nova de Foz Côa in 2018.

3.3.2. Determination of the Best Individual Environmental Variables

The developed application was executed, individually, for each variable. Land use clearly stands out as the best predictor of susceptibility to fire, while the slope is presented as the less meaningful predictor.

Table 9. Results of the validation of susceptibility models generated with training area 3, 6 years of occurrences and different individual environmental variables.

Variable	AUC (Training)	AUC (Prediction)
COS level 3	0.73	0.73
Distance to artificialized territories	0.56	0.57
Terrain slope	0.39	0.39
Terrain aspect	0.57	0.52
Winter NDVI	0.42	0.44
Summer NDVI	0.55	0.51
GVMI	0.44	0.52
Precipitation	0.60	0.35
Days with Tmin > 20 °C	0.56	0.46

presents the results of the model validation for each environmental variable, while

Table 10. Wildfire favorability scores generated by the developed script with training area 3, 6 years of occurrences and different individual environmental variables.

Variable	Class	Score
COS level 3	1.1.1 Continuous edified fabric	−3.88
	1.1.2 Discontinuous urban fabric	−5.65
	1.1.3 Empty spaces in built fabric	−1.36
	1.2.1 Industry	−1.05
	1.2.2 Trade	−1.2
	1.2.3 Agricultural facilities	−2.56
	1.3.1 Energy production infrastructure	−3.85
	1.3.2 Water infrastructure and waste treatment	−2.12
	1.4.1 Road and rail networks and associated spaces	−5.39
	1.5.1 Inert extraction areas	−1.74
	1.5.2 Waste deposition areas	−0.13
	1.5.3 Areas under construction	−3.51
	1.6.1 Sports equipment	−2.88
	1.6.2 Leisure facilities and campsites	−2.41
	1.6.3 Cultural equipment	−1.88
	1.6.4 Cemeteries	−1.27
	1.6.5 Other tourist equipment and facilities	−3.78
	1.7.1 Parks and gardens	−2.12
	2.1.1 Temporary rainfed and irrigated crops and rice fields	−0.5
	2.2.1 Vineyards	−3.31
	2.2.2 Orchards	−0.89
	2.2.3 Olive groves	−1.31
	2.3.1 Temporary crops and/or improved pasture associated with permanent crops	−4.38
	2.3.2 Complex cultural and partial mosaics	−1.2
	2.3.3 Agriculture with natural and semi-natural spaces	−0.36
	2.4.1 Protected agriculture and nurseries	−1.44
	3.1.1 Improved pastures	−0.05
	3.1.2 Spontaneous grazing	0.04
	4.1.1 Agroforestry	0.8
	5.1.1 Hardwood forests	−0.13
	5.1.2 Softwood forests	0.15
	6.1.1 Shrubland	0.67
	7.1.2 Bare rock	−1.73
	7.1.3 Sparse vegetation	1.17
	9.1.1 Water courses	0.41
	9.1.2 Water plains	−0.95
Distance to artificialized territories	Between 0 and 100 m	−2.26
	Between 100 and 200	−1.22
	Between 200 and 400	−0.73
	Between 400 and 800	−0.24
	Between 800 and 1600	0.14
	Over 1600	0.36
Terrain slope	Between 0 and 4°	−0.14
	Between 4 and 8°	−0.05
	Between 8 and 16°	0.09
	Between 16 and 32°	−0.01
	Between 32 and 90°	−0.04
Terrain aspect	Plain	−4.52
	North	−0.07
	East	0.27
	South	0.18
	West	−0.38
Winter NDVI	Between −1.0 and −0.3	−0.45
	Between −0.3 and −0.1	−0.58
	Between −0.1 and 0.1	−0.72
	Between 0.1 and 0.3	−0.54
	Between 0.3 and 1.0	0.16
Summer NDVI	Between −1.0 and −0.3	1.10
	Between −0.3 and −0.1	−0.69
	Between −0.1 and 0.1	−0.45
	Between 0.1 and 0.3	0.47
	Between 0.3 and 1.0	−0.09
GVMI	Between −1.0 and −0.3	0.37
	Between −0.3 and −0.1	0.60
	Between −0.1 and 0.1	0.22
	Between 0.1 and 0.3	0.53
	Between 0.3 and 1.0	0.40
Precipitation	Between 0 and 400 mm	−5.58
	Between 400 and 500 mm	−0.27
	Between 500 and 600 mm	−0.00
	Between 600 and 700 mm	−0.26
	Between 700 and 800 mm	−0.32
	Between 800 and 1000 mm	0.68
Days with Tmin > 20 °C	3	0.15
	4	0.19
	5	−0.20
	6	−0.83
	7	−2.63

presents the associated favorability scores.

3.3.3. Determination of the Best Combination of Environmental Variables

The susceptibility modeling plugin was executed for every pair of environmental variables with the collected data (

Table 11. Results of the validation of susceptibility models generated with different pairs of environmental variables. The value before the semicolon represents the training AUC, while the value after represents the prediction AUC.

	Level 3 COS	Distance to Artificialized Territories	Terrain Slope	Terrain Aspect	Winter NDVI	Summer NDVI	GVMI	Precipitation	Days with Tmin > 20 °C
Level 3 COS
Distance to artificialized territories	75; 75
Terrain slope	73; 74	60; 64
Terrain aspect	75; 74	65; 64	58; 55
Winter NDVI	73; 74	65; 67	56; 63	61; 59
Summer NDVI	75; 74	64; 63	56; 57	60; 53	61; 59
GVMI	75; 74	67; 67	61; 64	64; 62	64; 64	63; 61
Precipitation	75; 71	67; 59	61; 50	65; 51	63; 54	63; 47	68; 56
Days with Tmin > 20 °C	74; 72	67; 60	58; 50	62; 50	57; 50	61; 50	66; 59	64; 48

).

Of all combinations presented in Table 11, few have returned an AUC for the predictive model greater than or equal to the one obtained with the COS2018 variable individually. One can also verify that, indeed, the land use-slope combination (Figure 5 and Figure 6) returns one of the best predictive values. In Figure 6, the area of each category would be approximately 20% of the total area (39,804 ha), which represents 7960.8 ha each.

Furthermore, it was decided to cross the combinations, which have resulted in predictive models classified as good in groups of three—combining the level 3 COS variable with all others.

Table 12. Results of the modeling script combining level 3 COS with every other environmental variable in trios. The value before the semicolon represents the training AUC, while the value after represents the prediction AUC.

	Distance to Artificialized Territories	Terrain Slope	Terrain Aspect	Winter NDVI	Summer NDVI	GVMI	Precipitation	Days with Tmin > 20 °C
Level 3 COS + Distance (…)
Level 3 COS + Terrain slope	75; 76
Level 3 COS + Terrain aspect	76; 76	75; 74
Level 3 COS + NDVI Winter	75; 76	73; 76	75; 75
Level 3 COS + NDVI Summer	77; 76	75; 76	76; 75	75; 75
Level 3 COS + GVMI	77; 76	76; 75	76; 74	76; 75	76; 74
Level 3 COS + Precipitation	77; 74	75; 72	77; 72	75; 72	77; 72	77; 72
Level 3 COS + Days Tmin > 20	76; 74	74; 73	75; 72	73; 72	76; 73	77; 73	75; 71

shows how several of the prevision models generated using three variables exhibit a higher AUC than the one generated with the land use and slope variables, although it is a fairly small difference. Nevertheless, it was decided to cross the best trios (with AUC = 76%) with the variables involved in those crossings in combinations of four.

Table 13. Results of the modeling script combining groups of four environmental variables. The value before the semicolon represents the training AUC, while the value after represents the prediction AUC.

	Terrain Aspect	NDVI Winter	NDVI Summer	GVMI
Level 3 COS + Distance (…) + Terrain slope	76; 76	75; 77	77; 77	77; 76
Level 3 COS + Distance (…) + Terrain aspect
Level 3 COS + Distance (…) + NDVI Winter	76; 77
Level 3 COS + Distance (…) + NDVI Summer	78; 76	77; 77
Level 3 COS + Distance (…) + GVMI	78; 76	78; 77	78; 76
Level 3 COS + Terrain slope + NDVI Winter	75; 75			78; 77
Level 3 COS + Terrain slope + NDVI Summer	76; 75	75; 76		78; 76

presents the combinations of four variables mentioned before, where the highest AUC achieved for prediction models was 77%. This value resulted from several different combinations of variables.

4. Discussion

This work has resulted in a plugin that is able to generate susceptibility models that can be evaluated and optimized by the variation of the inputs, thus achieving the proposed objectives.

Different susceptibility models were generated by variating the shapefile used as the “Training area” input in the GUI (Figure 4). In theory, that shapefile could be a polygon with any limits insofar as it is fully covered by the raster data used as the respective occurrence inputs and environmental variable input. The same goes for the shapefile used as the “Modelation area”: this plugin can model the susceptibility of any area, as long as it is covered by both the occurrence and environmental data used.

Different susceptibility models were generated by variating the raster data just mentioned. The AUC obtained for the particular models evaluated in the results does not matter as much as the fact that it is possible to obtain such AUC and to save it, as well as the respective ROC curve graphic, for later comparison. No automatic decision algorithm was adopted for choosing the best models, mainly because the computational power required by QGIS software is substantial and programming the calculation of several models would make it excessively slow and prone to crashing, thus limiting its utility as a plugin.

The consideration of occurrence data had the logical presupposition of a single year of information being insufficient for building a good model, as well as that of decreasing utility of data with increasing age. Thus, the first model produced in that consideration involved the two earliest sets of information on risk occurrence, and the variation in the modeling was simply made by adding earlier data to the most recent data.

The consideration of environmental variables, on the contrary, involved the comparison of single variable models (achieved using a null raster as the “Variable 2 raster” input), and the variation in the number of variables used in the combination of the different variables in the double-entry Table 11, Table 12 and Table 13. In order to avoid excessive combinations in a work that is supposed to design and test the performance of a tool, the variables used in combinations of three and four were limited by imposing the presence of the best individual variable in the combinations (Level 3 COS) and later by imposing the presence of the best previous combinations.

As for the results of this particular case study of wildfires in Vila Nova de Foz Côa, several points may be highlighted.

In the analysis of the susceptibility of different land use classes to fire, the sparse vegetation registered the highest score, while spontaneous grazing, agroforestry, softwood forests and shrubland were also associated with positive susceptibility values. The classes associated with artificialized territories and agriculture were invariably associated with negative scores. Previous research corroborates the results in relation to the preference of wildfires for shrubland and for softwood over hardwood [53]. The water courses class (of COS) was surprisingly associated to high susceptibility, although it is safe to say it would be related to wildfires on the riverbank of the affluents of Douro, which are usually steep (inviting easier fire spread) and occupied by shrubland.

Susceptibility generally grows with the distance. The peak of susceptibility for slope seems to occur in the middle class (8° to 16°), although the difference between the highest and lowest scores is notably small. There were also no variations in susceptibility with respect to terrain aspect, except for the extremely low score associated with flat territory (probably related to large water surfaces and human occupation), corroborated by previous studies [54]. In the winter season, the NDVI values revealed a peak in the susceptibility score for high values (> 0.3)—although there were no considerable variations—which is likely due to areas with a lower NDVI being associated with artificialized territories or water bodies. In the summer season, the NDVI revealed a considerable peak in susceptibility score for low values (< −0.3), but otherwise do not seem to exhibit a pattern. The moisture index in summertime points out to greater susceptibility with lower values, as expected, since dry fuel is more prone to burning. Precipitation shows a peak susceptibility value in the highest class and a noticeable minimum in the lowest class, which may be related (as with NDVI) to the highest growth of vegetable fuel during winter, where most Mediterranean precipitation occurs. As for the number of annual days with a minimum temperature above 20 °C, susceptibility shows a clear tendency to decrease with an increasing number of such days, contrary to what would be expected (perhaps due to the method used for temperature estimation, strictly related to altitude).

All training areas considered, composed by different combinations of surrounding municipalities, proved to adequately constitute study areas, although the four adjacent municipalities in the northwest were determined to be the best by a small margin. The ideal wildfire occurrence period seems to involve 6 years of information, which may be the point where there is enough information to allow the generation of a good model and a further amount of information will not improve it significantly. Where environmental variables are concerned, land use seems to be best used with level 3 specification when using the Portuguese COS classification, as it was found to be the best predictor of wildfires by far. When combining two environmental variables to generate a susceptibility map, every model that involved land use returned a prevision AUC > 70%, while no model that did not involve it has done so. When combining three environmental variables, several models returned a prevision AUC > 75%. When combining four environmental variables, six combinations returned a peak prevision AUC of 77%. These results can be verified in the light of the scale discussed in Section 2.3.2, which would rate an AUC over 70% as good.

The use of AUC–ROC curves has been used to validate the performance of several models [55,56,57]; however, to our knowledge, the methodology presented in this work was not applied in other studies to generate a fire susceptibility model, including the development of a GIS plugin, that enforces the use of the methodology implemented. There are several susceptibility fire models using different techniques and methods, but not compared to this methodology. For instance, Hong et al. (2018) [55] used genetic algorithms (GA) to obtain the optimal combination of forest fire-related variables and apply data mining methods for constructing a forest fire susceptibility map in Dayu County (China), validating the model performance with AUC–ROC curves. Eskandari et al. (2021) [58] also predicted the variables to be used in a model using the random forest (RF) algorithm, in Golestan Province (Iran). The use of machine learning (ML) techniques has improved the efficiency of fire prediction [59,60,61]. Kalantar et al. (2020) [56] applied adaptive regression splines (MARS), support vector machine (SVM) and boosted regression tree (BRT) to estimate fire susceptibility in Chaloos Rood (Iran), and the results were also validated using AUC–ROC curves. Zhang et al. (2019) [62] used a deep learning algorithm, particularly the convolutional neural network (CNN), to estimate a spatial prediction model for forest fire susceptibility.

The main strength of this approach for modeling susceptibility is the flexibility it offers in adapting the models to different regions, applying the methodology to different risk types, and testing different variables in search of the best model. The main weakness would be the more simplistic approach presented for evaluating the resulting models, as it can leave out several possible scenarios and is still very time-consuming, although it was conceived to achieve quicker computational processes. There is also the threat of a lack of data, particularly in terms of the past occurrence of certain types of risk, limiting the application of this model. Furthermore, the high adaptability of this plugin means the opportunity for replacing older methodologies that use fixed environmental variables, where scores are defined at the national level, in the production of regional and local risk cartography. The results obtained in this study are mostly demonstrative of the methodology implemented in the GIS application. The selection of training areas would be different and based on their similarity both in the fire regime and in the biophysical characteristics. The presentation of a single performance indicator is not sufficiently informative, especially considering the spatial homogeneity of some of the variables used in this study. Many of the ecological processes or types of hazards are context dependent, so the training models to apply them in another context will not always give good results and this should be considered in future studies.

One of the main advantages of this work was to design a flexible tool with the potential to be applied in future works. The methodology implemented in the GIS application can be used for any study area, including several training areas, if the available dataset allows it, considering the same study area with past land use and older occurrences as a training area. More environmental variables, other than the ones experimented with, can be used in search of an optimal susceptibility model, considering that the different types of wildfires can be related to different variables, which means there should be different models built for susceptibility to small, shrubland fires and to larger, high fuel load fires. Moreover, the developed GIS plugin has a wide scope in both function and language, which means it can be applied in different contexts, other than wildfires.

Risk mapping is a crucial part of land management and land planning. The possibility of predicting which areas are susceptible to a specific type of disaster, including landslides or forest fires, is unquestioned. This work analyzed the susceptibility of different land use classes to fire. The GIS plugin developed in this study is an essential tool for land management and planning, which can be used by land use planners, foresters, wildfire risk analysts and policymakers, among others.

5. Conclusions

This work presents a new GIS application, free and open-source, for generating susceptibility models, running as a plugin under QGIS software. The developed application is capable of automatically generating susceptibility models and returning ROC curves with their respective AUC values, thus facilitating the selection of the best model for use in risk cartography.

The GIS application was used for the generation of several wildfire susceptibility models in the municipality of Vila Nova de Foz Côa, Portugal, exploring different training areas, occurrence periods and environmental variables. From the results obtained, it was possible to confirm the adequacy of the environmental variables adopted by the ICNF for the calculation of wildfire susceptibility, as the model obtained with the land use and slope pair returned a prevision AUC of 74%. Nevertheless, this work does not compare this model with the one that would be obtained with the methodology used by the ICNF.

In the future, it would be interesting to compare models obtained with this methodology with models obtained with the methodology supported by official documentation. Furthermore, as this tool was designed for application in a great variety of risk cartography, such comparison should not be limited to wildfire susceptibility. Such research could prove the developed application to be useful in various contexts.