Modeling historical deforestation on São Tomé

Marte Siebinga 01-07-2018

Supervisors: E.E. van Loon S.J. Norder

Modeling historical

deforestation on São Tomé

The application of a spatial explicit cellular automata model to

simulate historical land use dynamics

Abstract

São Tomé is a small African island located in the Gulf of Guinea. During the last four centuries the island’s environment has been drastically altered by deforestation. A broad overview of the deforestation trajectory of the island exists, however, it is lacking a spatial-explicit component. This study presents a simulation model of land use dynamics for simulating historical land use changes on São Tomé in a spatial explicit way. This was done in order to provide an insight into the most important drivers of historical deforestation. The spatial explicit model was made in Dinamica EGO, an environmental modelling platform based on CA. The model simulated land use dynamics on the island for two consecutive periods: 1891-1910 and 1910-1970. A total of five variables were included in the model: Humidity, Slope, Distance to Roads and Elevation and Distance to Deforestation. The model was calibrated and, ultimately, validated by using a Reciprocal Similarity Test. During calibration, the model achieved an acceptable spatial agreement of 87% for the first and 76% for the second period, both at minimum window size. For the independent validation of the simulated land use maps of 1947 and 1958, the map of 1958 achieved a higher similarity. The results show that Distance to Roads was the most important driver of historical deforestation on the island. Slope, Elevation and Humidity were also important drivers for historical land use change.

Table of Contents

3.4 Allocation of simulated land use changes...10

3.5 Calibration parameters...10 3.6 Validation...11 4. Results...14 4.1 Transition Rates...14 4.2 Weights of Evidence...14 4.3 Calibration...17 4.4 Output maps...20 4.5 Validation...21 5. Discussion...23 6. Conclusion...26 7. References...27 8. Appendices...30

Appendix A – Data table...30

Appendix B – Calculation transition rates...31

1. Introduction

During the last decades, tropical deforestation has been one of the biggest land use changes affecting the global environment (Geist & Lambin, 2002; Adhikari & Southworth, 2012). The decline of forest cover causes loss of biodiversity, soil degradation, influences the climate and increases the incidences of flooding. The disappearance of tropical forests is a result of several human activities, such as wood extraction, the extension of infrastructure and agricultural expansion. From those human activities, agricultural expansion is one of the leading land use change, associated with almost all deforestation cases (96%) (Geist & Lambin, 2002).

One African country whose environment is drastically altered by deforestation mainly resulting from agricultural expansion, is São Tomé e Príncipe. This former Portuguese colony is located in the Gulf of Guinea and consists of two small islands: São Tomé and Príncipe.Clearance of lowland forests on São Tomé began four centuries ago, when the Portuguese set up large plantations to cultivate sugar cane and cacao (Sayer et al., 1997; Eyzaguirre, 1986). During the peak of deforestation around 1910, the plantations claimed over 90% of the island’s surface and only the forests in the most inaccessible and wettest parts of the island remained undisturbed (Eyzaguirre, 1986).

A broad overview of the deforestation trajectory of the island exists, however, the information is very coarse and lacking a spatial explicit component. By making a spatial explicit land change model, changes in historical land use on the island can be described in a spatial explicit way. Spatial explicit land change models simulate the patterns of change in a landscape regarding coupled human-ecological dynamics. Therefore, spatial explicit land change models are essential for understanding human-environment interactions in the past and present (National Research Council, 2014; Petit & Lambin, 2002; Yang et al., 2014). The application of a spatial-explicit deforestation model on São Tomé might thus help to gain a better understanding of historical human impact on the ecosystems and environments of the island. Besides, modelling historical land use change on the island will provide an insight into the most important drivers of deforestation.

Dynamic spatial explicit models are mostly based on cellular automata (CA). In CA models space is represented as a grid of cells, each containing its own state. State change of a cell is dependent on a set of transition rules and on the states of the neighbouring cells. There are multiple advantages associated with CA models. First, dynamics can be modelled with a very high spatial resolution as CA are inherently spatial (White & Engelen, 1997). Besides, CA models

remote sensing or historical maps is relatively equal to the cellular data format, which makes processing straightforward (National Research Council, 2014). For those reasons, CA models have been widely used for modelling spatial dynamics in the past two decades (National Research Council, 2014).

This study aims to make a spatial explicit CA model of historical deforestation on São Tomé. The model will be made to reconstruct historical land use change with a spatial explicit component. Besides, the following research question will be answered: What were the most important drivers of land use change on São Tomé during the period 1891-1970? To model historical land use change, Dinamica EGO, a spatially explicit simulation model of landscape dynamics will be used (Soares-Filho et al., 2009).

2. Study Area

The island of São Tomé is located in the Gulf of Guinea just north of the equator (Figure 1). With an area of 857 km2_{, it is the main island of the country São Tomé e Príncipe. Due to its volcanic}

origin, the topography of the island is mountainous. The island has at least ten peaks of over 1000 m, of which the highest peak (Pico de São Tomé) is 2024 m (Sayer et al., 1997) The rugged terrain of the island makes a large area, mostly in the center and west, unsuitable for agriculture or forestry (de Lima, 2012). Furthermore, climate is influenced due to the rugged landscape. The mountains and the predominantly southerly winds create a rain shadow effect in the northeast and cover the southwest in a permanent cloud cover resulting in abundant rains throughout the year. Average rainfall in the southwestern parts of the island exceed 6000 mm annually while the average annual rainfall in the northeast stays below 600 mm (de Lima, 2012). The main dry period is between May and September. Besides, there can be a short dry period lasting a few weeks in January and February. The average annual temperature ranges from 23-25.5 °C at 300 m altitude to less than 13.5 °C at higher altitudes (de Lima, 2012). The soils of the volcanic islands are relatively deep, permeable and fertile, resulting in highly suitable conditions for agriculture (Eyzaguirre, 1986).

