Ivory-Coast: modelling and forecasting deforestation

• Notebook address: https://github.com/ghislainv/workshopReCaREDD/blob/master/defor_Ivory-Coast/notebook.ipynb
• GitHub repository for deforestprob Python module: https://github.com/ghislainv/deforestprob
• Author: Ghislain Vieilledent https://ghislainv.github.io

This notebook will show you how to use the deforestprob Python module. The module can be used to estimate the spatial probability of deforestation in the tropics depending on various spatial environmental variables.

We first import the deforestprob module and the other necessary modules to run the script.

We also create a directory to hold the outputs with the help of the function .make_dir() from the deforestprob module.

Python preambule

Python is a widely-used programming language and is excellent for working with spatial data. For an introduction to Python, have a look at the Python documentation. To see how you can use Python for Geoprocessing, I advice reading the book Geoprocessing with Python by Chris Garrard.

If you want to see how to use a function in Python, simply write help("function_name") in the Python console.

Downloading and preparing the data

To model the spatial probability of deforestation, we need a map of the past deforestation and maps of potential explicative environmental variables. Environmental variables can be derived from topography (altitude, slope, and aspect), accessibility (distance to roads, towns, and forest edge), deforestation history (distance to previous deforestation) or landscape policy (location inside a protected area) for example. In our example, we use the following variables :

Product	Source	Variable derived	Unit	Resolution (m)
Deforestation maps (2000-2005-2010-2014)	Hansen et al. 2013 (1)	forest/non-forest	-	30
		distance to forest edge	m	30
		distance to previous deforestation	m	30
Digital Elevation Model	SRTM v4.1 CSI-CGIAR (2)	altitude	m	90
		slope	°	90
		aspect	°	90
Highways	OSM - Geofabrik (3)	distance to road	m	150
Places		distance to town	m	150
Waterways		distance to river	m	150
Protected areas	Protected Planet UNEP-WCMC (4)	presence of protected area	-	30

1. http://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.2.html, Hansen Global Forest Change v1.2 (2000-2014), Hansen et al. 2013
2. http://srtm.csi.cgiar.org, SRTM 90m Digital Elevation Database v4.1
3. http://www.geofabrik.de, data extracts from the OpenStreetMap project
4. http://www.protectedplanet.net, March 2007 version

To obtain the data, we can either (i) execute a bash script which downloads the raw data from the sources and computes the variables, or (ii) directly used the pre-computed variables that we copy to the data directory.

In our example, fcc05_10.tif is a forest-cover change raster at 30m for the period 2005-2010. We can plot this raster and zoom on a region with function .plot.forest() in the deforestprob package to see the deforestation data. The remaining forest appears in green and the deforested areas appear in red.

We can also plot the environmental raster files using function .plot.var() in the deforestprob package.

Sampling the observations

We use the function .sample() from the deforestprob module to sample 10,000 points (pixel centers) in the deforested areas and in the remaining forest (20,000 points in total). The input_forest_raster argument defines the path to the forest raster including the deforested pixels (with value=0) and the remaining forest pixels (value=1) after a given period of time. The random seed can be set with argument Seed to reproduce the data of the random sampling.

The .sample() function also extracts information from environmental variables for each sampled point. Sampling is done by block to allow the computation on large study areas (e.g. country or continental scale) with a fine spatial resolution (e.g. 30m). The var_dir argument indicates the directory including all the environmental raster files (they must be GeoTIFF files with extension .tif) that we want to test.

The .sample() function identifies the spatial cell for each sample point (sample point are grouped within a spatial cell). Spatial cells and grouped observations are used to estimate the spatial autocorrelation of the deforestation process. The csize argument define the width (in km) of the square spatial cells. Each spatial cell will be attributed a parameter. To avoid estimating too many parameters, width of the square spatial cells must be sufficiently large. Both points sampling and extraction of environmental data are done by block to avoid memory problems for big datasets.

A pandas DataFrame is returned by the function (each row being an observation) and the sampled observations are also written to a text file defined by the argument output_file.

Sampled observations can be plotted using function .plot.obs() from the deforestprob module. Dark red dots indicate deforestation observations and dark green dots indicate forest observations.

Descriptive statistics

Before modelling the deforestation, it is important to look at the relationship between environmental variables and deforestation. Using formulas from the patsy Python module, we can specify the relationships that we want to look at. In the example below, we plot the relationships between some continuous environmental variables and the probability of deforestation using function .plot.correlation() from the deforestprob package. Note that -1 must be set at the end of the formula.

