Modeling the spatial evolution wildfires using random spread process

The study of wildfire spread and the growth of the area burned is an important task in ecological studies and in other contexts. In this work we present a model for fire spread and show the results obtained from simulations of burned areas. The model is based on probabilities of fire at different locations. Such probabilities are obtained from the intensity function of a spatial point process model fitted to the observed pattern of fires in the Valencian Community for the years 1993–2015. The models, applied to different wildfires in Spain, including the different temporal states, combines the features of a network model with those of a quasi-physical model of the interaction between burning and nonburning cells, which strongly depends on covariates. The results of the simulated wildfire burned areas resemble the burned areas observed in real cases, suggesting that the model proposed, based on a Markov process called Random Spread Process, works adequately. The model can be extended to simulate other random spread processes such as epidemics.


F I G U R E 1
A wildfire in the Valencian Province, august 2012 (left) and its effect on landscape (right) Modeling of fire spread in two dimensions can be done using vector and raster-based approaches.For vector based simulation of fire spread, it is necessary to simulate fire as a continuously expanding polygon (Anderson et al., 1982).This is the base of the FARSITE model (Fire Area Simulator, Finney & Andrews, 1994), of wide use in the area of fire ecology.However, the vector approach is highly dependent on the spatial scale used and requires the evaluation of a complex convex hull, which is computationally expensive.
An alternative approach is to use raster schemes, in which the domain of interest  is partitioned in square or hexagonal cells that form a lattice covering the domain.Fire spreads to neighboring cells using contact cell methods or using heath-accumulation methods (Berjak & Hearne, 2014;Clark et al., 1996;Frandsen & Andrews, 1979;Green et al., 1990;Hargrove et al., 2000;Kourtz & O'Regan, 1971;Yassemi et al., 2008).
Two main approaches to simulate fire spread in a rasterized domain are cell contact and heath accumulation (Green et al., 1983).The contact cell method has high computation efficiency.It interprets fire spread as a series of ignitions that span the area of an individual cell.(Frandsen & Andrews, 1979;Green et al., 1990;Kourtz & O'Regan, 1971) used this approach extensively, but their events are generated considering only the influence of the fastest spreading neighbor, neglecting that fire spread to a new unburned cell may be based upon the effect of multiple neigh-boring cells or prior heating (Green et al., 1983).
The heat accumulation approach considers the rate of spread of fire into a cell or pixel at a given time t as the sum of the contribution of ignited cells in its neighborhood during previous times.Despite its appealing to simulate fire spread, it is demanding computationally and requires the evaluation of some model parameters empirically (Yassemi et al., 2008).
In this article, we consider a random spread model, based in a raster contact approach.Fire ignitions are assumed to ignite at random within the interest domain  according to an underlying spatial point process  with intensity function (x) and density The ignition point x lies inside a pixel or cell v whose state changes from "unburned" to "burned".Note that this is equivalent to choose a pixel at random in  with probability proportional to f (x).Conditional on an ignition occurring at pixel v, neighboring pixels burn with probability proportional to the intensity function.We show how the model works using data from two wildfires occurring in 2016 and compare model performance using receiver's operation curve (ROC) and precision-recall (PRC) curve analysis.Simulations to test the model performance and to simulate the Artana and Catcaixent wildfires spread were done using a program code written in R-language.
The rest of the manuscript is organized as follows: In Section 2 we describe the model and used and provide all the details needed to clarify the methods used.In Section 3 we describe the data sets of this study.Section 4 presents the results and discussion with their implications in terms of computing efficiency in model fitting.Finally, Section 5 presents the conclusions, followed by references.

Model for fire spread
To model the propagation of wildfires in the study area, we will use a raster approach.Let  ∈ R 2 be a geographic area of interest, partitioned in N sub-areas or pixels v it ; i = 1, … , n, observed at discrete but arbitrary times t 0 , t 1 , … , T. We will assume that at any time t each pixel belongs to one of the following mutually exclusive categories: Burned (B t ), fire front (F t ), neighboring the fire front (V t ) and unburned (U t ) (Figure 2).Note that  t = B t ∪ F t ∪ V t ∪ U t and that for t + 1 some pixels in the complement of B t ∪ F t may join F t+1 .We assume that for each pixel there exist a probability probability of fire ignition p i which depends on local conditions through a set of covariates Z 1 , … Z p .This assumption allows to capture the spatial variability of fire probability often observed in real situations.For example, pixels covered by some vegetation types have higher fire probability compared to others (Boychuk et al., 2007(Boychuk et al., , 2009;;Díaz-Avalos et al., 2001).
Our interest is to model the transitions of pixels from V t to F t ∪ B t+1 .We will assume a null probability that pixels in V t at any time might caught fire from sources other than sparks and heat from contact with pixels in B t .
In an analogous way to Markov Random Fields, we assume that each pixel is surrounded by a set of neighboring pixels N i .This set of neighbors is defined in our case as the set of nearest neighbors of order k, k = 1, 2. Given that a fire has ignited at site s i ∈ W, it spreads to its neighbors s j ∈ V i according to their fire probabilities p(s j ).We assume that p(v) = 0, s ∈ B and that once a fire has ignited at a pixel v i at time t, v i becomes a member of B for t + k, k = 1, … , K. The p(v j ) may be obtained from the intensity function (x i ) as is described in the next subsection.

Estimation of the p(v i )
To obtain the fire probabilities at each pixel in , we first estimate the intensity function given the observed pattern of wildfires in  using methods for point spatial process modeling (Baddeley et al., 2016).Because we are considering only large sized fires, we fitted an area interaction point process model to the observed fire pattern (Baddeley & van Lieshout 1995).Area interaction point process models are useful to model spatial point patterns with moderate clustering to moderate repulsion (Aragó et al., 2016;Baddeley & van Lieshout, 1995).For any point pattern x = {x 1 , … , x n } ∈ , this class of models has a density function of the general form where A(x) = [∪ n i=1 B(x i ; r)] ∩  is the set formed by the union of discs of radius r centered at the events x i of the point process; ,  and r are model parameters and  is the normalizing constant.The interaction parameter  can be any positive number.If  = 1, then the model reduces to a homogeneous Poisson process, if  < 1, then the process is regular, while if  > 1 the process is clustered.Two points interact if the distance between them is less than r.For any finite point process, the Papangelou conditional intensity function is defined as (Coeurjolly et al., 2015) his is, for a given point pattern x, the conditional intensity function for a new event at u ∈  is the conditional density function of u given the point pattern x.Given a spatial point pattern x = {x 1 , … .xn with conditional intensity  ( u, |x ) , parameter estimation and statistical inference for the intensity function is done using the pseudo-likelihood function, and therefore, the log pseudo-likelihood is If the integral in the above equation is approximated by some quadrature rule, the approximated pseudo-likelihood becomes where u j are quadrature points and w j are quadrature weights.Details on the estimation method are beyond the scope of this work, but the interested reader may find a detailed technical description in (Diggle, 2021;Diggle et al., 1994).We included Land Use, elevation, Slope, Aspect and Combustion Index as covariates in our model.Land use is a categorical variable with the coding shown in Table 1.The covariate values for each pixel and each wildfire were allocated to the entries of a matrix Z(x).The covariate values for each wildfire in the data base were allocated using GIS software.The logarithm of the conditional intensity function for the Area-Interaction point process model takes the general form: (5) We point out that Z(x) is a matrix with the values of covariates at location x in its columns,  is a vector with regression coefficients (Renner et al., 2015;van Lieshout, 2009).We consider as wildfires those forest fires with final size above 10 ha, which correspond to the 20th percentile of the empirical distribution of the wildfires recorded in the study area.
With the estimates of the intensity function at each pixel v i , i = 1, … , N probabilities of fire occurrence at each pixel may be approximated using the standardization .
The area-interaction point process model was selected because its flexibility to capture different degrees of clustering and repulsion.The point process model for wildfires in Valencia included the covariates lithology, combustion index (Figure 3), land use, orientation, elevation and slope (Figure 4).Important covariates such as rainfall and wind speed and direction were not available for the study area and were not included in the model.

Simulation of the fire spread process
We simulated the process of fire propagation using simulated ignitions as follows: Identify the nonburning neighbors of F(t + 1) and repeat steps 2 and 3 for all of them and update the pixels in the fire front F t+1 5. Repeat step 4 until no extra pixels burn after K iterations or K = K max , where K max is the maximum number of iterations fixed by the user.In our case, K max = 200 The simulation scheme described above corresponds to the so called importance sampling, modified so that only unburned pixels neighboring the fire front are candidates to become burned.In applications, the point of origin of a wildfire is known so step 1 above may be switched to pick the wildfire origin.In most cases, fire spread is also associated to ground wind speed, which induces some degree of anisotropy in the spread process.Such anisotropy may be incorporated in the simulation algorithm by giving more weight to pixels of V t in the downwind direction of the burning area B t ∪ F t (Boychuk et al., 2009).The above described algorithm was implemented using a code written in R-language.
To asses model performance we simulated random polygons using the following scheme: 1. Select 100 points of an ellipse centered at the origin pixel v in the previous simulation process with fixed main and second axes and orientation.The orientation of the resulting ellipse will be considered the anisotropy direction.2. Add a random noise to the coordinates of the selected 100 points of the previous step and take a sample of 15 of the those points.This number of points provided a nonregular polygon that resembled the shape of the burned area caused by wildfires.3. Close the resulting polygon and use it as the simulated burned area.4. Locate the pixel of origin of the simulated fire near the major axis of the random polygon, upwind to the wind direction to be simulated.
Note that the way we simulated the fires is independent of the fire spread model.We simulated four random polygons corresponding to the orientations SW-NE, S-N, SE-NW, and E-W, all using the same random centroid.For each polygon, we simulated 49 realizations of the fire spread model for the two possible wind directions (e.g., N and S for the S-N polygon, and so on) with maximum time t = 200.For each t ∈ {1, … , 200} the proportion of times a given pixel was classified as "burned" estimates of the probability of fire spreading to it from its neighbors.A detailed description of the simulation scheme is shown in Table 2. Using the estimated burning probabilities for each pixel, for a given threshold value, a pixel may be classified as burned if  i >  and unburned otherwise.For the N pixels in the area of interest , the classification process may be summarized in a classification or confusion matrix (Fawcett & Tom, 2006; Table 3): With the fixed threshold the True Positive and False Positive rates (TPR or recall, and FPR, respectively) defined as (Egan, 1975) Using those quantities we computed ROC curves and the area under the curve (AUC) statistic (Egan, 1975).We also computed the precision and recall statistics, defined as where TP is the number of true positive (burned) pixels, FP is the number of pixels that were incorrectly classified as burned when they were not and FN is the number of pixels classified as nonburned when they actually burned.Precision is a measure of the proportion of burning pixels correctly identified, while recall is a measure of how many truly burning pixels are returned by the model.Precision-recall is useful to measure of success of binary prediction when there exists imbalance between the two classes.We present the results of this statistics in the form of precision-recall curves (PRC).

RESULTS AND DISCUSSION
To estimate the probability of fire ignition in the whole Valencian Community we used the methods described in Section 2.2.The results of the Area-Interaction point process model fitted to the point pattern of fires occurring in the Valencian Community between 2006 and 2015 are shown in Table 4. Except for the coefficients associated to sparsely vegetated area and slope, the rest of the linear term coefficients in the model for the conditional intensity function were asymptotically significant.Compared to the base category (urban areas) where fires in vegetated areas are rare, the values for the conditional intensity function and hence the probability of a wildfire ignition, is high in the rest of the land use categories.The highest probability of ignition occurs in areas covered by vegetated low-lying areas on regularly flooded soil.Elevation and aspect have a small negative effect on the conditional intensity function.The combustion index on the other hand has a positive effect.The resulting probability of wildfire ignition for the Valencian Community is shown in Figure 5.
To compare how the model behaves under isotropy with spatial variation in the fire probabilities against constant fire probability, we carried out 49 simulations with the same pixel of origin, using the fire probabilities obtained from the intensity function fitted to the observed point pattern a 67 × 67 pixel size square within the Valencian Community.The time interval for the spreading simulation was 200.The results are shown in Figure 6, where it may be seen that the shape of the burning probabilities under isotropy and spatial variation in the fire probabilities behave different when fire spreads in an environment with uniform fire probabilities.Thus, when the same anisotropic weights are used in the simulation process, fire spreads following the spatial pattern of the fire probabilities, as one would expect under the hypothesis that fire spreads depending on the fire risk of the pixels neighboring the fire front.
Figure 7 shows the results of the simulation of wildfire spreading on a random polygon within the Valencian Community using the anisotropy weights shown in Table 2.The time interval is [0, 200] and the number of realizations was 49.The figures represent the proportion of times a pixel was classified as burned at t = 200, which is equivalent to the   The burned areas resulting from the simulations are irregularly shaped and different amongst them, following patterns similar to those observed for the burned areas by real wildfires.
The simulated fire spread patterns do not necessarily matched the anisotropy angle selected, because some of the pixels in the fire front become neighbors with a zone where fire risk is low, this is, where the intensity function of the point pattern used to compute the fire risk is very low.In those cases the fire spreading in the anisotropy direction slowed and the result was a burned area different to the polygon simulated as the fire shape.
In Figure 8 we show the ROC curves and the AUC statistic for the eight images shown in Figure 7 In all cases the empirical ROC curves are well above the 1:1 line, showing that the fire spread model proposed in this works acts well as 1099095x, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2774by Universitat Jaume I de Castellon, Wiley Online Library on [28/10/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License a classifier of burned-nonburned pixels in the study area which corresponds random classification.In the context of the fire spread prediction the ROC curves as well as the AUC statistics show that the proposed fire spread model has a good classificatory performance and is a reliable predicting tool under the conditions of the simulation analysis.ROC curves have been used widely to test the quality of binary classifiers in medicine and machine learning problems.Despite their usefulness, ROC curves are not the best option to assess the quality of binary classifiers when data show imbalance, as is our case given that we have a high proportion of unburned pixels in the study area.In Figure 9 we show the PRC curves for the 8 simulations used to test the model.The PRC curves show that if one is willing to accept a sensitivity level of say, 80%, then the precision of the model we propose is above 80% in all cases.Thus, we may consider that the model shown here may be considered in the range of good to excellent classifier.This classification is consistent in terms of early retrieval, this is, the area under the curve the right quarter of the ROC curves (Saito & Rehmsmeier, 2015).

CASE STUDY: APPLICATION TO TWO WILDFIRES IN THE VALENCIAN COMMUNITY, SPAIN
During the last four decades, wildfires have become more frequent in the Mediterranean area of Europe (Moreira et al., 2011).The spatial pattern for wildfire occurrences in forests around the world has been related to external variables of the two wildfires and the estimated probabilities of fire inside each square.The plots at the center and right of Figure 10 are close ups of the squares around the Artana and Carcaixent 2016 wildfires.The true limits of the burned area for each wildfire are also shown with black lines.In the figure we may see that wildfires spread over areas with locally high probability of fire.The actual shape of the burned area however does not follow a circular-shaped pattern perhaps as a consequence of firefighting activities and factors such as winds and presence of fuel materials and different combustion index at each pixel.
In Figure 10 the upper panel of plots presents the burning probabilities obtained after 49 simulations of the fire spread model with t = 200.The corresponding ROC curves and AUC statistics are shown in the lower panel of the Figure 11.It can be seen that part of the predicted burned area (green) lies outside the actual perimeter of the Artana and Carcaixent wild fires.As mentioned earlier, this result in part due to the firefighting activities that block the limitless spread of a wildfire.However, notice that the final size of the simulated wildfires is bounded in all cases.The ROC curves for the simulated fire spread in Artana and Carcaixent show that the model has a good classificatory performance for the real data and we may say that the model proposed is adequate to describe the spreading pattern of these particular wildfires.This observation is also consistent with the PRC curves for the two wildfires shown in Figure 12.This process is similar to the one developed by (Boychuk et al., 2009), (Johnston et al., 2006), but in our case we do not rely on an autologistic approach.
The model proposed in this work has the advantage of being simple and easy to use.Incorporating factors that control fire spreading can be done by assigning appropriate weights for each direction.So if for instance wind speed increases in one direction at some section of the fire front, then that weight can be increased and run the model with the new weights.Also, if firefighting has been successful in some sections of the fire front, weights pointing to those direction can be modified accordingly.
The model has limitations as it does not include covariates explicitly.However, relevant covariates can be incorporated in the point process model fitted to obtain the fire risk at the pixel level, so covariate influence on fire spread is implicitly incorporated in the spatial fire risk estimation process.Although limited, the model is useful to obtain quick estimates of the potential burned area and to highlight areas where fire control efforts may decrease the rate of spreading of wildfires.

CONCLUSIONS
We have presented a simple model and an algorithm to simulate wildfire growth and spread using spatial point process methods to obtain probabilities of ignition in different locations at a pixel level.The fire spread model can be seen as a Markov transition process in the sense that there is a state space (burned, unburned, fire front) and transition probabilities.Albeit simplistic, the model provides a good approximation to the burned area and shape of the wildfires.In addition, the inclusion of covariates in the model improves the quality of the estimated intensity function used to obtain the local probabilities of fire used in the simulation and the detection of fire growth.
Finally, analysis and development of the random spread process methodology can be carried out in other fields in which the variables to be studied have a spatial growth in certain areas.The clearest example of the usefulness of this methodology presented here would be the study of epidemic spread in a geographic area.The inclusion of the variation in the study area of different covariates such as elevation, distance to fixed elements such as pollution sources, may improve the prediction of the epidemic spread.
right) in the Valencian Community.The small green dots show the wildfires occurring in the Valencian Community between 1993 to 2015.F I G U R E 4 Landuse, orientation, elevation and slope for the Valencian Community.The small green dots show the wildfires occurring in the Valencian Community 1993 to 2015.
Note: *, ** and *** means: significance at the 10, 5 and 1%.4e-08 6e-08 8e-08 1e-07 1.2e-07 F I G U R E 5 Spatial distribution for the probability of wildfire ignition in the Valencian Community, Spain.The values in the map were obtained from an Area-Interaction spatial point process model, fitted to the wildfire occurrence pattern observed between 2006 and 2016.
Spread patterns of simulated wildfires in a 67 × 67 pixel square in the Valentian Community using the same pixel of origin for the fire.The upper panel shows how fire would spread when fire probabilities are completely random in the square.The lower panel shows the pattern of fire spread under isotropy and constant fire probabilities in the square.In both cases the images show the burning probabilities obtained from 49 realizations of the fire spread model for a time interval from 1 to 200.number of times the fire reached a given pixel.The different shapes and sizes of the simulated wildfires highlight the spatial variability of the burning probabilities in the Valencian Community, consequence of the spatial variation of the covariates incorporated in the model to estimate the intensity function of the point pattern of fires in the study area.
Random polygons representing the burned area by wildfires (black lines) and burning probabilities obtained with 49 realizations of the fire spread model using the probabilities of fire obtained from the intensity function of the fire pattern observed in the Valencian Community between 2006 and 2016.ROC curves corresponding to the burning probabilities obtained with 49 realizations of the fire spread model using the probabilities of fire obtained from the intensity function of the fire pattern observed in the Valencian Community between 2006 and 2016.
U R E 9 PRC curves corresponding to the burning probabilities obtained with 49 realizations of the fire spread model using the probabilities of fire obtained from the intensity function of the fire pattern observed in the Valencian Community between 2006 and 2016.Close up for the zone of the Valencian Community where the wildfire of Carcaixent occurred in 2016.The left plot shows the estimated fire probabilities obtained by fitting an area interaction point process model to wildfires between 2006 and 2016, as well as the squares used in the simulation process and the actual shape of the Artana (Northern square) and Carcaixent wildfires (Southern square).Center plot: close up of the square for the Artana region, showing probabilities of fire and the actual polygon of the 2016 wildfire.Right plot: close up of the square for the Carcaixent region, showing probabilities of fire and the actual polygon of the 2016 Carcaixent wildfire.
1099095x, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2774by Universitat Jaume I de Castellon, Wiley Online Library on [28/10/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Land use in Castellón and codes assigned Heymann et al., 1994.ategoriesrepresent an aggregation of a more diverse categorization for land use published byHeymann et al., 1994.1099095x,0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2774by Universitat Jaume I de Castellon, Wiley Online Library on [28/10/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Origin coordinates, anisotropy weights and wind directions selected for the simulations of the random spread model proposed, in a random polygons within the Valencian Community, Spain Classification matrix obtained for a fixed probability threshold value used to classify pixels as burned or unburned after 49 simulations of the fire spread model with t = 200 TA B L E 2Note: The eight random polygons shared the same centroid.TA B L E 3 1099095x, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2774by Universitat Jaume I de Castellon, Wiley Online Library on [28/10/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Point estimates and confidence intervals for the parameters of the Area-Interaction point process model fitted to the observed wildfires occurring in the Valencian Community between 2006 and 2015 1099095x, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2774by Universitat Jaume I de Castellon, Wiley Online Library on [28/10/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Precision-recall plots for the wildfires of Artana and Carcaixent in 2016.The plots show that for sensibilities greater than 80%, the precision of the model proposed is well above 80%.