When the Portuguese discovered São Tomé around 1470, the island was most likely inhabited and fully covered by a dense rainforest known as obó. In the 1480s the cultivation of sugar cane began and the first period of deforestation started. Between 1550 and 1644 the island’s economy depended on sugar cane plantations. As a result most of the forests in the coastal lowland areas, particular in the drier northeast, were destroyed (de Lima, 2012). After a while, the Portuguese lost their commercial interests in the island and sugar cane production declined. Due to the decline of sugar cane plantations secondary forests (known as capoeira) regenerated the deforested areas in much of the areas above 400 m (Eyzaguirre, 1986).

After 1848, the Portuguese returned to the island and set up large plantations (roças) to cultivate the high-value export crops coffee and cocoa. This was the start of the second period of deforestation on the island. By the beginning of the 20th _{century the}

island became the world’s largest cocoa producer with an annual production reaching 36.148 metric ton (Eyzaguirre, 1989). In 1906, an approximate 70.000 ha of land was under cultivation by the roças. As a consequence of this intensive agricultural system deforestation peaked (Figure 2) (Eyzaguirre, 1986).

However, the period of the intensive cultivation of cocoa was only temporary. Several factors such as environmental degradation, increasing competitiveness of other colonies, and bad publicity regarding slavery led to the collapse of the roça system. Between 1910 and 1935, plantations were gradually abandoned and forest was regenerated. However, around 1935, the regeneration of secondary forest was interrupted by the world economic crisis and the accompanying reduced export values (Figure 2) (Eyzaguirre, 1986).

Figure 2: Changing agricultural practices over time on São Tomé (source: Eyzaguirre, 1986).

3. Methodology

In this research, a spatial explicit model of land use changes and GIS techniques were integrated in order to simulate historical land use dynamics in the study area. The study focusses on the second period of deforestation on the island, starting around 1850. To model the historical land use dynamics, the Dinamica software was used. This software is based on CA, available free of charge and has been extensively applied to model similar deforestation processes. Furthermore, Dinamica can model any number and type of transition over any time period. This time period can be divided into any number of time steps and phases with different, pre-defined transition rates. Despite, the software can still handle large amounts of data with a good performance (Soares-Filho et al., 2002). Therefore, Dinamica is considered an appropriate modelling software to simulate land use change. During the modelling process, different steps were carried out. All

these different steps are summarized in Figure 3. A more detailed description of the steps in the modelling process are given in the following sections.

Figure 3: Overview of the different steps carried out in the modelling process.

3.1

Input Data

For the model, three historical land cover maps of 1891, 1910 and 1970 were taken as input in order to determine the amount of change between different land use types. The land cover map of 1970 had already been digitized. The land cover maps of 1891 and 1910 still had to be digitized. This digitization was done manually by using ArcGIS 10.4.1.

Initially, the landscape map of 1970 presented four different land use types: Primary Forest, Secondary Forest, Shade Plantation and Non-forested Areas. However, the amount of land use types of the 1891 map was only limited to two. In order to prevent large data errors, the categorization of land use types of the different maps had to be the same. As a result, for this study, the land use types were reduced to Forest and Deforestation.

Besides the historical land use maps, a total of five landscape variables were used as input for the model to determine the location of change. Four of those variables were static, and thus kept constant over time, while one was dynamic and changed after each iteration. The static variables included in this study were: Humidity, Distance to Roads, Altitude and Slope. The dynamic variable used was Distance to Deforestation. As Forest cells may be converted into Deforestation cells during simulation, Distance to Deforestation had to be recalculated after each iteration. As a result Distance to Deforestation was considered a dynamic variable.

The selected variables were chosen based on their availability and data quality as well as on findings of preceding research and expert knowledge. Multiple studies showed that slope, altitude, climate, proximity to roads and proximity to previous deforested patches are important factors influencing tropical deforestation rates (Geist & Lambin, 2001; Kirby et al., 2006; Mertens & Lambin, 1997). Also, de Lima (personal communication, April 23, 2018) claims that the main predictors of deforestation on the island are accessibility and humidity. According to de Lima (personal communication, April 23, 2018), the most accessible and driest areas in the study area were deforested first, mainly because they were easy to burn down.

Dinamica only supports raster datasets in three formats: ER Mapper format, Geotiff and ArcView ASCII (Soares-Filho et al., 2009). Therefore, all the available data was prepared and exported in Geotiff format by using ArcGIS 10.4.1. Besides, all datasets had to be tied to the same coordinate system. In this study the WGS84 coordinate system was used and all input maps had the same spatial resolution of 133 m.

3.2

Transition Rates

To determine the amount of change between land use types in a certain simulation period, Dinamica calculates transition rates. Transition rates represent the total amount of changing cells for each transition in a given simulation period, without taking the spatial distribution of those changes into account (Soares-Filho et al., 2009).

For this study two different transitions were considered: Forest to Deforestation and Deforestation to Forest. As the transition rates varied through time (Figure 2), the model was split up into two parts and two different transition matrices were calculated. To calculate the transition rates for the first period, characterized by a high deforestation rate, the land use maps of 1891 and 1910 were used. The transition rates for the second period, characterized by the regeneration of forest, were calculated by using the land use maps of 1910 and 1970. However, as mentioned before, the regeneration of forest in this second period slowed down around 1935 (Figure 2). Therefore, the transition rates of this period were divided into a 60-40% ratio, with 60% of the transition occurring from 1910 to 1935, and 40% of the transitions occurring from 1936 to 1970.

3.3

Weights of Evidence

To identify the transition probability of each cell, Dinamica uses the weights of evidence method (Pérez-Vega et al., 2012; Stan et al., 2015). In this method, the probability of change is calculated by giving every variable a relative importance (weight of evidence). This weight of evidence is based on a Bayesian algorithm and represent each variable’s influence on the spatial probability of transition from one land use type to another (Stan et al., 2015). The weights of evidence for the transition Forest to Deforestation were determined by using the observed maps of 1891 and 1910. To determine the weights of the variables for the other transition, the observed maps of 1910 and 1970 were used.

