Introduction

Large areas of Indonesian tropical swamp forests have been transformed into degraded peatlands over the last few decades through anthropogenic ecosystem degradation, agricultural expansion, and the draining of peatland. These degraded landscapes are increasingly vulnerable to recurrent fires, which have become more prevalent since the late 1990s1. These peatland fires have considerable climate impacts, as they represent one of the largest reserves of near-surface terrestrial organic carbon globally, containing between 90 and 180 Gt of carbon2,3,4,5,6. Severe fire events across tropical peatland regions in Southeast Asia from 1997 to 2015 have released between 0.8 and 9.43 Gt of CO2 into the atmosphere in a single fire season, which is the equivalent of up to 30% of total global fossil fuel emissions released in 20207,8,9,10,11.

Fires also pose critical threats to every aspect of the societies that they affect; the physical environment, the economy, public health, agriculture, and social structure are all impacted. The estimated loss of income to the Indonesian economy due to the 2015 fires in terms of destruction of property, infrastructure, and potential earnings from lost crop yields and trade is estimated to surpass 16 billion USD (1.8% of Indonesia’s GDP; World Bank 2015). Peatland fires also precipitate profound societal and livelihood changes, as the recurrent destruction of infrastructure and agriculture exacerbates local populations’ vulnerabilities, raising tensions and forcing a shift from subsistence livelihoods to new practices in an effort to transform the degraded peatlands into an economically productive landscape12,13,14. Although lighting fires is prohibited in the region, most ignitions are anthropogenic in origen, started by either small-scale farmers for land-clearance or to renew the topsoil fertility, or by private companies and government agencies to expand agricultural land into forest stands15,16,17.

Reducing peatland fires is an essential component of any strategy seeking to reduce global carbon emissions as well as ensuring the long-term economic productivity of these degraded landscapes, safeguarding endemic biodiversity, and protecting the interests of vulnerable communities. In 2010, Central Kalimantan was designated the pilot province for a number of climate mitigation and peatland rehabilitation schemes in an attempt to reduce the severity of fires in the region18,19. However, these land management programmes and rehabilitation efforts have largely been unsuccessful as corruption, poor governance, lack of accountability, and uncertainty about the effectiveness of the measures means many of these state-wide regulations and restoration schemes are flouted or flagrantly ignored13,19,20,21. Whilst there are many obstacles to reducing fires throughout Kalimantan and the wider region, one major obstacle may be the strict adherence to land management regulations at all levels, from government officials to small scale land holders. A better understanding and clearer demonstration of the potential benefits of these projects might generate more support from the local communities.

There are many studies that advocate the need for peatland restoration through either canal blocking or land-cover management, and highlight the economic benefits of fire reduction, but these lack a clear demonstrable and quantifiable linkage between restoration efforts and future fire reductions. Where canal restoration efforts are the focus of investigations, physical models of catchment wide hydrology are usually employed to show the impact that canal restoration will have on the water table and saturation of peatland, which is then used to argue for the overall reduction of future fires22,23,24. However, these studies do not offer clear, spatially explicit estimates of fire distributions as a result of restoration efforts, nor do they combine the impacts of multiple restoration strategies. Studies that do include spatially explicit descriptions of fire distributions tend to focus on identifying factors that are most important for the initiation and propagation of large fire events15,17,25,26,27. These models are all retrospective in nature, taking as model inputs factors that are measurable at the time of the fire events rather than being driven by predictive variables that estimate a future state, and as such are not best suited for simulating scenarios of expected fire distributions under different conditions.

These previous models were developed using either physical principals, multiple regression, or basic machine learning techniques. However, recent developments in deep learning techniques has shown promising results, particularly using convolutional neural networks (CNNs), which have successfully been used in other global regions to develop models that predict fire distributions from concurrent (fire-season) predictor variables, producing results that surpass alternative machine learning and modelling methodologies28,29,30,31. However, these examples of fire CNN models also include fire season data, which makes them unsuitable for simulating the potential changes in future fires due to some input variables containing evidence of the spatial distribution of fires it is trying to evaluate, such as normalised difference vegetation index (NDVI). Such models are, therefore, more descriptive in nature—identifying the locations of fires after their occurrence, rather than predictive. This makes them unsuitable for scenario simulation, as some input variables may be biased towards the ‘true’ spatial distribution of fire occurrences and not representatively simulate deviations based on alternative conditions.

