Utilizing the Available Open-Source Remotely Sensed Data in Assessing the Wildfire Ignition and Spread Capacities of Vegetated Surfaces in Romania

Abstract

We bring a practical and comprehensive GIS-based framework to utilize freely available remotely sensed datasets to assess wildfire ignition probability and spreading capacities of vegetated landscapes. The study area consists of the country-level scale of the Romanian territory, characterized by a diversity of vegetated landscapes threatened by climate change. We utilize the Wildfire Ignition Probability/Wildfire Spreading Capacity Index (WIPI/WSCI). WIPI/WSCI models rely on a multi-criteria data mining procedure assessing the study area’s social, environmental, geophysical, and fuel properties based on open access remotely sensed data. We utilized the Receiver Operating Characteristic (ROC) analysis to weigh each indexing criterion’s impact factor and assess the model’s overall sensitivity. Introducing ROC analysis at an earlier stage of the workflow elevated the final Area Under the Curve (AUC) of WIPI from 0.705 to 0.778 and WSCI from 0.586 to 0.802. The modeling results enable discussion on the vulnerability of protected areas and the exposure of man-made structures to wildfire risk. Our study shows that within the wildland–urban interface of Bucharest’s metropolitan area, there is a remarkable building stock of healthcare, residential and educational functions, which are significantly exposed and vulnerable to wildfire spreading risk.

1. Introduction

Climate change and global warming are expected to affect natural hazards like flooding and wildfires worldwide. These may have multiplied domino effect consequences on other natural and urban systems, leading to severe disasters at local scales [1,2,3]. While the emergence of flooding events relies mainly on the weather conditions and the natural/artificial properties of the catchment area, wildfires implicate human behavioral activities. The multifaceted character of the wildfire phenomena is acknowledged in the literature. Chapin et al. [4] defined wildfire as a wicked problem. According to Levin et al. [5], wildfires are multi-layered phenomena that implicate diverse interacting cycles between causes and effects acting in certain territories. Identifying the relevant factors that significantly correlate with the wildfire regimes remains a critical challenge to scientists [6,7,8].

In the classical wildfire assessment approach, the interaction of favorable hydrometeorological conditions with the study area’s geophysical and vegetation properties is considered the core prerequisite of the fire environment triangle, which consists of three pillars: weather, topography, and fuel [9]. Lightning strikes are the primary natural igniters [10]. However, most wildfires are reported to be caused by human activity, either intentionally or accidentally [11]. Human activity patterns have become a determinant during the wildfire ignition phase [12,13].

Mansuy et al. [14] contrast the anthropogenic factors to the macro-environmental ones and report that the human footprint affects almost equal wildfire risk inside and outside the North-American protected landscapes. The consequences of human activities on fire regimes are reported to overshadow the effects of climate change significantly [15]. The effect of societal habits like the Daylight-Saving Time (DST) alterations have been acknowledged to increase wildfire ignitions. For example, Kountouris [16] reports that DST transition during the Spring season has increased the number of non-prescribed wildfire ignitions by about 30% in the US, relying on around 2 million wildfire ignition of 23 years records.

A more recent study presents also the impact of COVID-19 lockdown on the wildfire regimes in a wildfire-prone region like the Mediterranean [3,17]. The authors report a significant decrease in the total burned area during this period compared to the estimations that counted for similar drought-related circumstances to previous years. The decrease in social activities has resulted in a significant reduction of wildfire events. Thus, the integration of anthropogenic factors within the wildfire risk assessment tools has become indispensable to increase the models’ sensitivity [18].

The impact of anthropogenic factors on wildfire regimes has gained considerable attention in the literature. However, their combined usage alongside with hydrometeorological and biophysical factors in wildfire spreading capacity models, is rarely present in the available published research. Unlike our previous studies [19,20], we introduced population density as a new criterion within the wildfire ignition probability/wildfire spreading capacity index (WIPI/WSCI) model, considering that the current literature tightly correlates to the population density and the wildfire ignition risk and frequency [21,22].

This study aims to develop a comprehensive and practical GIS-based model for assessing the wildfire ignition and spread capacities based on freely available and remotely sensed geospatial data. The proposed model is aimed to be reproducible to other vegetated surfaces where the remotely sensed data are available. Another goal is to test the utility of the Receiver Operating Characteristic/Area Under Curve (ROC/AUC) method in weighting each criterion’s impact factor by comparing with the Analytic Hierarchy Process (AHP), which is widely used in previous studies. Here we aim to deliver tangible graphical (maps) and statistical results about the wildfire ignition and spread capacities in Romania, supporting disaster risk reduction agendas nationwide.

2. Materials and Methods

2.1. Study Area