Before the weight of evidence method was carried out, two pre-processing procedures were necessary. First, continuous variables, such as Altitude, Slope, Distance to Roads and Distance to Deforestation had to be categorized, since the weights of evidence method only applies to categorical data. This was done by using a function, incorporated in Dinamica, which calculates ranges according to a method adapted from Agterberg & Bonham-Carter (1990). For a more elaborate explanation of this method see Soares-Filho et al. (2009). Second, the correlation between pairs of variables was tested. This was done in order to assess the only assumption of the weight of evidence method that the variable maps have to be spatially independent. The correlation between pairs of variables was tested by using a function incorporated in Dinamica. This function performed pairwise tests for categorical maps, such as the Cramer test and the Joint-Uncertainty Information to assess the independence assumption (Soares-Filho et al., 2009). For this research, it was assumed that spatial dependence was low enough to maintain a particular variable, when values for Cramer’s Coefficient and the Joint Information Uncertainty were smaller than 0.5. This criterion was based on research conducted by Bonham-Carter (1994).

3.4

Allocation of simulated land use changes

Dinamica is a modelling software based on Cellular Automata (CA). CA models use discrete spatial units, arranged in a grid-like structure, as the basic units of simulation (National Research Council, 2014). Each cell is characterized by its state value, which may change after each iteration based on a set of transition rules or algorithm that is applied synchronously to all cells. This means that all cells change their state at the same time in accordance to the same set of rules. The rules consider the states of the cells in the neighborhood, which means changes in land use are dependent on the land use of the neighbouring cells (National Research Council, 2014; it et al., 2017; White & Engelen, 1997).

To model different scenarios, two CA functions or land use allocation algorithms were used by Dinamica: 1) the Expander and 2) the Patcher. The first function was dedicated only to the expansion or contraction of existing patches of a particular land use type. The second function accounted for the generation of new patches of a certain land use type (Soares-Filho et al., 2002; Mas et al., 2014). The combination of the two functions enabled the generation of numerous possibilities of spatial change patterns and was merged into the following equation (Soares-Filho et al., 2002):

Qij = r × (Expander Function) + s × (Patcher Function)

where Qij corresponds to the total amount of transitions of type ij, and r and s show the percentage of transitions executed by each function, with r+s = 1.

The values of the parameters associated with both functions were defined during the calibration of the model, which will be elaborated upon in the upcoming section.

3.5

Calibration parameters

Several parameters had to be calibrated before running the simulation. The first parameter was the percentage of transitions that will be executed by the Patcher and the Expander. Other parameters to be calibrated were mean patch size, patch size variance and patch isometry. These parameters can be independently altered for each of the two transition functions. By changing the parameters of the functions, the level of fragmentation or diversity in the landscape can be controlled (Mas et al., 2014). Increasing the mean patch size will result in the simulation of a less fragmented landscape. A more diverse landscape can be created by increasing the patch size variance. Finally, the patch isometry defines the aggregation of the simulated patches, with higher values corresponding to a more aggregated landscape (Soares-Filho et al., 2009).

To calibrate the parameters, the model was run several times to simulate the land use maps in 1910 and 1970. One parameter was changed at a time in order to get the best fit to the real spatial land use patterns on São Tomé. For every run, the simulated land use maps of 1910 and 1970 were compared with the observed land use maps of the same year. This was done by visual analysis as well as a Reciprocal Similarity Test (see next section for a more detailed description of this test). The parameters were adjusted until a maximum similarity between the simulated and observed land use maps was reached.

The values of mean patch size and patch size variance were calculated by using Dinamica and MATLAB. To calculate the parameters two maps of changes for the periods 1891-1910 and 1910-1970 were created. For each patch of change the area was calculated by using a function in Dinamica. Subsequently, the mean areas and the standard deviation of the patches were calculated in MATLAB. To calculate those values a distinction was made between areas generated by the Patcher function and areas created by the Expander function. This distinction was made visually and based on the maps of changes for the two periods and the locations of the changed patches. Since the model ran in annual steps, the parameters needed to be annual values. For instance, the mean

patch size of the Expander needed to be the annual mean patch size of all expanded areas. Therefore, the means and standard deviations were divided by the number of years for each of the two periods.

3.6

Validation

To assess the accuracy of the model, a Reciprocal Similarity Test was used. This test, which is already implemented in Dinamica, is a modification of the fuzzy similarity metrics established by Hagen (2003). The method developed by Hagen (2003) is based on the concept of fuzziness of location, which can be defined as the uncertainty and vagueness in the specification of the location of different categories of a map (Hagen, 2003). Each cell is assigned a fuzziness value depending on both the cell itself and the cells located in the neighbourhood of the central cell. The contribution of the neighbouring cells on the fuzziness is determined by an exponential decay function and decreases with increasing distance from the central cell. According to Hagen (2003), the choice of the size of the neighbourhood window depends on the vagueness of the data and the tolerance for spatial error. Besides an exponential decay function, a constant function can be applied to measure the accuracy of the model at various resolutions (Hagen, 2003).

By a cell-by-cell comparison of all the fuzziness values assigned to each cell, the similarity between two maps can be obtained. However, according to Hagen (2003), one-way comparison of a pair of maps tend to overestimate the spatial fit. Therefore, a two-way similarity comparison is applied as an alternative. To determine similarity, Soares-Filho et al. (2009) suggest to choose the minimum value from the two-way comparison, as random maps are likely to score higher. The similarity between the simulated and real maps can vary from 0-100%, where 100% indicates that both maps are identical (Soares-Filho et al., 2009).