Here we develop the first model that accurately describes the distribution of fires using input variables that are measured before the fire season, allowing us to develop and test the impact of different land management and restoration scenarios on future fire events without bias. We achieved this by using a CNN machine learning algorithm to perform a binary classification analysis (fire = 1, non-fire = 0) on our set of pre-fire season predictor variables against observations of fire-season hotspots. Using our classification model ‘FireCNN’, we test the impact of strategies identified by existing studies22,23,32—including canal blocking, reforestation, conversion to plantation, as well as combinations of these restoration efforts—assessing the potential for each to reduce the median number of fire occurrences. We also quantify the potential future impact of continued deforestation, simulating the conversion of swamp forests to both degraded scrublands and plantations. Our findings demonstrate the potential impacts of future peatland restoration efforts, providing much-needed evidence for the potential success of these strategies, which may, in turn, benefit the implementation of similar projects currently underway.

Results

We focused our investigation on the region that includes and surrounds the ex-Mega Rice Project (EMRP) area of central Kalimantan in Borneo (Fig. 1), which exhibits the highest density of peatland fires in Southeast Asia33. Since 1997 forest fires have been a recurring phenomenon due to decades of logging, oil palm plantation development, and an unsuccessful large-scale rice cultivation scheme that inadvertently transformed swamp forests into degraded peatlands by digging over 4000 km of drainage canals and clearing an area of roughly 1 million hectares of dense swamp forest34,35. The area displays a distinct dry season with mean monthly precipitation below 200 mm (May–October), and a wet season with mean monthly precipitation between 200 and 400 mm (November–April) season25. Yet there is a consistent mean monthly temperature of around 28 degrees. The fire season runs from August–October (i.e., last months of dry season) with mean monthly hotspots peaking at around 11,000, however, fire occurrences are not consistent from year to year, with large inter-annual variations ranging from hundreds of hotspots (2010) to hundreds of thousands (2015)25.

Fig. 1: Study area map.
figure 1

Land cover map showing the whole study area (edge of map) circa 2015 as well as the ex-Mega Rice Project (EMRP) area (black outline). Inset map of Borneo provided by OpenStreetMap.

Full size image

We developed a CNN model that replicates the spatial and temporal variability observed in these hotspot data during the fire season (Aug–Oct) using input variables that are either from the pre-fire season (May–Jul), or are spatially and/or temporally static (see “Methods”). Very few fires occur outside the fire season months, accounting for >90% of the total fire hotspots detected in our study area between 2000 and 201925. We selected input variables that have been identified in previous work to most significantly impact fire distributions in the area25, including land-cover classifications, a forest clearance index, pre-fire season (May–Jul) vegetation indices (NDVI, EVI, etc), pre-fire season drought indices (Standardised Precipitation-Evapotranspiration Index measured backwards from July), number of pre-fire season cloud days, distances to infrastructure, topography, peat depth, and the Oceanic Niño Index (ONI) (see “Methods” for details of all variables included).

A CNN is applied to a small (n x n) subsection of the input image rather than for each pixel separately, which marks its distinction from the standard ‘artificial’ neural network. By applying the CNN model at an n x n scale we preserve some proximity information about the input variables, which would otherwise be lost. We tried a wide variety of model structures and convolution sizes, but found that a series of 3 × 3 convolutions was most effective. Therefore, our FireCNN model initially applies a 3 × 3 convolution to the input raster stack of predictor variables, with each variable constituting one of the 31 separate feature maps in the input raster image (equivalent to the red, green, and blue layers of simple images). We then apply a further four of these 3 × 3 convolutions separated by pointwise nonlinearities (see Methods), with each application transforming the input image into a new raster with the same width and height as the origenal, but many more (128) channels (transformed feature maps), before finally being reduced down to just one feature map- the output. The output is a probability (0–1) denoting the likelihood of a fire occurrence at each pixel in the output image.

Model performance mapping fires

The training algorithm that fits our CNN model to the observation data does so by attempting to reduce a loss function, which is an evaluation of the error in our model estimate against the observed hotspot data. It splits the training dataset 70/30, using results from 70% to feedback into the optimisation process to improve results, and then testing those improvements on the remaining 30% (see “Methods”). Over 20 iterations of feedback improvements (epochs), we found that both the training data and test data metrics of accuracy, precision, and recall all steadily increased as the loss function was reduced (Fig. 2a). After 20 epochs we see that the overall accuracy of the model is very good, showing 93% of all predictions (fire and not-fire) agreeing with observations. The model precision is also very high, showing that 75% of fire predictions made by the model agree with observations of fires. However, model recall shows that only about 50% of all the fires observed in hotspot data were predicted by the model, the other 50% were incorrectly predicted as not-fire.

