Notebook for the deforestprob Python module

• Notebook address: https://ecology.ghislainv.fr/deforestprob
• GitHub repository: https://github.com/ghislainv/deforestprob
• Author: Ghislain Vieilledent https://ecology.ghislainv.fr

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 use Madagascar as a case study.

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.

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 (1990-2000-2010)	BioSceneMada (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 roads	m	150
Places		distance to towns	m	150
Waterways		distance to river	m	150
Protected areas	Rebioma (4)	presence of protected area	-	30

1. http://bioscenemada.cirad.fr, deforestation maps derived from Harper et al. 2007 and 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 for Madagascar,
4. http://rebioma.net, SAPM ("Système des Aires Protégées à Madagascar"), 20/12/2010 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 download the derived variables from a repository on Zenodo (http://doi.org/10.5281/zenodo.259582). The bash script is available here.

In our example, fordefor2010.tif is a forest raster at 30m for the year 2010 considering the deforestation on the period 2000-2010 in Madagascar. 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. Pixels with the highest probability of deforestation are deforested first. The function computes a probability threshold above which pixels are deforested.

In our example, we consider an annual deforestation of roughly 100,000 ha for Madagascar. Considering the period 2010-2050, this would correspond to 4 Mha of deforestation. This number can of course be more precise and refined considering various 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 red 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).

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 only explains 8% of the deviance while the model including spatial random effects explains 29% 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.

External validation

We can also compute these indices using independant observations. To do so, we will compare the predicted deforestation on the period 2010-2014 with the true observed deforestation on the same period. The forest cover for the year 2014 for Madagascar can be downloaded from the BioSceneMada website.

We compute the observed forest-cover change between 2010 and 2014 using gdal_calc.py.

We then need to compute the number of hectares that have been deforested between 2010 and 2014. We use gdalinfo to do so.

400601 hectares have been deforested between 2010 and 2014.

Finally, we compute the accuracy indices comparing the predicted and observed forest-cover change map between 2010 and 2014 using function validation() from the deforestprob package.

References

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

Harper, G. J.; Steininger, M. K.; Tucker, C. J.; Juhn, D. & Hawkins, F. 2007. Fifty years of deforestation and forest fragmentation in Madagascar. Environmental Conservation. 34:325-333. doi: 10.1017/S0376892907004262

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