However, in the abovementioned method, cells that did not change during simulation will cause an inherited similarity since simulated land use maps inherit the spatial patterns of initial maps (Almeida et al., 2008). To eliminate this inherited similarity, Dinamica uses a modified method which calculates similarity based only on the changed cells by using maps of changes (Soares-Filho et al., 2009). An overview of this method is shown in Figure 4.

For the validation of this model, two independent land use maps of 1947 and 1958 were used and compared to the observed land use maps by using a Reciprocal Similarity

Test with a constant decay function. Furthermore, multiple window sizes, varying from 1 x 1 to 25 x 25, were used. The simulated maps of 1910 and 1970 were validated as well, however, those validations are not considered independent since those specific land use maps were also used for calibration.

Figure 4: Overview of the Reciprocal Similarity Test used for validation of the model (source: Soares-Filho et al., 2009).

4. Results

4.1

Transition Rates

The average annual transition rates for the two different periods are presented in Table 1. As mentioned earlier, only two transitions were taken into account during this research. In the first period, the transition Forest to Deforestation was the only conversion observed. For the second period, both transitions occurred, however the transition from Forest to Deforestation was relatively small compared to the transition from Deforestation to Forest. For the second period, the annual transition rate was multiplied by the amount of years and divided in a 60-40 percent ratio. Then, the results were divided by the amount of years of the first (1910-1936) and second phase (1936-1970) of the second period. The results of the calculation are shown in Table 1, the full calculation can be found in Appendix B.

Table 1: Annual land use transition rates in the study area per time period.

Transition Annual rate (%)

1891-1910 Forest to Deforestation 6.30 1910 -1970 Forest to Deforestation 0.09 Deforestation to Forest 0.92 1910 - 1936 (60%) Forest to Deforestation 0.09 Deforestation to Forest 1.27 1936 - 1970 (40%) Forest to Deforestation 0.09 Deforestation to Forest 0.65

4.2

Weights of Evidence

Table 2 shows the results of the pairwise tests conducted to check for spatial dependence between the different pairs of variables. All the observed values of Cramer’s Coefficient and Joint Uncertainty were below 0.5. Therefore, based on the assumption of Bonham-Carter (1994), no

spatial dependence was observed between the used variables. Consequently, there was no need to exclude any of the variables from the model.

Table 2: Values of the Cramer’s Coefficient and Joint Uncertainty between the used variables.

The weights of evidence values (WoE) obtained for the different variables are shown in Figure 5. The higher and more positive the value, the greater is the influence of the particular variable on the transition considered. On the other hand, lower and more negative values depict a greater repelling effect on the particular transition. For the transition Forest to Deforestation, the highest value was found for the variable Distance to Roads. The lowest WoE value for this transition was found for the variable DEM. For the transition Deforestation to Forest, both the highest and the lowest WoE were found for the variable Distance to Roads. It is important to note that abnormal WoE values were automatically set to zero. For instance, for the variable Humidity, multiple weights were set to zero.

First variable Second variable Cramers coefficient Joint uncertainty

DEM Humidity 0.2500 0.1193

DEM Slope 0.0078 0.0009

DEM Dist. to Deforestation 0.2200 0.1561

DEM Dist. to Roads 0.1300 0.0806

Humidity Slope 0.0168 0.0006

Humidity Dist. to Deforestation 0.1739 0.0920

Humidity Dist. to Roads 0.3602 0.2200

Slope Dist. to Deforestation 0.0049 0.0005

Slope Dist. to Roads 0.0072 0.0006

Figure 5: WoE values obtained for each variable in relation to both transitions. The x axis show the upper limits of the assigned ranges (except for the categorical Humidity variable). Zero values are not shown in this

4.3

Calibration

According to the Reciprocal Similarity Tests and the visual analysis, the percentage of transitions from Forest to Deforestation executed by the Expander function was set to 90% for the 1891-1910 period. For the 1891-1910 – 1970 period, 60% of the transition from Forest to Deforestation and 100% of the transitions from Deforestation to Forest were executed by the Expander. Patch isometry in the first period was set to 1.5 for both functions. For the second period, patch isometry equaled 1.0 for all transitions and functions (Table 3). The parameters mean patch size and patch size variance were calculated by using various maps of changes. For the period 1891-1910 only one map of changes was used as simply one transition occurred (Figure 6). For the period 1910 – 1970, two maps of changes were utilized since both transitions were observed (Figure 6). The results of the calculations are shown in Table 3, for the entire calculation see Appendix 5.

Table 3: Values for the different parameters used in the model. The values indicated with a dash (-) did not have to be calibrated as the particular transition was only executed by the Expander function.

Parameter Patcher function Expander function

1891-1910 Fo re st to D ef or es ta ti

on Mean patch size (Ha) 857.37 857.37

Patch size variance (Ha) 1056.20 1056.20

Patch isometry 1.50 1.50

1910 -1970

Mean patch size (Ha) - 4.71

Patch size variance (Ha) - 1.23

Patch isometry - 1.00 D ef or es ta ti on to F or es

t 1910 - 1970 _{Mean patch size (Ha) 0.37 51.23}

Patch size variance (Ha) 0.65 125.91

Figure 6: Areas of land changed during the periods 1891 – 1910 (left) and 1910 -1970 (center and right). For the maps of change of 1910 – 1970, a distinction was made between the two different transitions and the areas changed by the Patcher and Expander function.

Figure 7 shows the results of the validation of the maps simulated during the calibration process. For the simulated map of 1910, the fuzzy similarity indices ranged from around 87% at minimum resolution to almost 95% at maximum resolution. For the 1970 map, a fuzzy similarity of 76% was reached for a window size of 1 x 1 and

Figure 7: Fuzzy similarity indices based on multiple window sizes for the simulated land use map in 1910 (left) and 1970 (right).

4.4

Output maps