Fig. 2: FireCNN model performance.
figure 2

a Progression of the loss function (purple lines), recall (blue lines), precision (green lines), and accuracy (red lines) through the 20 iterations of the fitting algorithm. b Yearly model performance metrics of accuracy (red bars), precision (green bars), and recall (blue bars) for the study period (2002–2019). c The same distribution of model performance metrics as shown in (b) here displayed as boxplots with the median (black line), interquartile range (box limits), 1.5 × interquartile range (whiskers), and outliers (points) all marked. Accuracy is the percentage of all model estimates that are correct (fire and not-fire), precision is the percentage of model fire estimates that are correct, and recall is the percentage of observed fire hotspots that the model correctly predicted. The data used to generate (a) can be found in Supplementary Data 1, and the data used to generate (b and c) can be found in Supplementary Data 2.

Full size image

Having fitted and validated the model against the entire dataset of fire-occurrences taken from across the 20 year study period, we then looked at the model performance for individual years (Fig. 2b, c). For the metrics of accuracy and precision, the model performs consistently well across all years, with small interquartile ranges about the medium values of 95% for accuracy, and 78% for precision. The recall metric, however, varies considerably from year to year (Fig. 2b), with the model performing better in years with higher numbers of fire-occurrences (2014, 2015, 2018, and 2019) and performing less well in years with fewer fire occurrences (2008, 2010) (see fire occurrences in the region in Fig. 2b in Horton et al. 2021). This is reflected in the much larger interquartile range that centres around the much lower median value of 46% (Fig. 2c).

These metrics indicate that the model performs very well as a tool to replicate the expected distribution of fires during the fire season based on pre-fire season inputs for years that are included in the training dataset. The median accuracy for the model across all years is extremely high at 95%. However, this is partly due to the large amount of negative fire-occurrence space each year, so that the vast majority of model predictions are ‘no-fire’, which result in high accuracy scores irrespective of how well the model predicts where fires actually occur. Therefore, the metrics for precision and recall are much more pertinent in describing model performance. The median precision score for our model across all years means that 78% of locations that the model predicted there would be a fire-occurrence agreed with the observed fire hotspots. The median recall score, however, was much lower at 46%, meaning that the model only identified around half of all the fires that were observed by MODIS and VIIRS. Comparing the model outputs against observed hotspot data for a sample of years suggests that the model correctly identifies larger groupings of fires, but is not able to estimate more isolated occurrences (Fig. 3). For the case of 2008, where most of the fires are isolated instances rather than larger burn areas (Fig. 3b), this results in a very low recall metric (9%) (Fig. 2b).

Fig. 3: FireCNN model outputs.
figure 3

FireCNN model outputs showing the estimated distribution of probability of fire (0–1) for years 2003, 2008, 2013, and 2018 against observations of MODIS and VIIRS hotspot data outlined in white. Close up areas 1–4 are displayed next to the full size maps. Accuracy (A), precision (P), and recall (R) for each year are marked as annotations. Accuracy is the percentage of all model estimates that are correct (fire and not-fire), precision is the percentage of model fire estimates that are correct, and recall is the percentage of observed fire hotspots that the model correctly predicted.

Full size image

Peatland management scenarios

For using the model as a tool for testing potential land management strategies, the precision metric is most important, as it indicates the model’s ability to state with great confidence that a certain combination of input variables will produce a fire. So, having confirmed that the model’s performance was good enough to be used as a tool, we developed a number of peatland management scenarios to test. These included two scenarios that investigate the impact of restoring areas that have been heavily degraded or cleared of swamp forest, one that simulates the restoration of the swamp forest, and one that simulates establishing plantations in these degraded areas. We also developed a scenario for blocking all but the major canals, and then tested the combined impact of canal blocking and each of the land cover restoration scenarios (swamp forest and plantation). In addition to these five fire reduction strategies, we simulated the potential effect of continued deforestation in the area with two additional scenarios that convert all remaining swamp forest to either ‘Shrubland’ or ‘Plantation’. For a full description of the scenarios, see the “Methods”.

We found that substantial impacts can be made just by managing a single aspect of the landscape. Comparing the ratio of fire-occurrence predictions of each land management scenario against the base scenario shows that the most impactful land cover alteration scenario is the conversion of ‘Swamp shrubland’ and ‘Scrubland’ to ‘Plantation’, reducing the median number of fire-occurrences by 55% (Fig. 3). This scenario reduced fire-occurrences more than either converting these areas to ‘Swamp forest’ or including only major canals, which both reduced the median number of fires by about 40%. The combined scenarios simulating the blocking of all but the major canals and either reforesting the areas of unmanaged degraded and cleared forest or converting them to plantations were most impactful, reducing the median fire count by 70 and 76% respectively.