The study area is represented by Romania’s territory, a country located in the central-eastern part of Europe with an area of 238,391 km² (Figure 1). Romania has a vast diversity of landforms, each with representative forestry variation, including the Transylvanian Depression located in the center of the Carpathian Mountains arc. This territory offers proper conditions for the extension of various types of deciduous forests (i.e., with species such as Carpinus betulus, Fagus sylvatica, Tillia tomentosa, Ulmus minor, Quercus petraea) over large areas. The Romanian Carpathian Mountains (Carpathians) occupy 57% of the country’s territory, where extensive forests are spread over large areas.

According to Romanian legislation, a forest is considered an area of at least 0.25 hectares of land occupied by forest vegetation. The mature specimens reach 5 m in height under normal vegetation conditions and have a coverage index (consistency) of more than 10% (0.1). These territories are added to the areas covered with junipers (Juniperus) in the high mountain area of over 1800–2000 m in altitude and forest protection curtains with more than 0.5 hectares and a width of more than 20 m. Forest protection curtains are accepted crucial interventions in forest protection policy. They are projected to have a significant protective effect on Romanian forest cover on the brink of climate change [23].

At the territorial level of Romania, 29.9% of the surface is covered by different forests, covering 7.13 million hectares [24]. Romania is one of the countries with the highest percentage of occupied forest areas, with the latest estimates having a significant growth rate (i.e., 19.3 million m³/year for conifers, 19 million m³/year for beech, 8.1 m³/year for quercinea, 8.6 million m³/year for hardwoods, and 3.4 million m³/year for softwoods), to which are added old forests and virgin forests in different stages of conservation. Large forest areas are predominantly in the mountainous and hilly areas and areas with lower altitudes. There is a higher density of human settlements, a crucial aspect considering the present study’s objectives. We have included further details about Romania’s forest structure in

Table A1. Distribution of forests in Romania by consistency categories at the level of significant relief units.

Major Relief Unit	Consistency Categories
	Arboretum with Full Consistency (C = 100)		Arboretum with Almost Full Consistency (C = 70–90)		Arboretum with Sparse, Lighted, or Braced Consistency (C = 40–60)			Arboretum with Degraded Consistency (C = 10–30)
	Km2	%	Km2	%	Km2	%	Km2	%
Romanian Plain	45,459.89	92.9	881.28	1.8	1874.65	3.8	703.66	1.4
Getic Plateau	8358.95	60.5	1647.01	11.9	3580.86	25.9	234.7	1.7
Western Hills south of Mures	2458.28	54.3	394.81	8.7	1559.46	34.4	115.32	2.5
Banat Mountains	1429.78	20.4	612.92	8.7	4568.48	65.2	399.24	5.7
Transylvanian Depression	18,164.37	71.8	1970.75	7.8	4593.99	18.2	568.38	2.2
Apuseni Mountains	2319.16	21.7	1102.65	10.3	5275.04	49.4	1973.78	18.5
Moldavian Plateau	19,059.32	83.1	802.81	3.5	2766.84	12.1	300.88	1.3
Subcarpathian	8206.62	49.5	1941.2	11.7	4867.73	29.3	1574.77	9.5
Dobrogea Plateau	9136.14	60.5	160.5	11.9	765.5	25.9	90.79	1.7
Danube Delta	4289.01		72.52		87.99		10.18
Mehedinti Plateau	333.11	41.6	112.58	14.1	350.02	43.7	5.46	0.7
Western Plain	15,337.49	97.8	236.43	1.5	504.67	3.1	99.38	0.6
Eastern Carpathians	8753.47	30.9	2738.64	9.7	12,664.1	44.8	4128.64	14.6
Curvature Carpathians	954.29	15.6	429.88	7.0	2605.22	42.5	2141.88	34.9
Southern Carpathians	3852.61	27.2	1092.47	7.7	7488.87	52.8	1741.48	12.3
Western hills north of Mures	5258.13	63.3	649.76	7.8	2090.52	25.2	305.84	3.7
Total (km2)	15,3370.62		14,846.21		55,643.94		14,394.38

(see Appendix A).

Changes over time in the areas covered by forest are under the direct impact of natural factors such as the influence of climate change on the consistency and composition of forests, the migration of forest species beyond known ecological limits, the negative influence of floods with short return periods, the degradation of physical and chemical properties of soils due to soil erosion, and vegetation fires with natural causes. Anthropic changes are also present due to deforestation caused by logging, legal and illegal, whose rate increased after 2000 [25,26] and contamination of soils and groundwater supplies [27]. However, some territories showed a forest gain due to the afforestation of large areas of abandoned pastures [26]. Furthermore, the spatio-temporal evolution of forest cover in Romania is tightly correlated with the forest management regimes affected by socio-political fluctuations from the early 19th century [28].

2.2. WIPI/WSCI Model and the Current Updates