Four of the simulated land use maps are shown in Figure 8. The initial land use map used for the simulation of the first period (1891 – 1910), was the observed land use map from 1891. For the second period, the observed land use map from 1910 was utilized as input map.

4.5

Validation

The simulated and observed land use maps used for the independent validation are shown in Figures 9 and 10.

Figure 9: The observed (left) and simulated (right) land use maps for 1947.

Figure 11 shows the fuzzy similarity indices per spatial resolution for the simulated maps in 1947 and 1958. For the simulated map in 1947, fuzzy similarity ranged from 35% for a window size of 1 x 1, to around 44% for a window size of 25 x 25. The fuzzy similarity indices for the 1958 map were higher and varied from 70% for minimum window size, to 84% for maximum window size.

Figure 11: Fuzzy similarity indices based on multiple window sizes for the simulated land use maps in 1947 (left) and 1958 (right).

5. Discussion

In general, the predictors of land use change found in this research are consistent with outcomes of previous research. For example, according to Nepstad et al. (2001) deforestation occurs close to roads as they provide physical access to forests. Also, Kaimowit & Angelsen (1998) and Kirby et al. (2006) claim that roads are the strongest driving factors of tropical deforestation. This complies with the results from this study as the WoE value for the variable Distance to Roads was the highest observed value and decreased with increasing distance.

Furthermore, an increase in slope and elevation results in lower and eventually more negative values of the WoE, as shown in the results (Figure 5). Since the variables slope and elevation essentially represent terrain accessibility, it can be concluded that the transition Forest to Deforestation was more likely to occur in accessible areas. Besides, the results show that dry areas favored the transition Forest to Deforestation as WoE values for Humidity are higher in dry areas. Both of these conclusions are in accordance to the findings that accessible and dry areas on São Tomé were deforested first (de Lima, personal communication, April 23, 2018).

According to Soares-Filho et al., (2001), deforestation processes generally radiate from previously deforested patches due to spatial autocorrelation. Based on this study, however, deforestation processes are not completely spatially auto-correlated. As shown in the results, WoE of Distance to Deforestation decreases with increasing distance. Nevertheless, this decrease only continues until a certain point, after this point WoE starts to increase again. It is likely that this increase is caused by a third variable that was not taken into account during this research. This third variable might have had a spatial structure in accordance with the particular distance, for instance patches within the same scale, influencing the results from the point WoE starts to increase again.

During the calibration of the first part of the model (1891-1910), similarities ranging from 87% to 95% were found (Figure 7). The fuzzy similarity indices of the second part of the model (1910-1970) ranged from 76% to 85% (Figure 7). The found similarity indices in this study are considered satisfactory based on results obtained in similar studies carried out with Dinamica. For example, Soares-Filho et al. (2002) achieved values ranging from 63.6% to 80%, and Stan et al. (2015) found values of 85% at minimum window size to values over 90% at maximum window size. The overall higher similarity of the first part compared to the second part of the model could be explained by the time span of both simulations. According to Khatibi et al. (2018), who

also used a Dinama EGO model to simulate land use change, the results of a simulation of 15 years were significantly more accurate and reliable compared with simulation results of the same model for a period of 29 years. As the time span of the first part of the model is 21 years and of the second part is 69 years, similarity might differ due to the length of the different periods.

For the independent validation of the 1958 map a spatial agreement of around 84% was achieved at maximum resolution (Figure 11). However, for the 1947 map a similarity of only 44% was reached at maximum resolution (Figure 11). Besides, it is notable that the plot of 1947 is flatter compared to the 1958 plot, and does not flatten out with increasing window size. According to Constanza (1989), the rapid increase of the plot in the beginning, such as for the 1958 plot, indicates that the patterns between the two compared maps are well matched, but the precise boundaries at maximum resolution are slightly off. On the other hand, if the plot is flatter, which is the case for the 1948 plot, the pattern is not well matched.

An explanation for the difference in similarity between the pairs of maps of 1947 and 1958 might be a wrong division of transition rates. As explained in the methodology, the transition rates of the second period were divided into a 60-40% ratio, with 60% of the transitions occurring during the period 1910-1936 and 40% after 1936. However, this ratio was solely based on the visual interpretation of the amount of transitions shown in Figure 2. As a result, transition rates during the first phase might have been too high causing lower similarity of the 1947 map. Other transition rates were used during calibration, however, this resulted in less accuracy of the 1970 map.

Another explanation for the lower similarity could be the inaccuracy of the observed 1947 land use map. According to Yang et al. (2014), maps, particularly historical maps, may contain errors, inaccuracies or might be completely wrong due to military or political intelligence reasons of the time period. Furthermore, historical maps may have geometrical irregularities compared to modern maps as a result of the purpose of the surveyor and the methods used at the time they were made (Yang et al., 2014). For those reasons, it is likely that the historical 1947 map, and the other historical maps used, contained some inaccuracies or small errors. Besides, by comparing the historical 1947 and 1958 maps, a lack of detail was observed for the 1947 map. This might have caused the lower similarity of the 1947 map compared to the 1958 map. Inaccuracies were also observed for the historical 1891 map, which was used as input. This map contained some geometrical irregularities compared to more modern maps of 1970 and

simulation results in general. Overall, the lack of accurate and detailed historical maps makes accurate spatial analysis of land use change more difficult and can thus be considered as limitation of this research.

Another consequence of the limited availability of accurate and detailed data is the simplification of the model. Initially, three land use types were modelled: Primary Forest, Secondary Forest and Deforestation. However, the simulation of those three land use types was only possible for the period 1958-1970 due to the lack of data. Consequently, land use types had to be reduced to Forest and Deforestation. Also, improving the model by additional variables was not possible since data was lacking. For instance, in comparable research by Yanai et al. (2012) and Mertens & Lambin (1997) variables such as distance to settlements and soil fertility were used, however, the addition of these variables was unable in this study.