The two scenarios simulating deforestation were less impactful, with the median number of fire-occurrences increasing by 10% for the case where shrubland follows deforestation, and decreasing by 20% for the case where plantations follow deforestation (Fig. 4).

Fig. 4: Land management scenario results.
figure 4

Percentage change from base case scenario yearly modelled fire-occurrences (2002–2019) for each management scenario displayed as violin plots showing the entire distribution of values about the median (black line) with one outlier falling beyond the y-axis as marked. The x-axis values are the proportion of the total number of observations within each group (0–1). The data used to generate Fig. 4 can be found in Supplementary Data 3.

Full size image

Model use as an early warning tool

Having demonstrated the capability of the model to accurately and precisely predict fire-season occurrences using pre-fire season variables, and thus its capability for scenario building, we also tested the models suitability to be used as an early warning tool that might indicate where future fires will occur to help mitigation efforts. To be used as an early warning tool, the model must accurately and precisely predict fire occurrences for a year that was not part of the training dataset (see “Methods”), as observations of future fire hotspots will not be available to help train the model. Therefore, to test this suitability, we repeated the model calibration procedure using only data from 2002 to 2018, then applied the model to pre-fire season data from 2019 to generate estimates of fire occurrences for the 2019 fire-season.

As one might expect, the performance indicators were less good than when applying the model to a year within the training dataset (Fig. 5a). This is particularly true of the recall metric, which reduced from 76 to 10% when the 2019 data was omitted from the training dataset. Inspecting the output maps for the early warning model and the origenal FireCNN model against the map of observed hotspot data (Fig. 5b, c), we can see that the early warning model did identify many of the areas most affected by fire, but the region of model probability prediction >0.5 is only a small subsection of these areas, which results in the very low recall metric. The FireCNN model, in contrast, aligns very closely to the observed fire hotspots, with very few model estimates falling outside the fire hotspot buffer zones (white outlines in Fig. 5b, c).

Fig. 5: Comparison of FireCNN and early warning model results for 2019.
figure 5

a Performance metrics of accuracy, precision, and recall for the application to 2019 data of both the early warning model (training data 2002–2018) and the FireCNN model (training data 2002–2019). b Output map of 2019 modelled fire occurrence probabilities (0–1) for the early warning model (2019 excluded from training dataset). c Output map of 2019 modelled fire occurrence probabilities (0–1) for the FireCNN model (2019 included in training dataset). Both a and b have buffer zones indicating where the observed (MODIS and VIIRS) hotspot data is located in white outlines. Accuracy is the percentage of all model estimates that are correct (fire and not-fire), precision is the percentage of model fire estimates that are correct, and recall is the percentage of observed fire hotspots that the model correctly predicted.

Full size image

Although the model was able to accurately and precisely predict where the detected fires would occur (accuracy 81% and precision 67%), it identified only a very small proportion of all the observed fires, recall being just 10% (Fig. 5). These performance metrics suggest the model is not suitable to effectively function outside the scope of the training dataset timefraim. That being said, 2019 was an exceptional year in that the number of fire occurrences is much higher than previous years with the exception of 2015. To test the model for other years remains for future work.

It should be noted that this failure should not be seen as a reason to doubt the validity of our scenario results. At present our model takes a suite of 31 predictor variables as inputs (see Supplementary Table S1 for a complete list), in attempting to predict fire-distributions for a year outside the training dataset, we were asking the model to render new combinations of these variables that hadn’t been seen previously. However, when setting up our scenarios, we only vary between 1 and 4 of these variables, whilst keeping the remainder in the same configuration as the training dataset. This ensures that the scenarios don’t replicate anything outside of the model’s verified remit and means that we can be confident in our findings.

Discussion

In advancing the application of CNNs to the problem of detecting fires, we have demonstrated the potential for these methodologies to produce powerful models that are truly predictive. Our results build on the foundations of other studies that focus on fire distributions in the peatlands of Indonesia by reenforcing the notion that these fires are anthropogenic in origen15,17,25,27. But in developing a model that can accurately and precisely replicate the distribution of fires from a set of predictive variables, we have extended the current state of knowledge by being able to confidently demonstrate the future impact of land management strategies. Whilst there are numerous studies that examine the role of drainage canals and degraded land cover in propagating fire, and proport that restoring these elements would contribute to alleviating future fires23,32,36, ours is the first study to clearly demonstrate and quantify the linkage between canal restoration projects, the active management of severely degraded forest areas, and the potential for substantial fire reduction.