This study methodologically relies on the Wildfire Ignition Probability/Wildfire Spreading Capacity Index (WIPI/WSCI). Initially, the method defines criteria that have proven relation with either the wildfire occurrence or behavior. The number of criteria varies according to the available data and the specifics of the study area. In this study, we shortlisted 16 criteria about the anthropogenic (S-social), hydrometeorological (E-environmental), geophysical (P-physical), and fuel (F) properties of the study area (Romania) following our earlier GIS-based method [19]. Through a literature review and evaluating the available open access geospatial data, we considered the following criteria: population density (S1), distance to settlements (S2), distance to transportation network (S3), distance to main roads (S4), agriculture distance (S5), solar radiation (E1), precipitation (E2), maximum temperature (E3), wind speed (E4), slope (P1), aspect (P2), altitude (P3), distance to water sources (P4), fuel type (F1), Tree Cover Density (TCD) (F2), and Normalized Difference Vegetation Index (NDVI) (F3).

The method starts by defining a regular grid of points. According to the inventory phase results, there are 70,410 reference point locations within Romania’s vegetated surfaces. The distance between points is 1 km, and each point represents a vegetated surface of 1 km². Each reference point location within the vegetated surface is loaded with unique absolute values through a multi-criteria inventory procedure [19]. Each criterion’s relative weighted factor was initially assigned via AHP pairwise comparison method. The sensitivity of the model has been assessed via ROC/AUC method in another study focusing on the case of Montenegro [29].

Figure 2 presents the methodical workflow of this study. It includes the updates that we push forward as improvements of WIPI/WSCI, applied in Romania’s case. At this stage, the workflow consists of seven sequential stages. Besides the inventory procedure, the first stage includes defining the vegetated surfaces within the study area and the reference points that spatially represent the vegetation surfaces. The reference points serve as data collecting pivots loaded with all 16 criteria’ unique values, as shown in Figure 2. Among a regular points grid with a spatial distance of 1 km, we selected by location the points that overlapped with the vegetated surfaces derived from CORINE Land Cover (CLC) data.

Following the inventory phase, the unequal range of inventory values necessitates a normalizing procedure before indexing calculations. This stage equalizes the range of inventory values of each criterion into a gradient between 0 and 1. The max/min normalizing procedure is selected as it is accepted as the most suitable and straightforward method for well-known sets of records [30].

The third stage consists of subgrouping the criteria into two sets according to their relationship with either wildfire ignition or spreading (see Figure 2, third stage). This division is based on a literature review shown in our earlier work [19]. Moreover, a relevancy indicator is given to each criterion according to their direct or indirect relationship with wildfire regimes. This is explained in detail in

Table 1. The weighted impact factors (WIPI and WSCI) of each criterion via AHP and ROC/AUC methods.

				WIPI				WSCI
		ROC/AUC	AUC/Norm	+ −	AHP	ROC	AUC/Norm	+ −	AHP	ROC	AUC/Norm
E1	Solar radiation	0.823	0.110	+	0.032	0.823	0.110	+	0.011	0.823	0.103
E2	Precipitation	0.201	0.027	−	0.097	0.799	0.107	−	0.048	0.799	0.100
E3	Max. Temp.	0.811	0.109	+	0.032	0.811	0.108	+	0.022	0.811	0.101
E4	Wind speed	0.575	0.077					+	0.155	0.575	0.072
F1	Fuel type	0.593	0.080	+	0.056	0.593	0.079	+	0.033	0.593	0.074
F2	TCD	0.388	0.052					+	0.299	0.388	0.048
F3	NDVI	0.405	0.054	−	0.125	0.405	0.054	−	0.170	0.595	0.074
P1	Slope	0.278	0.037					+	0.033	0.278	0.035
P2	Aspect	0.511	0.069	+	0.049	0.511	0.068	+	0.013	0.511	0.064
P3	Altitude	0.192	0.026	−	0.016	0.808	0.108	−	0.006	0.808	0.101
P4	Dist. To water	0.411	0.055					+	0.070	0.411	0.051
S1	Pop. density	0.501	0.067	+	0.026	0.501	0.067
S2	Dist. Settlements	0.381	0.051	−	0.076	0.619	0.083	+	0.017	0.381	0.048
S3	Dist. roads	0.551	0.074	−	0.140	0.449	0.060	+	0.006	0.551	0.069
S4	Dist. main roads	0.477	0.064	−	0.045	0.523	0.070	+	0.047	0.477	0.060
S5	Dist. agriculture	0.354	0.048	−	0.305	0.646	0.086

. The first three methodical stages are borrowed from our previous studies.