The function .correlation() returns a serie of graphics that can be analyzed to choose the right relationship for each continuous variable (linear or polynomial for example).

In this example, we can see that a linear model should be sufficient to represent the relationship between the probability of deforestation and the standardized distance to the nearest road or town. On the contrary, it could be interesting to fit a polynomial model for the the standardized distance to previous deforestation (dist_defor variable) for which the relationship seems non-linear. Several models can be fitted and compared to see if a second-order or third-order polynomial relationship is justified.

Spatial deforestation model

We propose to use the Binomial iCAR model (Vieilledent et al. 2014) to estimate the deforestation probability of a pixel given a set of environmental variables. The Binomial iCAR model is a linear Binomial logistic regression model including an intrinsic Conditional Autoregressive (iCAR) process to account for the spatial autocorrelation of the observations. Parameter inference is done in a hierarchical Bayesian framework. The .model_binomial_iCAR() function from the deforestprob module includes a Metropolis-within-Gibbs algorithm written in pure C code to reduce computation time.

Figure caption: Parameter inference is done in a hierarchical Bayesian framework, which relies on the Bayes' theorem named after Reverend Thomas Bayes. Each parameter has a prior and an approximated posterior probability distribution from which we can compute the mean, standard deviation, credible intervals at 95%, etc.

For the deforestation process it is very important to take into account the spatial autocorrelation of the process with spatial random effects. Indeed, the selected fixed environmental variables are not able to fully explain the spatial variability of the deforestation process, especially when working at large geographical scales, such as the national or continental scale. Spatial random effects allow estimating a higher/lower probability of deforestation in a particular region (associated to unmeasurable or unknow factors) that is different from the mean probability of deforestation derived from the environmental factors included in the model. The Binomial iCAR model can be discribed as follow:

Ecological process

\begin{equation} y_i \sim \mathcal{B}inomial(\theta_i,t_i) \\ \text{logit}(\theta_i) = X_i \beta + \rho_{j(i)} \end{equation}

$y_i$: random variable for the deforestation process (0 if no deforestation, 1 if deforestation)

$\theta_i$: probability of deforestation

$t_i$: number of trials (always 1 in our example)

$X_i$: vector of values for environmental explicative variables

$\beta$: vector of fixed effect parameters

$\rho_j$: spatial random effect

$j(i)$: index of the spatial entity for observation $i$.

Spatial autocorrelation

The spatial autocorrelation is modelled with an intrinsic conditional autoregressive (iCAR) process:

\begin{equation} \rho_j \sim \mathcal{N}ormal(\mu_j,V_{\rho} / n_j) \end{equation}