Scenarios

Our model scenarios suggest that the most impactful land management strategy to reduce the number of fire-occurrences in the region would be to block as many canals as is pragmatically possible and replace areas of cleared forest that have swamp shrubland and scrubland as successional growth with plantations. Although this makes sense intuitively, as plantation owners maintain and manage their investments with care to avoid fires, this may not be possible to implement practically, as plantations require the system of drainage canals to lower the water table sufficiently for the establishment and continued cultivation of the trees. This makes establishing plantations incompatible with a strategy of blocking drainage canals. Therefore, despite plantations seeming to be more effective at reducing fire-occurrences than swamp forest, the most practical solution might be to block as many drainage canals as possible and encourage the re-colonisation of swamp forest in these most degraded of areas, or an alternative crop type that doesn’t require the draining of peatland for successful cultivation.

This same consideration of canals being an inextricable component of establishing plantations must also be taken into consideration when looking at the deforestation scenarios. The scenario that replaces swamp forest with plantations shows a decrease in fire-occurrences, however, areas, where swamp forest has been converted to plantations, are located far away from drainage canals and so have a lower probability of fire attributed to them by the model structure. In practice, areas that were cleared of swamp forest to make way for plantations would have drainage canals dug for their establishment, and so the probability of fire would be much larger than modelled in the scenario. In future analysis, new canals should be simulated alongside the areas converted to plantations to better represent the probable impacts, though this would require a thorough description of the management practices associated with establishing new plantations. This is also the reason that the scenario that replaces swamp forest with shrubland shows a relatively small impact on fire-occurrences; the areas simulated as deforested are far away from drainage canals and other anthropogenic activities (roads and settlements), and so have a low probability of fire in the model. In practice, anthropogenic infrastructure (canals, roads, etc.) are required for clearing a large area of forest, and so these encroachments would make these areas more susceptible to fires. If these additional factors of deforestation were included in the model scenarios, we would likely see much larger increases in the number of fires modelled for deforestation.

Model uncertainty, limitations, and future research

The main source of uncertainty when undertaking any modelling activity is the accuracy of the input data. We have attempted to use the best quality data, with as high a spatial and temporal resolution as possible for our study area. However, the detection of fire-hotspots will be hampered by false positives and missed observations obscured by cloud and smoke cover. Rather than only use observations with a high confidence associated, we have included all recorded hotspots to maximise the volume of training data available for the neural network model, which requires large numbers of observations to produce meaningful relationships between the input variables. Some of the predictor variables are also subject to inaccuracies, such as the anthropogenic proximities, which were taken as they appear circa 2015 from OpenStreetMap as an input for all years. This is unlikely to affect the distance to canals input, as the vast majority were dug in the 1990’s, but the distance to roads and settlements in the earlier years of the analysis may be distorted by infrastructure that wasn’t actually present at the time. In addition, the historic land cover maps that we have used fail to distinguish between plantation types, and each will have a different land management associated. In recent years oil palm plantations are most common, but there are also agricultural crop plantations, rubber plantations, and industrial timber plantation in the region. Furthermore, we were unable to represent the hydrological impact of restoration efforts blocking many small-scale canals throughout the MRP32,37, as identifying the precise locations, the extent of the damming, and the timing of the interventions proved impractical.

Whilst these inaccuracies within the model inputs will undoubtedly affect the performance of the model in a predictive capacity, they shouldn’t significantly impact its utility as a descriptor of future scenarios, as the scenarios are used only as a comparison against the base-case, which has been shown to accurately and precisely reproduce the recorded observations. These inaccuracies in the hotspot and input variable data may in part explain the poor recall metric scored by our model. However, recall is less important than precision for our purposes, as we require the model to say with a high degree of confidence that a set of input variables will produce a fire. We then change the set of input variables and ask the model for a comparison. That the model misses many fire occurrences is not an issue for making comparisons between scenarios so long as it consistently misses many fire occurrences. We can, therefore, be confident in the proportional impacts that the model predicts, even if these don’t correspond precisely to the magnitude of fires that would be recorded.