In the fourth stage, we propose ROC/AUC analysis (via SPSS software) as a weighting method among criteria, besides the AHP pairwise comparison method. This relies on the specific characteristics of the study area and historical data on fire regimes. Indexing values are calculated as the sum of the products between inventory value and each criterion’s impact factor, as shown in Equations (1) and (2). Inventory values of criteria are freely accessible as explained in Data Availability Statement at the end of this article. Table 1 delivers the impact factors of criteria as calculated via AHP and ROC analysis.

$$WIPI={\sum }_{i=1}^n{\alpha }_i{C}_i$$

(1)

where WIPI is the normalized wildfire ignition probability index, $C_i$ is the inventory value of criterion i, and ${\alpha }_i$ is the weighted impact coefficient of criterion i.

$$WSCI={\sum }_{j=1}^m{\beta }_j{C}_j$$

(2)

where WSCI is the normalized wildfire spreading capacity index, $C_j$ is the inventory value of criterion j, and ${\beta }_j$ is the weighted impact coefficient of criterion j.

We compare the earlier model results (WIPI/WSCI) and the updated one (WIPI_ROC/WSCI_ROC) as applied in Romania’s case. During the sixth stage, the ROC/AUC method assesses both models’ accuracy, leading to a comparative discussion. At the final stage, the WSCI_ROC model results are used in vulnerability assessment of protected areas and exposure analysis of urbanized zones.

2.3. Data Acquisition

This study depends on a variety of free access to remotely sensed geospatial data. These data are acquired from various sources. We have included detailed information in

Table A2. The open-source raw data utilized in this study.

	Data	Type	MMU	Utility within the Method	Source
1	CORINE land cover	Vector	25 ha	Vegetated surfaces; Settlements (S2); Fuel type (F1); Distance to agriculture (S5)	Copernicus/EEA
2	DEM	Raster	25 m	Slope (P1); Aspect (P2); Altitude (P3)	Copernicus/EEA
3	NDVI	Raster	250 m	NDVI (F3)	Earth data/NASA
4	TCD	Raster		TCD (F2)	Copernicus/EEA
5	Weather data	Raster	30 s	Solar radiation (E1); Precipitation (E2); Maximum Temperature (E3); Wind speed (E4)	Worldclim 2.0
6	Burned area fraction	Raster		ROC/AUC analysis	Climate Change Service
7	Burned areas	Vector		ROC/AUC analysis	EFFIS/JRC
8	OSM	Vector		Distance to water (P4); Distance to any road (S3); Distance to the main road (S4); Exposure analysis	OSM Geofabrik
9	Population Density	Raster		Population density (S1); Exposure analysis	National Institute of Statistics (RO)
10	Administrative/NUTS	Vector		Exposure analysis	EUROSTAT
11	Protected areas	Vector		Vulnerability analysis	Eionet/EEA

(Appendix B), which presents the complete list of the data name, data type, Minimum Mapping Unit (MMU), source, and utility within the method’s workflow. CLC is a pan-European data provided by the European Environment Agency (EEA), which delivers a hierarchical classification of 44 land cover types [31]. The original classification method simultaneously relies on both the Sentinel-2 satellite imagery (i.e., the 1st dedicated European satellite for land monitoring) and Landsat-8 images for gap-filling.

In this study, we rely on the CLC data of 2018 to gather geospatial information about vegetation surfaces, settlements (S2), fuel type (F1), and agricultural areas (S5). The precise data extraction method regarding fuel type and the relative weighting is explained in our previous work [20]. Land cover classes of vegetated surfaces were ranked according to their relative implication with either wildfire ignition or spread. For example, broad-leaved forests (CLC-311) lead the list of correlation with wildfire spreading capacity. EEA supplies other data such as Digital Elevation Model (DEM) and Tree Cover Density (TCD). DEM is delivering information about slope (P1), aspect (P2), and altitude (P3) in raster format of 25 m in resolution. The population density information is produced based on population records at the smallest local administrative unit level, as shown in Figure A1 (Appendix B) based on the data provided by the National Institute of Statistics of Romania. All maps in this study are projected according to ETRS89/ETRS-LAEA (ESPG: 3035) projection reference system in QGIS.

Normalized Difference Vegetation Index (NDVI) data is extracted from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) Version 6. These data are produced in a resolution of 250 m by choosing the most reliable pixel value among daily values within 16 days. The low percentage of cloud coverage, low view angle, and the highest NDVI value are among the applied selection criteria [32,33]. In this study, we used the period between 28 July and 12 August within the fire season of 2018.