An additional limitation of this study is the calculation of the parameters mean patch size and patch size variance. One drawback of the method used to calculate the parameters is that it is difficult to identify mean patch size and patch size variance for both the Patcher and Expander function separately. For the period 1910-1970, a distinction was made between the patches created by the Expander and the patches generated by the Patcher. It was assumed that the new patches adjacent to the main patch of forest in 1910 were created by the Expander. Patches not adjacent to this main patch of forest were created by the Patcher. However, this distinction was based merely on the visual observation of the map of change and can, thus, be considered arbitrary. For the period 1891-1910, it was not possible to make this distinction and it was assumed that the patches generated by both the Patcher and Expander had the same size and variance. The use of this method made the values for the parameters less precise, and might have influenced the accuracy of the simulation in general.

Additionally, due to the lack of spatial data on road development, this study assumed that existing roads did not change and had an equal influence overtime. Nevertheless, in reality, it is likely that infrastructure improves overtime, providing an improved accessibility and impacting deforestation processes.

6. Conclusion

To conclude, the outcomes generated by the model show land use change on São Tomé during the period 1891-1970 in a spatial explicit way. The model used was made in Dinamica EGO and based on CA. To model historical land use change the simulation was divided into two different periods: 1891-1910 and 1910-1970. Transition rates were calculated for both periods based on input maps from 1891, 1910 and 1970. The transition rates for the second period were split up to simulate the different rates of the regeneration of forest during that period. To determine the location of land use changes during simulation, several variables were given a relative importance by using the Weights of Evidence method. The variables used in this study were: Humidity, Slope, Elevation, Distance to Roads and Distance to Deforestation. Before the simulation, several parameters had to be calibrated. This was done by evaluating the outcomes visually as well as with a Reciprocal Similarity Test. The values of the parameters mean patch size and patch size variance were based on a calculation. After simulation, two simulated land use maps (1947 and 1958) were used for independent validation by comparing them to the observed land use maps of the same years. Based on the results of the validation, the model was considered acceptable.

The results of the Weights of Evidence method show that Distance to Roads was the strongest driving factor of historical deforestation on the island. During the period 1910 – 1970, areas close to the roads were more likely to be deforested than areas further away. Slope and elevation were important predictors as well. Land at high elevation or on steeper slopes was less likely to be deforested as accessibility of the land was lower. Humidity also influenced the deforestation processes on São Tomé. The results show that dry areas were more probable of being deforested. Overall, those findings are consistent with the results of previous, similar studies.

7. References

Adhikari, S., & Southworth, J. (2012). Simulating forest cover changes of Bannerghatta National Park based on a CA-Markov model: a remote sensing approach. Remote Sensing, 4(10), 3215-3243.

Agterberg, F. P., & Bonham-Carter, G. F. (1990, September). Deriving weights of evidence from geoscience contour maps for the prediction of discrete events. In Proceedings of the 22nd APCOM Symposium, Berlin, Germany (Vol. 2, pp. 381-395).

Almeida, C. D., Gleriani, J. M., Castejon, E. F., & Soares‐Filho, B. S. (2008). Using neural networks and cellular automata for modelling intra‐urban land‐use dynamics. International Journal of Geographical Information Science, 22(9), 943-963.

Bonham-Carter, G. F. (1994). Geographic information systems for geoscientists-modeling with GIS. Computer methods in the geoscientists, 13, 398.

Costanza, R. (1989). Model goodness of fit: a multiple resolution procedure. Ecological modelling, 47(3-4), 199-215.

de Lima, R. F. (2012). Land-use management and the conservation of endemic species in the island of São Tomé (Doctoral dissertation, PhD Dissertation. Lancaster University, Lancaster).

Eyzaguirre, P.B. (1986) The Ecology of Swidden Agriculture and Agrarian History in São Tomé. Cahiers d’Études africaines, 26, 113–129.

Eyzaguirre, P.B. (1989) The Independence of Sao Tome e Principe and Agrarian Reform. The Journal of Modern African Studies, 27, 671–678.

Geist, H. J., & Lambin, E. F. (2001). What drives tropical deforestation. LUCC Report series, 4, 116. Geist, H. J., & Lambin, E. F. (2002). Proximate causes and underlying driving forces of tropical deforestation: Tropical forests are disappearing as the result of many pressures, both local and regional, acting in various combinations in different geographical locations. BioScience, 52(2), 143-150.

Ghosh, P., Mukhopadhyay, A., Chanda, A., Mondal, P., Akhand, A., Mukherjee, S., ... & Hazra, S. (2017). Application of Cellular automata and Markov-chain model in geospatial environmental modeling-A review. Remote Sensing Applications: Society and Environment, 5, 64-77.

Hagen, A. (2003). Fuzzy set approach to assessing similarity of categorical maps. International Journal of Geographical Information Science, 17(3), 235-249.

Kaimowitz, D., & Angelsen, A. (1998). Economic models of tropical deforestation: a review. Cifor. Khatibi, A., Pourebrahim, S., & Danehkar, A. (2018). A cellular automata model for monitoring and simulating urban land use/cover changes toward sustainability. Journal of Environmental Engineering and Landscape Management, 26(1), 1-7.

Kirby, K. R., Laurance, W. F., Albernaz, A. K., Schroth, G., Fearnside, P. M., Bergen, S., ... & Da Costa, C. (2006). The future of deforestation in the Brazilian Amazon. Futures, 38(4), 432-453.

Mas, J. F., Kolb, M., Paegelow, M., Olmedo, M. T. C., & Houet, T. (2014). Inductive pattern-based land use/cover change models: A comparison of four software packages. Environmental Modelling & Software, 51, 94-111.

Mertens, B., & Lambin, E. F. (1997). Spatial modelling of deforestation in southern Cameroon: spatial disaggregation of diverse deforestation processes. Applied Geography, 17(2), 143-162. National Research Council. (2014). Advancing land change modeling: opportunities and research requirements. National Academies Press.