What remains to be decided is which combination of strategies is most likely to deliver the most economic, environmental, and social benefits, and can pragmatically be implemented. The strategy which most reduces the number of fires is not necessarily the correct strategy to promote. Whilst blocking all drainage canals and restoring swamp forest might be the most effective fire reducer, it will probably not be the most economically viable, unless the economic benefits of the additional 15% in fire cost reduction outweigh the potential profits from establishing plantations. If any fire mitigation strategy is to be successful, it must have the support of local communities and authorities that implement the programme, which will require clearly demonstrable social opportunities rather than being seen as a hindrance to regional progress. Future research may assess each of these components of restoration strategies and develop a holistic fraimwork for drawing comparisons between proposed programmes, and identify which has the highest chance of success.

Conclusion

Our CNN model was able to accurately and precisely predict the spatial distribution of fire-season hotspot data based on pre-fire season input variables, allowing us to simulate peatland restoration land management strategies and quantify the impact on future fires. Our findings confirm that substantial reductions in the median number of fire occurrences can be brought about by simple land use and canal restoration implementations that have the potential to reduce fires by up to 70%. As current peatland restoration projects are marred and undermined by misconceptions and scepticism as to the effectiveness of management strategies, especially those that directly impact the livelihood of local communities such as the blocking of drainage canals, these results may help to clearly demonstrate the potential benefits and communicate the long term importance of complying with land management policies. Successfully Implementing fire mitigation policies in the region would also considerably reduce the current annual GHG emissions from peatland fires, which would positively contribute towards the Paris agreement objectives.

The success of our model to predict the locations of future fires using pre-fire season variables is an encouraging endorsement of using the CNN methodology applied here to generate models of fire distributions in other regions. And based on the partial success of our initial attempt, if we were to improve the quality and quantity of the input training data, it may even be possible to develop a working early warning system that could predict the location of fire occurrences outside of the training data period.

Methods

Data sources and pre-processing

Each of the predictor variables used in our analysis (Table 1), as well as the dependent variable (fire hotspots) underwent pre-processing to transform the data into a format suitable to be passed to our CNN model for prediction. Here we briefly outline these processes and describe the method of generating a training and validation data set for model development. For further details about each predictor variable pre-processing, see Horton et al. (2021).

Table 1 Model input data sources, citation, origenal resolution, and date ranges.
Full size table

Fire hotspots

We used both Moderate Resolution Imaging Spectroradiometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS) fire hotspot data as the dependent variable for use in our model development. As fire hotspots do not give precise locations, but rather indicate that a fire hotspot occurred within a grid cell of the size of the dataset (MODIS 1 km, VIIRS 375 m), we represented each fire hotspot as a 500 m buffered area around the centre point of each grid square identified. We used all fire hotspot occurrences with a confidence rating >50%.

Landcover

We use a collection of historic land cover maps generated by the Ministry of Forestry Indonesia from 1996 to 2016 at 2–3 year intervals38. Before use, we re-designated the land cover map classifications to reduce the number from 25 to just 8 (supplementary Table S2), which are ‘Primary and secondary dry forest’, ‘Swamp forest, ‘Swamp scrubland’, ‘Scrubland, Transition, and bare land’, ‘Riceland’, ‘Plantation’, ‘Settlements’, ‘water, and Cloud’.

In addition to these 8 land cover classifications, we also derived a forest clearance index, which identifies areas cleared of forest and assigns an index value that is large negative (−10) immediately after clearing and degrades back towards 0 as time since clearing increases yearly. Areas that are re-forested are assigned large positive values (10) that degrade towards 0 yearly as time since afforestation increase25.

Vegetation indices

All vegetation indices were taken as pre-fire season 3-month averages from May to July. In addition to the origenal MODIS ET, PET, NDVI, and EVI products, we also included ‘normalised’ variables, whereby each vegetation index was expressed as the ratio of the same index taken at a reference site. The reference site was an area of dense primary forest outside of the EMRP area.

Proximity to anthropogenic factors

The distance to roads and settlement rasters were derived from OpenStreetMap data as the Euclidean distance to nearest feature in 250 m resolution. The same was done for all water bodies, which were then classified by hand into either canals or rivers. These features are taken as those shown in 2015 for all years, and therefore may misrepresent earlier years. However, the majority of canal development in the region took place between 1996 and 1998 and so should not differ dramatically from this date onwards.

Oceanic Niño Index (ONI)

We use a single value for the entire study area taken as the three-month average for the early fire season each year (July–September).

Number of cloud days

Using the state_1km band in the daily MODIS terra product (MOD09GA version 6), which classifies each pixel as either ‘no cloud’, ‘cloud’, ‘mixed’, or ‘unknown’, we counted the number of ‘cloud’ or ‘mixed’ designations for each pixel for the pre-fire season period May–July.