$\mu_j$: mean of $\rho_{j'}$ in the neighborhood of $j$.

$V_{\rho}$: variance of the spatial random effects.

$n_j$: number of neighbors for spatial entity $j$.

Figure caption: Representation of the neighborhood for the intrinsic conditional autoregressive (iCAR) process. Target spatial cell $j$ has 8 neighbors in this case. Several observations (black points, equivalent to pixel centers in our case) can be located in each spatial cell. Deforestation probability in one spatial cell $j$ depends on deforestation probability in neighboring cells.

Before running the model, we add a column indicating the number of trials for each observation (1 in our case as we are considering a Bernoulli process). We then remove any observation with non-available data (NA) from the data-set. We also compute the number of neighbors (nneigh) and the neighbor identifiers (adj) for each spatial cell using function .cellneigh from the deforestprob module.

A model formula must also be defined to specify the explicative variables we want to include in the model. The formula allows specifying some variable transformations (such as standardization in our case). See the patsy module for more information. In our model, we included the following variables: location inside a protected area, altitude, distance to past deforestation (with a degree two polynomial), distance to forest edge, distance to nearest road and distance to nearest town. The formula must end with the name of the variable indicating the spatial cell for each observation (cell in our case).

Results of the model

Once the model has been fitted, we can print a summary of the model showing the parameter estimates. The 95% credible intervals obtained from the posterior distribution of each parameter, except distance to nearest town (dist_town), do not include zero, indicating that parameters are significantly different from zero. The variance of the spatial random effects (Vrho) is given together with the deviance value, which can be used to compare different statistical models (lower deviance is better). Looking at the parameter estimates, we can see that the deforestation probability is much lower inside protected areas and that deforestation probability decreases with altitude, slope, distance to past deforestation, forest edge, roads and towns. Parameter values are then coherent regarding the deforestation process and easy to interpret.

To check for the convergence of the Markov chain Monte Carlo (MCMC), we can plot the traces and the posterior distributions of the estimated parameters using method .plot() associated to the hSDM_binomial_iCAR class defined in the deforestprob module. This method returns the figures showing the traces and posterior distributions.

Predicting the spatial probability of deforestation

We use the model to predict the spatial probability of deforestation at the national scale for Madagascar. Before, doing so, we smooth the spatial random effects which have been estimated at a coarse resolution (10km in our example). To do so, we use the function .resample_rho() from the deforestprob module to resample the results at a finer resolution using a bilinear interpolation. The function writes a raster file to the disk with a resolution of the raster specified in the argument input_raster of the function (1km in our case). The function .resample_rho() returns a figure of the spatial random effects that can be plotted.

The .predict_hSDM() function of the deforestprob module can be used to predict the spatial probability of deforestation from an hSDM_binomial_iCAR model (i.e. an object of class hSDM_binomial_iCAR). The function writes a raster of predictions to the disk and returns a figure with the predictions that can be plotted. The prediction is done by block to avoid memory problems for big datasets. Functions will return NA for pixels with no forest or for pixels with missing environmental variables.

We can then plot the spatial probability of deforestation.

Predicting the future forest cover

Given the spatial probability of deforestation and the number of hectares that should be deforested, we can predict the future forest cover using function .deforest() from the deforestprob package. The number of hectares are converted into number of pixels to be deforested. The function computes a probability threshold above which pixels are deforested. Pixels with the highest probability of deforestation are deforested first.

Using GlobalForestWatch data and a tree cover threshold of 50%, we estimate an historical annual deforestation rate of roughly 568,270 ha for the Democratic Republic of Congo on the period 2001-2015. Considering 40 years on the period 2010-2050, this would correspond to 22.73 Mha of deforestation. This number can of course be more precise and refined considering various historical periods and deforestation scenarios (demographic growth, economic development, etc.).

Because deforestation probability is scaled on the interval [1, 65535] including 65535 values, there might be a difference between the deforested hectares set as input and the true deforested surfaces. The magnitude of the difference depends on the distribution of the deforestation probabilities. Usually, the distribution of the deforestation probabilities should be asymmetric with low frequencies for higher probabilities of deforestation.

As a consequence, the error should be negligible compared to the number of hectares set as input.

We can plot the predicted future forest cover in 2050. The purple areas represent the deforestation on the period 2010-2050. The green areas represent the remaining forest in 2050. Most of the remaining forest in 2050 are inside the protected areas or located in remote areas, at high altitudes and far from roads and big cities (for example in the Tsaratanana mountain region and around the Masoala peninsula, north-east Madagascar).

Carbon emissions

Carbon emissions associated to future deforestation can be computed using function .emissions() in the deforestprob package in combination with an aboveground forest biomass map. For example, we can use the pantropical 1km aboveground biomass map of Avitabile et al. (2016).

Model validation

Deviance explained

We can compute the percentage of deviance explained by the model with or without accounting for spatial autocorrelation. To do so, we need to compute the deviance of the full (or saturated) and null model. The full model assumes a perfect fit of the data and has a deviance equals to zero. The null model is a simple mean-only model ($\text{logit}(\theta_i)=\mu$). The deviance of the null model and of the model without spatial autocorrelation can be computed using a classical Binomial regression, without hierarchical structure. We will use the statsmodels Python module to do so.

We see that the model without spatial random effects explains 28% of the deviance while the model including spatial random effects explains 41% of the deviance. These percentage might appear rather low, but one should remember that we are trying to predict the deforestation at the pixel scale with a precision of 30 meters.

Predictions vs. observations

Another way of assessing model performance is to compute some accuracy indices based on a confusion matrix. We select the following accuracy indices: the Overall Accuracy (OA), the Specificity (Spe), the Sensitivity (Sen), Cohen's Kappa (K), the True Skill Statistics (TSS) and the Figure Of Merit (FOM). These indices are classical indices used to estimate the accuracy of deforestation models (see for example Dezécache et al. 2017, Vieilledent et al. 2013 and Pontius et al. 2011).

Internal validation

We first compute these indices on the sampling data-set. Here we are interested in seeing if our model is predicting well the location of the deforestation (relative deforestation probability), not the intensity of deforestation (absolute deforestation probability). We then transform the predicted probability of deforestation into binary values using the deforestation probability "observed" in the sampling data-set. This last value should be closed to 0.5 as we sampled as many forest and non-forest pixels, at the exception of the pixels with no-data value for environmental variables. The same proportion of deforested pixels will be applied to the predictions selecting the pixels with the highest predicted deforestation probabilities.

We use function accuracy_indices() from the deforestprob package to compute the accuracy indices.

The accuracy indices show that the model has good predictive abilities. Their value is relatively high (close to 1). Once we know the number of pixels to be deforested, the model correctly identifies the pixels that should be deforested first.

Conclusion

We showed here an approach to model and forecast tropical deforestation at the country level. The approach allows us to estimate the effect of each environmental factor on the probability of deforestation. Effects are easy to interpret: deforestation probability increases with accessibility and deforestation is much less important inside protected areas. We also showed that there is a high spatial variability in the deforestation process (see variance of the spatial random effects). Accounting for this additional variability is necessary to be able to forecast deforestation in a realistic way at the national scale.

Nonetheless, several improvements could be made to refine the modelling and forecasting approach:

• Predictions of the future forest-cover change should be validated with observations (for example using the 2010-2014 period).
• The choice of the historical forest-cover change map is really important. Several maps can be used (Hansen, Roadless or regional maps for example). Regarding the Hansen tree-cover map, different tree-cover thresholds can be tested.
• Other environmental explicative variables (such as soil types) could be added to the model.
• The statistical model can be refined removing some variables with non-significant effects or using higher order polynomials.
• Regarding the intensity of deforestation, several scenarios can be considered using various historical means on diferent periods, a deforestation intensity increasing with population growth or considering extreme scenarios (both in terms of mitigation or acceleration of the deforestation).

References

Avitabile, V., Herold, M., Heuvelink, G. B. M., Lewis, S. L., Phillips, O. L., Asner, G. P., Armston, J., Ashton, P. S., Banin, L., Bayol, N., Berry, N. J., Boeckx, P., de Jong, B. H. J., DeVries, B., Girardin, C. A. J., Kearsley, E., Lindsell, J. A., Lopez-Gonzalez, G., Lucas, R., Malhi, Y., Morel, A., Mitchard, E. T. A., Nagy, L., Qie, L., Quinones, M. J., Ryan, C. M., Ferry, S. J. W., Sunderland, T., Laurin, G. V., Gatti, R. C., Valentini, R., Verbeeck, H., Wijaya, A. and Willcock, S. 2016. An integrated pan-tropical biomass map using multiple reference datasets. Global Change Biology. 22: 1406–1420. doi: 10.1111/gcb.13139

Dezécache, C.; Salles, J.-M.; Vieilledent, G. & Hérault, B. 2017. Moving forward socio-economically focused models of deforestation. Global Change Biology. in press. doi: 10.1111/gcb.13611

Hansen, M. C.; Potapov, P. V.; Moore, R.; Hancher, M.; Turubanova, S. A.; Tyukavina, A.; Thau, D.; Stehman, S. V.; Goetz, S. J.; Loveland, T. R.; Kommareddy, A.; Egorov, A.; Chini, L.; Justice, C. O. & Townshend, J. R. G. 2013. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science. 342:850-853. doi: 10.1126/science.1244693

Pontius R. G. Jr., J. & Millones, M. 2011. Death to Kappa: birth of quantity disagreement and allocation disagreement for accuracy assessment. International Journal of Remote Sensing. 32:4407-4429. doi: 10.1080/01431161.2011.552923

Vieilledent, G.; Merow, C.; Guélat, J.; Latimer, A. M.; Kéry, M.; Gelfand, A. E.; Wilson, A. M.; Mortier, F. & Silander Jr., J. A. 2014. hSDM CRAN release v1.4 for hierarchical Bayesian species distribution models. Zenodo. doi: 10.5281/zenodo.48470

Vieilledent, G.; Grinand, C. & Vaudry, R. 2013. Forecasting deforestation and carbon emissions in tropical developing countries facing demographic expansion: a case study in Madagascar. Ecology and Evolution. 3:1702-1716. doi: 10.1002/ece3.550