Raster data about solar radiation (E1), precipitation (E2), maximum temperature (E3), and wind speed (E4) are derived from WorldClim 2.0 database (http://www.worldclim.com/version2, accessed on 3 June 2020). It consists of raster images of 30 s (1 km²) resolution that provide monthly average values recorded between 1970 and 2000. These products rely on an interpolation method using historical data collected by 60,000 weather stations around the globe. Furthermore, remotely sensed data like maximum and minimum land surface temperature and cloud cover, as obtained via the MODIS satellite platform, were used as satellite-derived covariates for accuracy improvement [33]. Fick and Hijmans provide further details about the production process of remotely sensed data via Google Earth Engine in their previous work [34]. In this study, we use August’s average records as the weather conditions for wildfire spread are highest. The remaining criteria, like distance to water (P4) and transportation network (S2 & S3), stand on Open Street map (OSM) data enabled for free via the Geofabrik portal (https://www.geofabrik.de, accessed on 27 December 2020).

Another crucial data set used in this study is the Burned Area (BA) products acquired from Copernicus Climate Change Service (2019). They provide information about the total BA at the pixel level (250 m). The results are prepared via reflectance change analysis of medium resolution sensors like Terra MODIS, Sentinel-3 OLCI combined with the thermal data by MODIS. These products are vital raw data for research that focus on themes like climate change, land use and land cover dynamics, wildfire risk assessment, among others [35]. MODIS products are the most widely used global dataset by the scientific community [36,37,38].

According to the Burned Area Fraction (BAF) maps reported by Copernicus Climate Change Service during the fire season between 2015 and 2019 (Figure 3), 1980 points overlap with surfaces marking a cumulative BAF above 100%. In this study, this threshold is assigned considering that a cumulative BAF value above 100% (the sum of 5 years records) is considered a high probability that a particular location experienced a fire during the 5 years. Apart from BAF, we utilized the recorded large fires as derived from the EFFIS database. The wildfire occurrence and size are recorded based on the MODIS burned area product (MCD64A1) following the Globfire method [39]. MCD64A1 products use Terra and Aqua satellite imagery, resulting in dependable identification of large fires in different locations [40,41,42]. Nevertheless, Moreno-Ruiz et al. [43] conclude that there is still space for further improvement in the MODIS sensor-based burned area detection algorithms. Further details about the data used in this study are included in Table A2 (Appendix B).

3. Results

3.1. Multi-Criteria Inventory of Wildfire-Related Factors in Romania

First, the method delivers inventory results on an individual level per each criterion. Figure 4 presents the relative wildfire proneness map of vegetated surfaces (Figure 4q) in Romania based on each criterion. The color palette is set as the gradient of red-yellow-green, where red shows the highest risk areas while green the least risk. The range of inventory values are reclassified into 10 classes following Jenks natural break normalizing procedure. The gradient direction is assigned according to the relative indicator, as explained in Table 1. For example, criteria like solar radiation (E1) and precipitation (E2) are shown under reversed color gradient. In other words, the highest solar radiation values indicate the highest risk. In contrast, the highest precipitation records correlate with the lowest risk.

According to Figure 4, the Carpathians that cross the Romanian territory in the central region from north to south-west direction stand as a determinant of the spatial distribution of the relative risk for the majority of the criteria. First, the hydrometeorological criteria visually correlate with the topography of the study area. Figure 4a shows that solar radiation (E1) is higher at lower altitudes, especially in Carpathians’ southeast, and lower at high altitudes. Similarly, the recorded maximum temperatures (E3) are on the same line as the altitude (P3) values presented in Figure 4j.

Wind speed (E4) is the only environmental criterion that is not tightly correlated with altitude (Figure 4d). The slopes of Carpathians face south and southeast and remain an exposed area to winds flowing from the Black Sea and the Mediterranean. Similarly, the remaining plain territories in Romania’s south-eastern region facing the Black Sea are exposed to considerable average wind speeds compared to the north-western plains (Figure 4d). Among geophysical criteria, slope (P1) is the only criterion that correlates with the altitude values (Figure 4h). Most of the sloped surfaces are found along with the Carpathian’s layout; whereas the aspect (P2) values are more uniformly dispersed in the territory (Figure 4i) and distance to water surfaces (P4) follows the spatial distribution of water elements in the landscape (Figure 4k), independently to the altitude values.

The visual results discussed earlier are on the same line as the correlogram presented in Figure 5. We performed a correlation analysis among all 16 criteria and BAF values for 70,410 point locations to determine which criteria correlate more with wildfire records of 2015–2019. We used the “corrplot” package in RStudio [44], to produce the correlogram. The default calculation is based on the Pearson correlation method. The following correlogram delivers both visual plotting and numerical evidence of the correlation coefficient of each pair of criteria. The larger the circle and the darker the blue color, the higher the positive correlation among criteria. Similarly, the larger the circle and the darker the red color, the higher the inverse correlation among the two variables.

On the other hand, the smaller the circle, the less the correlation between the two variables. For example, there exists a significant positive correlation between precipitation (E2) and altitude (P3) (0.84). Similarly, there exist significant negative correlations between solar radiation/maximum temperature (E1/E3) and altitude (P3) (−0.81/−0.96). Cumulative BAF values (monthly, 2015–2019) positively correlate highest to solar radiation (E1) and maximum temperature (E3), respectively, 0.43 and 0.47; while they inversely correlate highest to precipitation (E2) and altitude (P3), respectively, −0.43 and −0.44.

3.2. Calibrated Weighting of Criteria

We adopted the weighted factors as assigned via an AHP pairwise comparison method from our previous study [20], as listed in Table 1. AHP was applied in two different levels and was based on experts’ evaluation. First, we compared four major categories of criteria (anthropogenic, hydrometeorological, geophysical, and fuel). Then we assigned a second level coefficient to each criterion by comparing within each category. The current coefficient is shown in Table 1 (AHP) and is the product of two coefficients. According to our previous results, distance to agriculture (S5) is the criterion that has the highest impact on wildfire ignition probability, almost one-third of all wildfire ignition coefficient distribution. On the other hand, the criterion that affects the most wildfire spreading capacity of vegetated surfaces is tree cover density (F2). Our previous study suggested that sensitivity analysis via ROC/AUC method should be tested in assigning weights of the impact factors.

According to each criterion’s sensitivity analysis concerning the historical fire regimes, the revised weighted impact factors are assigned. Figure 6 presents the ROC analysis values per each criterion under four groups; hydro meteorological-environmental (E), fuel (F), geophysical (P), and anthropogenic-social (S). The sensitivity analysis is performed in SPSS software via the ROC analysis tool. One-thousand-nine-hundred-and-fifty-six points out of 70,410 have a 5-year cumulative burned area fraction (2015–2019) value above 100%. These points are considered positive samples in the ROC analysis procedure in finding each criterion’s sensitivity. On the other hand, there are 41,683 points (about 60%) that have a BAF value of 0, indicating the areas that have been safest during the fire season of 2015–2019.

The results presented in Figure 6 reveal the hypothetical models’ sensitivity that has a single determinant, being each criterion. It is a way to find the sensitivity of the model to each criterion based on the burned area fraction in Romania (see Figure 3). According to Figure 6, the highest AUC value belongs to E1 (solar radiation) and E3 (maximum temperature), respectively, 0.823 and 0.811. While the lowest AUC values are recorded for E2 (precipitation), P1 (slope), and P3 (altitude), respectively, 0.201, 0.278, and 0.192. In principle, the lower the AUC value, the higher the inverse correlation between the criterion and the wildfire recorded burned area fraction. The criteria that score an AUC value less than 0.5 have an inverse correlation with the wildfire records. These results converge with the correlation matrix analysis presented in Figure 5.

The first two columns of Table 1 present the absolute and normalized AUC values of all criteria. The following columns deliver each criterion’s relative impact factors as calculated via AHP pairwise comparison and ROC/AUC analysis. Besides, each criterion is assigned an indicator for either direct (+) or inverse (−) relation with the wildfire risk. This indicator is assigned based on assumptions inferred from the literature review on the relationship between wildfire regimes and driving factors [20]. For example, the higher the solar radiation, maximum temperature, fuel type, and aspect values, the higher wildfire ignition and spread probability. On the other side, the lower the precipitation, NDVI, and altitude values, the higher the wildfire risk. Simultaneously, anthropogenic criteria like distance to settlements and transportation networks are unevenly related to wildfire ignition and spreading phases.

3.3. Comparing between WIPI/WSCI and WIPI_ROC/WSCI_ROC Results

We calculated the WIPI and WSCI index values of each reference point according to Equations (1) and (2) using the weighted values via the AHP method. The results of the ROC/AUC analysis show the relatively low sensitivity of the model. An accumulative BAF record point higher than 100% (see Figure 3) was accepted as a positive case in the model. According to Figure 7, the AUC of WIPI is 0.705 and has an overall model quality of 0.69. It is significantly higher than the WSCI model sensitivity marking an AUC value of 0.587 and an overall model quality of 0.57. It can be inferred that the model, which relies on the weighted values calculated via AHP, is more accurate for predicting wildfire occurrence events rather than the wildfire spreading process.

Later, we recalculated each reference point’s index values as the sum of the products between normalized inventory value and the weighted factor via the ROC/AUC method (Figure 2). The revised weights, as presented in Table 1, led to improved model sensitivity. Figure 8 presents a comparative ROC/AUC analysis between the former WIPI/WSCI and the revised WIPI_ROC/WSCI_ROC models.

The updated models’ curves are shown in green, while the previous versions are red. According to Figure 8a, the AUC value of WIPI_ROC increased to 0.778, marking a sensitivity increase of 0.073. The overall model quality of the WIPI model was improved by 12% (from 0.69 to 0.77). The improvement is more visible in the case of the WSCI model (Figure 8b). The revised model (WSCI_ROC) recorded an AUC value of 0.802, 37% higher than the earlier version. A similar escalation is recognized in the overall model quality. We rely on the WSCI_ROC results during the final stage of vulnerability and exposure analysis.

Beyond statistical analysis about the model accuracy, the results of the WSCI_ROC model deliver essential findings of the spatial distribution of the wildfire, spreading the capacity risk of vegetated surfaces in Romania. According to Figure 9, the highest wildfire spreading risk is concentrated in its eastern and southern regions. A secondary area under wildfire spreading capacity risk is along the western borders. Simultaneously, the central and the northern regions appear to be safer from wildfire, spreading the risk. These regions coincide with the surfaces that are high in elevation and located away from human activities. Among them, the sub-regions that record the lowest WSCI_ROC values are surfaces that are oriented towards the northern direction and gaining a minimum of solar radiation. The gradient of green color indicates these surfaces.

Furthermore, Figure 9 includes the administrative boundaries of the third level (Classification of Territorial Units for Statistics, NUTS-L3) of Romania. Referring to the Eurostat data, Romania has 42 local administrative units at the third level. According to the results presented in Figure 9, 41 units consist of at least 1 km² of vegetated surface that has been indexed here. The only unit that has no wildland vegetated surface is the capital city of Romania. As highlighted in Figure 9, Bucharest is the smallest in surface area compared to the other units. However, it does not mean that it is safer. On the contrary, when jointly considered with the metropolitan area of Ilfov, which envelopes the urban area of Bucharest, the wildfire ignition and spreading risk in the wildland-urban interface (WUI) is critical and worth investigating. We further discuss this issue in the following sections while assessing human-made structures’ exposure to wildfire spreading risk.

The overlapping of wildfire spreading capacity risk and the local administrative units highlight the municipalities that need to enhance their wildfire prevention measures. Figure A2 (Appendix C) presents the box plots that show the WIPI_ROC and WSCI_ROC values distribution by local administrative units. According to the results plotted in Figure A2a, the administrative units with the highest WIPI_ROC index values are located in the country’s south-eastern and southern regions.

4. Discussion on the Vulnerability of Protected Areas and Exposure of Settlements

Wildfires are native events on earth estimated to have happened during the last 350 million years [45]. They are accepted to significantly contribute to vegetation recovery with their natural schedule, the biogeochemical cycles of carbon and nitrogen, and the atmosphere’s chemical properties [46,47,48]. However, they are also reported to have considerable consequences on the territory’s ecological and social systems. Protected areas are among the land surfaces where the ecological and socio-cultural interests converge with each other. Uncontrolled fires in the vegetated protected areas may cause unrecoverable consequences on the native vegetation structure. On the other side, the WUI zone is boldly highlighted in scientific reports. Wildfire events threaten human activities. We expand our findings by discussing further concerning protected areas’ vulnerability and exposure to human-made structures.

4.1. The Vulnerability of the Romanian Protected Areas to Wildfires

Globally speaking, the protected areas are under consistent threat caused by climate change processes, land-use alterations, provisioning of raw materials, socio-cultural activities, and flourishing of invasive species [49,50]. Jones et al. [51] conclude that only two-thirds of the protected areas are safe from globally intensive human activities. While Schulze et al. [52] list fire and fire suppression activities as the third out of 36 threats that protected areas usually face according to the list of level two threats included in the IUCN-CMP Threats Classification Scheme. Forest fires can lead to an invasive plant expansion in disturbed sites [53]. The native species in Romania’s protected areas have been at risk of several natural and human-induced hazards [54]. Thus, assessing the vulnerability of the protected areas to wildfires is of great concern in Romania.

This assessment relies on the European inventory of nationally designated protected areas, as acquired from the EEA open-source datasets. According to these data, Romania has 946 protected sites, covering a total area of 13,985 km². Figure 10 presents the protected areas overlapping the WIPI_ROC results. Most of these areas are found in the alpine lands along with the Carpathians. About 13% of the 70,410 reference points within the vegetated surfaces we analyzed are located within the protected areas. Consequently, we may infer that 65% of Romania’s protected surfaces are vegetated and potentially vulnerable to wildfire risk.

According to the box plot in Figure 10a, the protected areas have lower WIPI_ROC and WSCI_ROC values than unprotected surfaces. Furthermore, referring to Figure 10b, most vegetated surfaces within the protected patches are greenish, implying a relatively low wildfire spreading capacity (WSCI_ROC). These values’ main reasons are the high elevation and remote location of protected surfaces to human activities (settlement and transportation network). However, some cases consist of a gradient of WSCI_ROC values within the same protected surface (see the enlarged protected patch in Figure 10).

The relation between vegetated landscape patches and wildfire spreading risk has been a questionable literature topic [55]. O’Donnell et al. [56] report that the fragmentation among vegetated landscapes in unmanaged Australian semi-arid shrublands and woodlands directly impacted the reduction of wildfire intervals between 1940 and 2006. Thus, the connectivity among vegetated surfaces within the same protected area may boost the wildfire, spreading greenish areas’ risk. Future studies must focus on protected patches to assess the vulnerability to wildfire spreading risk at a finer spatial scale.

4.2. Wildfire Exposure of Populated Areas within the Metropolitan Area of Bucharest

The urban fringes are critical hybrid areas where the human-made structures are exposed to different environmental hazards [57]. Studies from developing countries report that uncontrolled urban expansion increases inhabited surfaces’ exposure to natural hazards like floods, landslides, fire, and sinkholes, among others [58]. Sestras et al. [58] report landslide assessment at a local scale as an inherent threat in Romania’s newly developed suburban zones. The wildland–urban interface represents an area of contradiction where both the settling interest and wildfire risk are significantly high [59,60,61].

We performed an exposure assessment of built structures to the wildfire spreading capacity of vegetated surfaces within the Romanian capital city’s metropolitan area (Figure 9). The wildfire exposure analysis relies on the juxtaposition between the WSCI_ROC results and existing building stock, as shown in Figure 11. We bring a demonstrative example from Ilfov metropolitan area, which includes the capital city, Bucharest. The hazard map of wildfire spreading capacity relies on the results reported in this article (see Figure 9). The WSCI_ROC point’s layer is utilized for preparing the hazard heatmap (kernel density estimation) with a selection radius of 1 km.

OSM data provide the building stock within the focal study area. Nevertheless, OSM data accuracy is still debatable due to the lack of professional backgrounds [62]. A substantial number of building features do not include building-use information. However, we bring this discussion forward as an exposure analysis method that can deliver critical information about the human-made structures under wildfire risk within WUI areas at a metropolitan scale if further improved by introducing validated building stock data.

According to Figure 11b, 9596 structures overlap with the WSCI_ROC heatmap, exposing the minimum wildfire spreading capacity. Furthermore, 1734 out of the total exposed structures are building types of significant socio-economic. This building stock’s exposure to wildfire spreading risk may have critical consequences on socio-economical processes and their users’ life security. Figure 11c presents the box plot of WSCI_ROC heatmap values distribution per building type. We included just the most critical building types like hospitals, industrial, residential, religious, leisure, and educational use. While other buildings like warehouses, greenhouses, and abandoned structures were ignored at this stage.

According to the box plot, just one hospital is located within the critical wildfire spreading capacity heatmap. Nevertheless, it has a WSCI_ROC heatmap value of 2.37, which indicates a significant exposure of its users to wildfire spreading risk. According to the outlier values (dots) shown in Figure 11c, there are four industrial, one house, and one school building highly exposed to wildfire, which can spread the risk.

5. Conclusions

This study presented the indexing of vegetated surfaces in Romania by Wildfire Ignition Probability and Spreading Capacity Index (WIPI/WSCI). The method offered here relies on open-source and available remotely sensed data, which supplies the analytical process with geospatial information about the anthropogenic, hydrometeorological, geophysical, and fuel properties of the study area. We identified 16 criteria that significantly impact either the wildfire ignition or spreading phase of the wildfire event in Romania. The impact of each criterion on wildfire is weighted via ROC/AUC analysis. During the analysis, the positive cases rely on the burned area fraction records between 2015 and 2019 (5 years).

According to our results, solar radiation (E1), precipitation (E2), and maximum temperature (E3) have the highest correlation to the cumulative burned area fraction reported for the 2015–2019 fire season. The vegetated surfaces in Romania’s eastern and southern regions face the highest wildfire spreading capacity index values. Considering that these regions make home to urbanized lands of high population density, the high WSCI records indicate an elevated risk. We performed the wildfire spreading risk exposure analysis of the building stock within the Bucharest metropolitan area, Ilfov. These results imply that critical structures like hospitals and residential and educational units are found at locations of significant risk.

On the other side, Romania’s central areas scattered along the Carpathians have the lowest index values. This is generally driven by high altitude values, which are directly correlated with other climatic criteria such as solar radiation, precipitation, and maximum temperature. These regions include most protected areas. According to our results, some protected vegetated surfaces in Romania hold a gradient of wildfire spreading risk within the same protected area geometry. In such cases, the whole area within the protection borders must be considered under risk, as the connectivity among vegetated surfaces in wildfire risk analysis is considered a weakness.

The method we presented in this study is reproducible in other wildfire-prone geographies. It is also flexible enough to integrate the most up-to-date and most reliable remotely sensed geospatial data available to the community in the future. The results presented in this study can help the institutions at the national and local levels responsible for wildfire risk reduction in Romania.