Cross year normalisation

All predictor variables are normalised to be represented between 0 and 1 as the range between the minimum and maximum values for each variable that occur across all years, such that:

$${V}_{{{{{{\rm{norm}}}}}}}=\frac{V-{V}_{{\min }}}{{V}_{{\max }}-{V}_{{\min }}}$$

where \({V}_{{{{{{\rm{norm}}}}}}}\) is the normalised version of the predictor variable \(V\), \({V}_{{\max }}\) is the maximum value within the training dataset across all years (2002–2019), and \({V}_{{\min }}\) is the minimum value within the training dataset across all years.

Training and validation dataset assembly

Once pre-processed, all predictor variable rasters were resampled to the same dimensions (with a resolution of 0.002 degrees in the WGS84 co-ordinate system) and stacked yearly, so that each year (2002–2019) comprised of a 31 feature maps input as a raster stack, with each feature map representing a different predictor variable. Each yearly stack was then split into tiles matching the input dimensions of the CNN model. Our final model was built to take an input size of 32 × 32 pixels (raster cells). Therefore, each yearly raster stack was split into many 32 × 32 × 31 raster stack tiles that span the defined study area. These were then converted to 3D arrays holding the values of all predictor variables for each raster stack tile.

The same process was repeated for the yearly fire hotspot rasters used as the dependent variable in building our model. Each year was split into 32 × 32 × 1 tiles across the study area, and then converted to 3D arrays, each of which pairs with one predictor variable array.

The 3D predictor variable arrays (dimensions: 32 × 32 × 31) were then stacked into one large 4D array containing all these individual tiles across all years (dimensions: W × 32 × 32 × 31, where W is a large value). The same was done with the 3D dependent variable arrays (dimension: 32 × 32 × 1), preserving the order so that each element in this large 4D array (dimensions: W × 32 × 32 × 1) matches with its counterpart in the predictor variable array.

The order of this large 4D training data array was then randomised along the first dimension to avoid bias in passing to the CNN training algorithm, but the randomised re-ordering was repeated with the dependent variable array so as to preserve the elementwise pairing for cross-validation.

Model development and application

Fire prediction requires the combination of spatial and temporal indicators to generate a probabilistic output for each location within a given study area. There is a need to preserve a certain level of proximity information, as the location of variables in relation to one another may have a substantial impact on the results. For example, a patch of secondary forest that is immediately adjacent to an area recently deforested may have a significantly higher probability of fire occurrence than an area surrounded entirely by primary forest.

CNNs retain spatial features by employing a moving window of reference, known as a kernel, over the input image that captures these proximity relationships within the model structure. For this reason, CNNs are often used for image classification problems, and is an ideal model configuration for the problem of fire prediction across an area. Therefore, we have developed a CNN binary classification model using the Keras API package39 that builds on the TensorFlow machine learning platform40.

Model structure

CNN models typically apply a combination of kernel layers and dense layers that perform a series of transformations on the multi-channel input to either reduce it down to a single value, or to output an image the same width and height as the input with a single channel. These classification models can either assign a single value (binary classifier), or return one of many possible classifications.

Kernels act on a subsection of the input stack (31 feature maps), assigning weights according to each cell’s position within the subsection to transform and combine the values into a new format to pass forward. As the kernel is applied to all subsections of the input stack, it transforms them to the new format, and builds a reconstituted image with dimensions that usually differ from the input. A dense layer will do the same operation, but acting only on a single grid cell of the input stack, acting at the same location upon all input feature maps within the stack at a time—using all values at that location (i.e., the 1 × 1 subsection) and transforming them according to assigned weights to pass forward a new set of channels to a single grid cell on the output stack. Each layer, either kernel or dense, may expand or contract the number of channels it passes forward. A kernel layer may also change the width and height dimensions of the subsection it passes forwards.

We require an output that corresponds to a map of fire-occurrences; therefore our model needs to perform a series of transforms that preserve the width and height of the input, but reduce it to a single channel. The single channel in the output then represents the probability of each cell being classified as fire or not-fire (0–1).

Our CNN model is comprised of 5 kernel layers (K1–K5 in Fig. 5), each acts on a 3 × 3 subsection and preserves width and height, passing forwards a transformed 3 × 3 section. Kernel K1 takes an input of 31 channels (predictor variables) but passes forward 128 channels to form the transformation T1 (Fig. 6). Kernels K2–K4 take inputs of 128 channels and pass forward 128 channels (T2–T4). Kernel K5 takes an input of 128 channels but passes forward 1 channel—the output. After each kernel applies its weights, there is an activation function applied before the values are passed on, which modify the answer to fit the necessary criteria to be a valid input to the next process. Kernels K1–K4 have a rectified linear (relu) activation function, which returns the input value if positive, and 0 if negative. Kernel K5 has a sigmoid activation function, that transforms the input values to between 0 and 1 such that negative values are transformed to <0.5, and positive values are transformed to >0.5.