Nepstad, D., Carvalho, G., Barros, A. C., Alencar, A., Capobianco, J. P., Bishop, J., ... & Prins, E. (2001). Road paving, fire regime feedbacks, and the future of Amazon forests. Forest ecology and

management, 154(3), 395-407.

Pérez-Vega, A., Mas, J. F., & Ligmann-Zielinska, A. (2012). Comparing two approaches to land use/cover change modeling and their implications for the assessment of biodiversity loss in a deciduous tropical forest. Environmental Modelling & Software, 29(1), 11-23.

Petit, C. C., & Lambin, E. F. (2002). Long‐term land‐cover changes in the Belgian Ardennes (1775– 1929): model‐based reconstruction vs. historical maps. Global Change Biology, 8(7), 616-630. Sayer, J. A., Harcourt, C. S., & Collins, N. M. (1997). Conservation Atlas of Tropical Forests: Africa. Environmental Conservation, 24(1), 90.

Soares-Filho, B. S., Assunção, R. M., & Pantuzzo, A. E. (2001). Modeling the Spatial Transition Probabilities of Landscape Dynamics in an Amazonian Colonization Frontier: Transition probability maps indicate where changes may occur in the landscape, thus enabling better evaluation of the ecological consequences of landscape evolution. AIBS Bulletin, 51(12), 1059-1067.

Soares-Filho, B. S., Cerqueira, G. C., & Pennachin, C. L. (2002). DINAMICA—a stochastic cellular automata model designed to simulate the landscape dynamics in an Amazonian colonization frontier. Ecological modelling, 154(3), 217-235.

Soares-Filho, B. S., Rodrigues, H. O., & Costa, W. (2009). Modeling environmental dynamics with Dinamica EGO. Centro de Sensoriamento Remoto. Universidade Federal de Minas Gerais. Belo Horizonte, Minas Gerais, 115.

Stan, K., Sanchez-Azofeifa, A., Espírito-Santo, M., & Portillo-Quintero, C. (2015). Simulating deforestation in Minas Gerais, Brazil, under changing government policies and socioeconomic conditions. PloS one, 10(9), e0137911.

White, R., & Engelen, G. (1997). Cellular automata as the basis of integrated dynamic regional modelling. Environment and Planning B: Planning and design, 24(2), 235-246.

Yanai, A. M., Fearnside, P. M., de Alencastro Graça, P. M. L., & Nogueira, E. M. (2012). Avoided deforestation in Brazilian Amazonia: simulating the effect of the Juma Sustainable Development Reserve. Forest Ecology and Management, 282, 78-91.

Yang, X., Zheng, X. Q., & Chen, R. (2014). A land use change model: Integrating landscape pattern indexes and Markov-CA. Ecological Modelling, 283, 1-7.

8. Appendices

Appendix A – Data table

File/article name File format Description Source

N00E006 HGT Digital Elevation Model São Tomé de Lima * The Ecology of

Swidden Agriculture and Agrarian History in São Tomé

PDF Map of expansion of rainforest: 1910, 1947,

1982 Eyzaguirra (1986)

Landuse1970 Shapefile (.shp) Shapefile containing polygons of land use types on São Tomé in 1970. Four different land use types can be distinguished: Primary Forest, Secondary Forest, Shade Plantation and Non-forested Areas

de Lima *

Cartas Militares Folder with

JPEG-files Maps of land use and infrastructure on SãoTomé in 1958 de Lima * Picture05 JPEG Vegetation map of the island in 1957 de Lima * cartografia de s_tome3 JPEG Map of land use (cultivated area,

non-cultivated area) of the island in 1891 de Lima * ST_Roads Shapefiles(.shp) Shapefile containing locations of the roads,

on the island (Date unknown) de Lima * Picture01 JPEG Map showing the humidity on São Tomé and

Príncipe (1957) de Lima *

Appendix B – Calculation transition rates

Multi Step Transition Matrix showing the annual transition rates derived from Dinamica for the period 1910-1970:

From \ To | 1 2 1 | XXXX 0.0009253 2 | 0.0091933 XXXX

The transition rates for the period 1910-1970 were divided into a 60-40% ratio. The calculation is shown below.

Annual transition rate for transition Deforestation to Forest times the amount of years of the period:

0.00919 ×60=0.5514

60% of the transitions occurred during the period 1910 – 1936:

0.5514 × 0.6=0.33084

0.33084 ÷ 26=0.01272

40% of the transitions occurred during the period 1936 – 1970:

0.5514 × 0.4=0.22056

Appendix C – Calculations mean patch size and patch size variance

Table 1, 2 and 3 show the areas of all the new patches created during the first (1891-1910) and second period (1910-1970). The areas were calculated by using a function incorporated in Dinamica. Mean patch sizes and the standard deviations were calculated by using MATLAB. The MATLAB script is shown below.

Table 1: Areas of the different patches formed by the transition Forest to Deforestation during the period 1910 – 1970.

Patch ID Area_In_Cells Area_In_Hectares Area_In_Square_Meters

1 1 1.782385204 17823.85204 2 3 5.347155611 53471.55611 3 115 204.9742984 2049742.984 4 2 3.564770407 35647.70407 5 11 19.60623724 196062.3724 6 32 57.03632652 570363.2652 7 31 55.25394131 552539.4131 8 11 19.60623724 196062.3724 9 104 185.3680612 1853680.612 10 4 7.129540814 71295.40814 11 3 5.347155611 53471.55611

Table 2: Areas of the different patches formed by the transition Deforestation to Forest during the period 1910 – 1970.

Patch ID Area_In_Cells Area_In_Hectares Area_In_Square_Meters