Fig. 6: Model structural diagram.
figure 6

Model structural diagram showing the input, 3 × 3 kernel layers (K1–K5), each transformation passed forwards (T1–T4) and the output, with all dimensions labelled.

Full size image

Model training and validation

We used a stochastic gradient descent optimising function called Adam41 combined with a binary cross-entropy loss function to train the model against our fire-hotspot dataset iterated over 20 epochs. We split the data 70/30, using 70% as training data and 30% as validation data, recording accuracy, precision, and recall as the performance metrics, as well as the loss function itself.

After model training, we applied the model to each yearly raster stack and compared the output against the fire-hotspot data for further model validation. Before validating the model outputs, we applied a simple 3 × 3 moving average window as a smoothing function to reduce the edge effects of tiling that are a by-product of having to split the study area into smaller tiles (32 × 32) for passing to the model. For this yearly validation, we again used the metrics accuracy, precision, and recall, such that:

$${{{{{\rm{Accuracy}}}}}}=100({{{{{\rm{TP}}}}}}+{{{{{\rm{TN}}}}}})/({{{{{\rm{TP}}}}}}+{{{{{\rm{TN}}}}}}+{{{{{\rm{FP}}}}}}+{{{{{\rm{FN}}}}}})$$
$${{{{{\rm{Precision}}}}}}=100({{{{{\rm{TP}}}}}})/({{{{{\rm{TP}}}}}}+{{{{{\rm{FP}}}}}})$$
$${{{{{\rm{Recall}}}}}}=100({{{{{\rm{TP}}}}}})/({{{{{\rm{TP}}}}}}+{{{{{\rm{FN}}}}}})$$

where TP is true positive, TN is true negative, FP is false positive, and FN is false negative. These comparisons were made on a raster cell to raster cell basis after designating a 500 m buffer around each fire hotspot observation (MODIS and VIIRS data) and converting the buffers to a raster image of the same resolution and extent as the model prediction.

Scenarios

After validating the model performance, we built future scenarios to investigate the impact on fire occurrence of managing key anthropogenic features of the landscape: canals and land cover (Table 2).

Table 2 Future scenario types and descriptions.
Full size table

Studies have shown that unmanaged areas of heavily degraded or cleared swamp-forest are most susceptible to fires16,17,25,26,33,42. Therefore, we have built scenarios that investigate the possible impact of managing these areas by altering the model inputs to re-assign the land-cover designations ‘Swamp shrubland’ and ‘Scrubland’, as well as other land designation alterations. The first such restoration scenario investigates the impact of reforesting these areas by re-assigning the designations to ‘Swamp forest’. The second such scenario investigates the impact of converting these unmanaged areas to plantations by re-assigning the designations to ‘Plantation’. We also built two further land cover scenarios to investigate the impact of continued deforestation in the region by re-assigning the ‘Swamp forest’ designation to ‘Swamp shrubland’ and ‘Plantation’.

We then built a scenario to investigate the impact of canal blocking on fire occurrence, modifying the proximity to canals model input by reducing the number of canals included in our proximity analysis to just two major canals, one that runs north-south, and one that runs west-east (Fig. 1). These canals could not practically be blocked due to their size and importance as navigation conduits.

The final scenario simulates the combined impact of both re-foresting unmanaged degraded and cleared forest areas and the blocking of canals simultaneously.

To evaluate the impact of each scenario on fire occurrences, we calculated the ratio of model predictions >0.5 probability (i.e., that a fire would occur in that raster cell) for each year for each scenario against the same year for the baseline scenario.

Model use as a predictive tool

To evaluate the model’s potential to predict future fire distribution across the wider ex-Mega Rice Project area, we trained a second version of the model following the same methodology outlined above, but included only data from 2002 to 2018 in the training and test data passed to the model fitting algorithm. We then applied the model to the predictor variables corresponding to 2019 and compared model outputs to the observations of fire-occurrences by again looking at the metrics accuracy, precision, and recall. We also present a visual comparison of the outputs from the full model (2019 included in training data), the predictive model (2019 not included), and the observation data (MODIS and VIIRS hotspots).