1 159 283.3992474 2833992.474 2 12 21.38862244 213886.2244 3 3 5.347155611 53471.55611 4 82 146.1555867 1461555.867 5 5 8.911926018 89119.26018 6 8 14.25908163 142590.8163 7 5 8.911926018 89119.26018 8 2 3.564770407 35647.70407 9 11 19.60623724 196062.3724 10 39 69.51302294 695130.2294

11 12 21.38862244 213886.2244 12 761 1356.39514 13563951.4 13 44 78.42494896 784249.4896 14 6 10.69431122 106943.1122 15 7 12.47669643 124766.9643 16 1597 2846.46917 28464691.7 17 5 8.911926018 89119.26018 18 8 14.25908163 142590.8163 19 2 3.564770407 35647.70407 20 1 1.782385204 17823.85204 21 10 17.82385204 178238.5204 22 7 12.47669643 124766.9643 23 26 46.34201529 463420.1529 24 10 17.82385204 178238.5204 25 6 10.69431122 106943.1122 26 7 12.47669643 124766.9643 27 4 7.129540814 71295.40814 28 10 17.82385204 178238.5204 29 8 14.25908163 142590.8163 30 1 1.782385204 17823.85204 31 4 7.129540814 71295.40814 32 13 23.17100765 231710.0765 33 6 10.69431122 106943.1122 34 12 21.38862244 213886.2244 35 19 33.86531887 338653.1887 36 1 1.782385204 17823.85204 37 3 5.347155611 53471.55611 38 20 35.64770407 356477.0407 39 32 57.03632652 570363.2652 40 4 7.129540814 71295.40814 41 1 1.782385204 17823.85204 42 12935 23055.15261 230551526.1 43 7 12.47669643 124766.9643 44 16 28.51816326 285181.6326 45 11 19.60623724 196062.3724 46 4 7.129540814 71295.40814 47 10 17.82385204 178238.5204 48 1 1.782385204 17823.85204 49 13 23.17100765 231710.0765 50 1 1.782385204 17823.85204 51 1 1.782385204 17823.85204 52 13 23.17100765 231710.0765 53 24 42.77724489 427772.4489 54 5 8.911926018 89119.26018 55 3 5.347155611 53471.55611 56 2 3.564770407 35647.70407 57 2 3.564770407 35647.70407 58 5 8.911926018 89119.26018

59 36 64.16586733 641658.6733 60 4 7.129540814 71295.40814 61 50 89.11926018 891192.6018 62 5 8.911926018 89119.26018 63 2 3.564770407 35647.70407 64 6 10.69431122 106943.1122 65 50 89.11926018 891192.6018 66 13 23.17100765 231710.0765 67 2 3.564770407 35647.70407 68 11 19.60623724 196062.3724 69 2 3.564770407 35647.70407 70 155 276.2697066 2762697.066 71 4 7.129540814 71295.40814 72 4 7.129540814 71295.40814 73 86 153.2851275 1532851.275 74 3 5.347155611 53471.55611 75 5 8.911926018 89119.26018 76 9 16.04146683 160414.6683 77 5 8.911926018 89119.26018 78 6 10.69431122 106943.1122 79 1 1.782385204 17823.85204 80 9 16.04146683 160414.6683 81 5 8.911926018 89119.26018 82 1 1.782385204 17823.85204 83 3 5.347155611 53471.55611 84 1 1.782385204 17823.85204 85 1 1.782385204 17823.85204 86 2 3.564770407 35647.70407 87 1 1.782385204 17823.85204 88 44 78.42494896 784249.4896 89 2 3.564770407 35647.70407 90 2 3.564770407 35647.70407 91 17 30.30054846 303005.4846 92 6 10.69431122 106943.1122 93 9 16.04146683 160414.6683 94 2 3.564770407 35647.70407 95 13 23.17100765 231710.0765 96 3 5.347155611 53471.55611 97 8 14.25908163 142590.8163 98 3 5.347155611 53471.55611 99 1 1.782385204 17823.85204 100 4 7.129540814 71295.40814 101 58 103.3783418 1033783.418 102 2 3.564770407 35647.70407 103 1 1.782385204 17823.85204 104 7 12.47669643 124766.9643 105 103 183.585676 1835856.76

107 8 14.25908163 142590.8163

Table 3: Areas of the different patches formed by the transition Forest to Deforestation during the period 1910 – 1970. Patch ID Area_In_Cell s Area_In_Hectare s Area_In_Square_Meter s 1 1178 2099.64977 20996497.7 2 17101 30480.56937 304805693.7

% Script used to calculate the parameters mean patch size and patch size

% variance used for calibration of the land use change model clc

clear all close all

%% 1910 - 1970 transition 2 --> 1

areas = csvread('Areas_2_to_1_1910_1970.csv',1);

nr_yrs = 60; % number of years expander = [12,16,35,40,42,59,60,70,76]; % ids patches formed by expansion

areas_ha_exp = areas(expander,3);

mean_patch_size_exp = (sum(areas_ha_exp)/size(expander,2))/nr_yrs patch_size_var_exp = std(areas_ha_exp)/nr_yrs

patcher = [1:104];

patcher(expander) = []; % ids patches formed by patcher

areas_ha_pat = areas(patcher,3);

mean_patch_size_pat = (sum(areas_ha_pat)/size(patcher,2))/nr_yrs patch_size_var_pat = std(areas_ha_pat)/nr_yrs

%% 1910 - 1970 transition 1 --> 2

nr_yrs = 60; % number of years areas_ha_exp = areas(:,3);

mean_patch_size_exp = (sum(areas_ha_exp)/2)/nr_yrs patch_size_var_exp = std(areas_ha_exp)/nr_yrs

%% 1891 - 1910

areas = csvread('Areas_of_change_1891_1910.csv',1); nr_yrs = 19; % number of years areas_ha_exp = areas(:,3);

mean_patch_size_exp = (sum(areas_ha_exp)/2)/nr_yrs patch_size_var_exp = std(areas_ha_exp)/nr_yrs