Evaluating Dominant Land Use/Land Cover Changes and Predicting Future Scenario in a Rural Region Using a Memoryless Stochastic Method

Abstract

The present study used the official Portuguese land use/land cover (LULC) maps (Carta de Uso e Ocupação do Solo, COS) from 1995, 2007, 2010, 2015, and 2018 to quantify, visualize, and predict the spatiotemporal LULC transitions in the Beja district, a rural region in the southeast of Portugal, which is experiencing marked landscape changes. Here, we computed the conventional transition matrices for in-depth statistical analysis of the LULC changes that have occurred from 1995 to 2018, providing supplementary statistics regarding the vulnerability of inter-class transitions by focusing on the dominant signals of change. We also investigated how the LULC is going to move in the future (2040) based on matrices of current states using the Discrete-Time Markov Chain (DTMC) model. The results revealed that, between 1995 and 2018, about 28% of the Beja district landscape changed. Particularly, croplands remain the predominant LULC class in more than half of the Beja district (in 2018 about 64%). However, the behavior of the inter-class transitions was significantly different between periods, and explicitly revealed that arable land, pastures, and forest were the most dynamic LULC classes. Few dominant (systematic) signals of change during the 1995–2018 period were observed, highlighting the transition of arable land to permanent crops (5%) and to pastures (2.9%), and the transition of pastures to forest (3.5%) and to arable land (2.7%). Simulation results showed that about 25% of the territory is predicted to experience major LULC changes from arable land (−3.81%), permanent crops (+2.93%), and forests (+2.60%) by 2040.

1. Introduction

In many areas of the globe, the scale and extent of land use/land cover (LULC) changes have been affected by socioeconomic and biophysical factors (e.g., [1,2]). More specifically, environmental conditions, internal and external policies, the functioning of local and national markets, and the demographic situation have been pointed out as the main causes directly responsible for the LULC changes [3,4,5]. Therefore, over the last decade, for dynamic visualization and quantification of the spatial patterns of the LULC changes, Geographic Information Systems (GIS) have been widely used, since they allow for the use of vector or remote sensing data, expressed spatially and temporally [6,7,8]. Particularly, LULC change monitoring and assessment is a practice that assists a more conscious and thrifty management of this natural resource, which could result in effective environmental planning and informed policy decisions.

Commonly, in LULC change analysis, the LULC transitions are quantified and visualized based on a transition matrix by computing the inter-class percentage and by measuring, e.g., the persistence, gain/loss, net, and rate of the changes, in order to achieve comprehensive monitoring of physical changes over time [9,10]. Although these are the most commonly used measures, since they are very instructive for understanding the nuance of LULC changes, there are other measures that are equally informative and can contribute to more insights regarding LULC change processes [11] and to the design of anticipatory measure responses to changes. For instance, the quantification of both net change and swap (which comprises both the gain and loss of an LULC class in different locations of the study area), as well the gross gains/losses, allows for the detection of the inter-class transitions that have occurred as a systematic or as a random process of change [11,12]. Thus, the dominant signals of LULC changes can be distinguished.

Hence, whereas a systematic transition occurs by a continuous or general process of change, having a gradual tendency for progressive increment influenced by different factors (e.g., population growth or internal and external policies) [13], a random transition occurs as a single abrupt episode process of change influenced by factors that occur suddenly (e.g., floods or wildfires) [12,13]. Accordingly, distinguishing systematic (dominant) processes within a pattern of LULC changes, which tends to evolve in a progressive way driven by established and continued processes (see, e.g., [12,14]), allows one to focus on the mechanisms of change and further in-depth assessment of the potential change’s main drivers. Therefore, a better understanding of the LULC dynamics and better policies and measures in line can reduce the adverse impacts of dominant LULC changes [12,14,15].

In addition, apart from LULC change dynamics analysis, since the end of the 20th century, the prediction of future scenarios at different spatial and temporal scales has been extensively carried out [16,17,18]. In fact, the study of past and future LULC changes has been very promising for understanding past trends [2] and potential future developments, being a decision-making basis for environmental management and the sustainability of the ecosystems. Therefore, as a complex process, LULC dynamics is regularly modeled through several methods that are efficient in the simulation and prediction of future LULC [19], including the discrete-time Markov chain (DTMC) model (more commonly known as the Markov chain [20]) [21,22], the cellular automata (CA) model [18,23,24], and neural networks [25], among others. Of these approaches, the DTMC model is of countless importance for modeling LULC changes [26,27,28], as it has the advantages of a time dimension and is able to represent the LULC change data, since it has the ability to quantify not only the states of conversion between LULC classes but also the rate of conversion among them. Thus, it is a stochastic process, where space is discrete (i.e., mainly based on probabilities, not certainties).

Context and Background

The last decades of the 20th century and the first decades of the 21st century have shown that European and Mediterranean rural areas have been undergoing significant changes [29,30]. Indeed, one of the most remarkable changes in LULC in Europe, identified in Fuchs et al.’s [2] research, is the increase in forest to the detriment of the agricultural land. Likewise, the same trend has been observed in Portugal [30]. Indeed, the LULC changes in Portugal have been the most dynamic among the European Union (EU) countries [30,31].

For instance, despite the existence of numerous studies regarding the LULC changes in Portugal [32,33,34], the primary research focus has been on urban areas rather than rural areas. However, rural areas and agriculture have an important role in the regional and national economy of individual countries, as they constitute an important source of income for part of the population. Therefore, it is urgent to monitor the dynamics of rural LULC changes to promote and ensure a weighted and parsimonious use of the available resources.

Accordingly, given the stark changes of the Portuguese landscape, this study seeks to unravel the LULC changes situation in the Beja district, the largest Portuguese district in terms of area and the second in terms of percentage of agricultural land (about 47% in 2010); a rural region where we can find the Montado (named Montado in Portugal and Dehesa in Spain), an agrosilvopastoral agricultural heritage system [35,36], that has been indicated as a globally important agricultural system according to the Globally Important Agricultural Heritage System (GIAHS) program promoted by the Food and Agriculture Organization of the United Nations (FAO) [37].

Moreover, this research is an effort to evaluate the spatiotemporal LULC changes with a particular focus on the evaluation of the dominant LULC changes and to simulate the future changes, in order to determine the main LULC processes’ transformations and suggest new strategies for formulating better land use policies. The analysis was performed in a 10,229.05 km² area in the southeast of Portugal, which is experiencing marked landscape changes [32]. Thus, the main objectives of the present study were threefold: (i) identify and analyze the spatiotemporal dynamics of LULC over the past 23 years (1995–2018); (ii) quantify and analyze in depth the inter-class dominant (systematic) LULC transitions; and (iii) quantify the states of conversion between LULC classes and simulate and analyze the future LULC changes by 2040 using the DTMC model.

Our research is organized as follows: Section 2 details the study area, the used data, and the applied methodology. The results of the LULC changes and future LULC developments are described in Section 3. Finally, in Section 4, we discussed the main findings of our results and the implications of the applied methodology, presenting also the main conclusions of our research.

2. Materials and Methods

2.1. Study Area

The study area, the Beja district, is located between 07°37′58″ W and 07°41′17″ W longitude and between 38°08′35″ N and 38°10′08″ N latitude in the southeast of Portugal in the Nomenclature of Territorial Units for Statistics (NUTS) II Alentejo litoral and Baixo Alentejo (Figure 1). The district of Beja is bordered to the north by the Évora district, to the south by the Algarve region, to the east by Spain, and to the west by the district of Setúbal and the Atlantic Ocean. Beja is the largest Portuguese district, with an area of 10,229.05 km² and a resident population of 152,758 inhabitants in 2011 [38].

This region is characterized by a vast landscape of wheat, cork oaks, and olive trees, where the dominant land use is mixed agrosilvopastoral. The landscape of this region is distinctive, with both fragmented parcels and more compact parcels of land due to the different cropping calendars and field geometries. Dispersed settlements are the predominant urban form in this region. The climate is influenced by its distance from the coast, with a Mediterranean climate characterized by hot and dry summers and wet and cold winters. To the southeast, the territory is flat; to the north and west, the extensive plains are cut by tiny hills, the valley of the Guadiana River is the district’s main geographical feature, which crosses from north to south in the eastern part of the district. The Alqueva dam, the largest in the country, is located in this district.

2.2. Official Portuguese Land Cover (COS) Data

The official spatial land cover data (COS) produced by the Portuguese General Directorate for Territorial Development (DGT) for 1995, 2007, 2010, 2015, and 2018 were used. These maps have the highest cartographic detail for the study area and are available for free on the DGT website (http://mapas.dgterritorio.pt/geoportal/catalogo.html). The COS maps are produced by photointerpretation (i.e., over orthophoto maps), with vector (polygons) structure data, a cartographic accuracy of 0.5 m, a minimum mapping unit (MMU) of 1 ha, and a scale of 1:25,000 [39]. The five datasets have a hierarchical nomenclature system (a DGT report describes in detail the COS nomenclature characteristics [39]).

Spatiotemporal LULC change analysis was developed considering three major classes of the first hierarchical level of the COS nomenclature, namely, artificial surfaces, forest, and water bodies. In order to specify the major changes in croplands and characterize the type of agricultural activity in this region, we used the second hierarchical level of the COS nomenclature of the agriculture mega class, namely arable land, permanent crops, pastures, and heterogeneous areas. Therefore, a total of seven classes were used (

Table 1. Land use/Land cover (LULC) class description.

LULC Classes	Description
Artificial surfaces	Urban fabric; industrial, commercial, and transport units; mine, dump, and construction sites; artificial non-agricultural vegetated areas.
Agricultural-Arable land	Lands under a rotation system used for annually harvested plants and fallow lands, which are rain-fed or irrigated. Includes flooded crops such as rice fields and other inundated croplands.
Agricultural-Permanent crops	Lands not under a rotation system. Includes ligneous crops of standard cultures for fruit production such as extensive fruit orchards, olive groves, chestnut groves, walnut groves shrub orchards such as vineyards and some specific low-system orchard plantation, espaliers, and climbers.
Agricultural-Pastures	Lands that are permanently used (at least 5 years) for fodder production. Includes natural or sown herbaceous species, unimproved or lightly improved meadows and grazed or mechanically harvested meadows. Regular agriculture impact influences the natural development of natural herbaceous species composition.
Agricultural-Heterogeneous areas	Areas of annual crops associated with permanent crops on the same parcel, annual crops cultivated under forest trees, areas of annual crops, meadows and/or permanent crops which are juxtaposed, landscapes in which crops and pastures are intimately mixed with natural vegetation or natural areas.
Forest	Forests; shrub and/or herbaceous vegetation associations; open spaces with little or no vegetation.
Water bodies	Inland wetlands; coastal wetlands; marine waters.

). All LULC maps were clipped to the extent of the case study area, converted to a raster format of 100 × 100 m, and reclassified into seven classes. The cartographic information treatment and the spatiotemporal analysis were developed in QGIS software.

2.3. Spatiotemporal LULC Change Analysis

To identify and analyze the spatiotemporal dynamics of LULC in the Beja district, we first presented a brief evolution and distribution analysis of each class between 1995 and 2018. We quantified the changes for each year and reported them spatially and temporally. Secondly, we computed four change matrices in order to understand the main transition dynamics within each LULC class during the periods 1995–2007, 2007–2010, 2010–2015, and 2015–2018. The matrices were presented in graphical form. Each graph for each period shows the changes within each class from the initial time to the final time.

2.4. LULC Change Detection

For the LULC change detection in-depth analysis, we compared two maps (1995 and 2018) and produced a cross-tabulation matrix for the overall period (1995–2018). We decided to only analyze the overall period since each period comprises different years. While the first period has 12 years (1995–2007), the second (2007–2010) and the last (2015–2018) have 3 years, and the third has 5 years (2010–2015), so the change detection analysis could be influenced by the year’s discrepancy between periods. Thus, the transition matrix represents the percentage of LULC changes within each class (j), where the diagonal values indicate the percentage of LULC that persists $\left(P_{jj}\right)$ for each LULC class from initial time $\left(T_1\right)$ (1995) to the final time $\left(T_2\right)$ (2018), rows indicating the transitions from initial time $\left(P_{j+}\right)$, and columns indicating the transitions in the final time $\left(P_{+j}\right)$.

By following Pontius et al.’s [11] study, we obtained several measures from the transition matrix result and grouped them into two categories. The first category includes the (i) persistence, (ii) gain, (iii) loss, (iv) total change, (v) absolute value of net change, and (vi) swapping tendency. Since net change cannot detect the swap between each class, and therefore underestimates the LULC total change, the swapping tendency is calculated, since it comprises both the gain and loss of an LULC class in different locations of the study area [11]. Additionally, we presented the LULC vulnerability to transition by calculating the gain-to-loss, lost-to-persistence, and the gain-to-persistence ratios. A description of each measure is in

Table 2. Description of each measure for LULC change detection analysis.

Persistence $\left(P_{jj}\right)$	The amount of LULC for each class which remains from initial time and the later time.
Gain $\left(G_j\right)$	The amount of gains is the difference between the total value of each LULC class for the later time $\left(P_{+j}\right)$ and the value of persistency $\left(P_{jj}\right).$
Loss $\left(L_j\right)$	The amount of losses is the difference between the total for the initial time $\left(P_{j+}\right)$ and the value of persistency $\left(P_{jj}\right)$.
Total change $\left(C_j\right)$	The total change (overall change) is the sum of gain $\left(G_j\right)$ and loss $\left(L_j\right)$.
Net change $\left(D_j\right)$	The amount of net change is the maximum value of the gains and losses $\left(P_{+j}\right)$ minus the minimum value of gains and losses $\left(P_{j+}\right)$.
Swap tendency $\left(S_j\right)$	The amount of swap is two times the minimum value of the gains and losses $\left(G_j,L_j\right)$ for each LULC class.
Gain-to-Loss $\left(G/L\right)$	The tendency to experience more gain or loss is the division of gain $\left(G_j\right)$ and loss $\left(L_j\right)$. Values above 1 suggest larger tendency to gain than loss.
Loss-to-persistence $\left(L_p\right)$	The tendency to experience more loss or persistence is the division of loss $\left(L_j\right)$ and persistency $\left(P_{jj}\right).$ Values above 1 suggest larger tendency to loss than to persistence.
Gain-to-persistence $\left(G_p\right)$	The tendency to experience more gain or persistence is the division of gain $\left(G_j\right)$ and persistency $\left(P_{jj}\right).$ Values above 1 suggest larger tendency to gain than to persistence

.

The second category includes the measurement of systematic and random LULC changes by using the off-diagonal entries of the transition matrix (which indicates the transitions from one LULC class to the other) for the significant inter-class transitions detection. To detect if the inter-class transition is systematic or random, the expected LULC value was compared with the observed LULC value in the matrix. Following the statistical definitions outlined by Pontius et al. [40], a random transition indicates that the gain from other LULC classes was equivalent to the size of those other losing classes, or the loss to other LULC classes was equivalent to the size of those other gaining classes. Therefore, values close to zero or at zero indicate a partly systematic or random transition, respectively, while numbers far from zero indicate a systematic inter-class transition [11].

As such, the measurement of the systematic and random LULC transitions will result in three values in terms of gain $\left(G_{ij}\right)$ and in three values in terms of loss $\left(L_{ij}\right)$. Firstly, we compute the expected LULC gains percentage if the gain in each class were to occur randomly:

$$G_{ij}=\left(P_{+j}-P_{jj}\right)\left(\frac{P_{i+}}{{\sum }_{i=1, i\ne j }^j{P}_{i+}}\right)$$

(1)

where $\left(P_{+j}-P_{jj}\right)$ is the total gross gain of LULC class j, and $P_{i+}$ is the size of LULC class i in the initial time (1995). $\sum _{i=1, i\ne j }^j{P}_{i+}$ is the sum of the sizes of all LULC classes, excluding LULC classes of j, in the initial time. Secondly, we compute the expected LULC losses percentage if the loss in each class were to occur randomly:

$$L_{ij}=\left(P_{i+}-P_{ii}\right)\left(\frac{P_{+j}}{{\sum }_{j=1,j\ne j }^j{P}_{+j}}\right)$$

(2)

where $\left(P_{i+}-P_{ii}\right)$ is the observed total gross loss of LULC class i, and $P_{+j}$ is the size of LULC class j in the final time (2018). $\sum _{j=1,j\ne j }^j{P}_{+j}$ is the sum of the sizes of all LULC classes, excluding LULC classes of j, in the final time.

Third, we compute the differences between the observed LULC percentage and the expected LULC percentage under a random process of gain:

$$P_{ij}-G_{ij}$$

(3)

where $P_{ij}$ is the observed LULC, and $G_{ij}$ is the expected LULC if the gain in each class were to occur randomly. Fourth, we compute the differences between the observed LULC percentage and the expected LULC percentage under a random process of loss:

$$P_{ij}-L_{ij}$$

(4)

where $P_{ij}$ is the observed LULC, and $L_{ij}$ is the expected LULC if the loss in each class were to occur randomly.

2.5. The Discrete-Time Markov Chain (DTMC) Model

2.5.1. States Transition Probability of the Markov Chain

To quantify the states of conversion between the LULC classes and simulate and analyze the future LULC changes by 2040, we used the Discrete-Time Markov Chain (DTMC) model. Specifically, the DTMC model use matrices to model the LULC transition probabilities among a set of discrete states. Therefore, a DTMC is characterized by an array of random variables $X_1, X_2, \dots , X_n$ showing a memoryless property (i.e., a Markov property). The memoryless property postulates that the transition probabilities are independent of the chain state before its current state. Thus, the $X_{n+1}$ state distribution relies only upon the current state $X_n$ regardless of the $X_{n-1}, X_{n-2}, \dots , X_{n1}$ states.

$$P_r\left(X_{n+1}=x_{n+1}|X_1=x_1, X_2=x_2, \dots , X_n=x_n\right)=P_r\left(X_{n+1}=x_{n+1}|X_n=x_n\right)$$

(5)

Strictly, a Markov chain is a probabilistic automaton, where $X_n$ has a set of possible states $S=\left\{S_1, S_2, \dots , S_r\right\}$ entitled as the “state space of the chain.” The procedure begins in one of $S$ states and changes consecutively from one to another. Each change is named either step or transition, and the transition probability $p_{ij}$ is set as the probability of change from state $s_i$ to state $s_j$ in just one step.

$$p_{ij}=P_r\left(X_1=s_j|X_0=s_i\right)$$

(6)

whereas the probability to change from state $s_i$ to state $s_j$ in n steps is given by $p_{ij}^{\left(n\right)}=P_r\left(X_n=s_j|X_0=s_i\right)$.

If Equation (7) proves to be true, the DTMC is time-homogeneous, which indicates that the original transition probabilities remain unchanged through time.

$$P_r\left(X_{n+1}=s_j|X_n=s_i\right)=P_r\left(X_n=s_j|X_{n-1}=s_i\right)$$

(7)

For $k>0$, if $p_{ij}=P_r\left(X_{k+1}=s_j|X_k=s_i\right) \wedge p_{ij}^{\left(n\right)}=P_r\left(X_{n+k}=s_j|X_k=s_i\right)$, then the Markov chain is said to be time-homogeneous. The transition probability can be represented as a transition probabilities matrix (or simply transition matrix) $P={\left(p_{ij}\right)}_{i,j}$, where each element $\left(i,j\right)$ denotes the position of $p_{ij}$.

$$P=\left[\begin{array}{ccc}{p}_{ii}& \cdots & p_{ij}\\ ⋮& \ddots & ⋮\\ p_{ji}& \cdots & p_{jj}\end{array}\right]$$

(8)

The dispersion through the states can be transcribed as a stochastic $\left(\sum _i{x}_i=1,x_i\ge 0\right)$ row vector $x$. Hence, $x^{\left(1\right)}$ relates to $x^{\left(0\right)}$ as $x^{\left(1\right)}=x^{\left(0\right)}P$ and, recursively, we have $x^{\left(2\right)}=x^{\left(0\right)}{P}^2 \wedge x^{\left(n\right)}=x^{\left(0\right)}{P}^n$.

2.5.2. Classification States and Simulating Future States

The states can be classified by means of their periodicity. A state that is not absorbing is a transient one. A state $s_i$ is said to be transient if there is a state $s_k$ that is attainable from $s_j$, but j not being attainable from $s_k$. In other words, if $s_j$ is transient, there is no guarantee that, after leaving $s_j$, the system will return there. Suppose that, for a chain with initial state $x$, the number of periods required to return to that state is given by $T^{x\to x}$. If $P\left(T^{x\to x}<+\infty \right)=1\equiv P\left(T^{x\to x}=+\infty \right)>0$, the $x$ state is transient.

A state $s_j$ is said to be recurrent if it is not transient, i.e., if it is always possible to return to $s_j$. As such,

$$\begin{array}{ccc}\begin{array}{l}statexisrecurrentif\\ P\left(T^{x\to x}<+\infty \right)=1or\\ P\left(T^{x\to x}=+\infty \right)=1\end{array}& inadditionis& \{\begin{array}{l}nullrecurrentifP\left(T^{x\to x}\right)=+\infty \\ positiverecurrentifP\left(T^{x\to x}\right)<+\infty \\ absorbingifP\left(T^{x\to x}=1\right)=1\end{array}\end{array}$$

(9)

The time necessary to reach a specific state can also be computed. The so-called hitting time, i.e., the first passage $\left(s_i\to s_j\right)$, is given by the number of steps $\left(T_{ij}\right)$ that the chain requires until it hits the state $s_j$ for the first time, considering that $X_0=s_i$.

$$h_{ij}^{\left(n\right)}=P_r\left(T_{ij}=n\right)=P_r\left(X_n=s_j, X_{n-1}\ne s_j, \dots , X_1\ne s_j|X_0=s_i\right)$$

(10)

where $h_{ij}^{\left(n\right)}=\sum _{k\in S-\left\{s_j\right\}}{p}_{ik}{h}_{kj}^{\left(n-1\right)}$ and $h_{ij}^{\left(n\right)}=p_{ij}$.

A measure interconnected with $\left(h\right)$ h that is frequently used is the mean first passage time, i.e., the expected hitting time. In an ergodic Markov chain, the mean first passage time from $s_i$ to $s_j$ is the estimated number of steps necessary to go from the initial state $s_i$ to $s_j$ for the first time. This measure is the $\left(h\right)$ average value, explicitly ${\overline{h}}_{ij}=\sum _{n=1\dots \infty }n{h}_{ij}^{\left(n\right)}$. When delineating the first passage time, one can assume that $s_i=s_{ij}$ and thus obtain the first recurrence time $T_i=inf\left\{n\ge 1: X_n=s_i|X_0=s_i\right\}$.

Therefore, having the recurrent states, we calculated the mean recurrence time, i.e., the expected number of steps to return to the initial state:

$$r_i={\displaystyle \sum _{k=1}^{\infty }k\cdot P\left(T_i=k\right)}$$

(11)

Finally, as the period under analysis presents different cycles, that is more dynamic cycles (1995–2007 and 2010–2015) interspersed with less dynamic ones (2007–2010 and 2015–2018), we chose to perform future simulations based on trends of the whole set (1995–2018). Thus, based on the year 2018 and using the trends recorded over 23 years, we created a scenario of LULC for 2040 (instead of 2041 for rounding purposes).

2.5.3. Uncertainty of Future States

Uncertainty can then be characterized by the probability vectors that are produced as a sub-product of simulation, made available by most stochastic modeling procedures [41]. These probabilities provide useful information about the quality of the resulting classification in terms of the uncertainties involved. To fully explore the information of the probability vector, additional measures of uncertainty are needed [42]. The greater the uncertainty, the lower the probability that a pixel will be associated with a given class, so uncertainty is the complement of the probability. Formally, if we considerer $P_{ij}$, the probability that a given pixel $j$ is associated with class $i$, the uncertainty can be described by $U_{ij}=1-P_{ij}$.

However, normally the true class of the pixel is not known and, consequently, the amount of information required to indicate the pixel class is also unknown. The pixel entropy is therefore defined as the expected information index of a piece of information that reveals your true class. For this purpose, the entropy measure combines the uncertainties of the various pixel classes, weighing them by their probabilities. In this way, the global entropy can be estimated from the normalized values of probabilities by class, using the expression $E_{ij}=-{\sum }^{\text{}}{P}_{ij}\ast \log {}_2{}^{}P{}_{ij}$.

3. Results

3.1. LULC Spatiotemporal Evolution and Distribution

The spatiotemporal analysis results showed that the Beja district has undergone major changes between 1995 and 2018 (

Table 3. LULC from 1995 to 2018 (in %).

LULC Class	1995	2007	2010	2015	2018
Artificial surfaces	0.79	1.04	1.09	1.21	1.24
Agricultural areas	69.43	63.30	62.88	63.12	63.96
Arable land	26.64	23.69	22.83	22.77	20.70
Permanent crops	7.43	8.80	9.87	10.25	12.56
Pastures	14.13	10.63	10.11	10.68	11.50
Heterogeneous areas	21.23	20.18	20.07	19.42	19.20
Forest	28.35	33.68	33.94	33.53	32.66
Water bodies	1.44	1.99	2.09	2.14	2.14

). At first sight, the predominance of agricultural and forestry areas over the entire study period was observed. However, the occupancy tendency of both classes showed large fluctuations. Examining the results, between 1995 and 2010, there was a loss of almost 6.6% of the agricultural area, but this mega class continued to occupy the largest portion of the territory (in 2010 about 62.88%). Subsequent years (2015–2018) saw an increase in agricultural area (+1.08%).

Indeed, except for the permanent crops, there has been a significant decrease in remaining croplands. Particularly, arable land, the cropland with the most significant predominance in the Beja district (in 2018 about 20.7%), presented a decrease of almost 6% between 1995 and 2018, while heterogeneous areas and pastures decreased by about 2% and 2.6%, respectively. Heterogeneous areas are mainly concentrated near forestland, and the pastures class almost disappears in the central region of the Beja district. Permanent crops increased by about 5%, particularly in the northern part of the district, where the patch correspondent to this class is noticeable (Figure 2).

By contrast, the forestland presented an area increase of almost 5.6% between 1995 and 2010 (in 2010 about 34%), but in subsequent years (2010–2018) there was a loss of almost 1.3% (in 2018 about 32.7%). It is noticeable that the green area corresponds to the forestland class in the SSW of the study area, probably due to its being closer to the Atlantic Ocean and to its orography (plateaus). However, the increase in forest was especially near the wetlands, untangling the southeast part of the district. As for the artificial surfaces, although substantially reduced in this region, they present a linear growth increment between the period under analysis (close to 0.5%), especially in the north of the district near the Beja airport. The same occurred with the wetlands class with an overall increase of 0.7%, with some visible expression in the northeast part of the district (Figure 2).

3.2. Main LULC Transitions by Period

The analysis of the inter-class transitions by period is present using the cross-tabulation matrices in graphical form for the four periods under analysis (Figure 3). Specifically, during the first period (1995–2007), almost 20% of the LULC changed. The results revealed a clear tendency for each type of agricultural use to transition to another agricultural use or to transition to forest. Additionally, the forest class and the pastures class were the two that presented the highest LULC transition values. Indeed, the first period experienced a significant decrease in croplands (−6.3%), mostly due to forestland, with transitions from the pastures (3.8%), heterogeneous areas (1.8%), and arable land (1.7%) classes. However, about 1.8% of forestland also transitioned to agricultural areas. By contrast, arable land and pastures exchange between themselves.

As for the second period (2007–2010), about 3% of the LULC changed. This period experienced the same trend: a decrease in croplands, albeit much less significant (−0.42%). The arable land, pastures, and permanent crops classes transitioned between each other, the highest area transitioning from arable land to permanent crops (0.8%).

The subsequent period (2010–2015) experienced about 8.5% of LULC changes. However, this period counters the decreasing trend in cropland, albeit modestly, as croplands increased but only by about 0.24%. Thus, a transition tendency similar to the first period can be observed, since the major recorded transitions were between agricultural areas and forestland. Accordingly, the heterogeneous areas and forestland classes transitioned between each other, with heterogeneous areas transitioning to forest (1.18%) and forest transitioning to heterogeneous areas (1.55%).

During the last period (2015–2018), almost 4.7% of the LULC changed. This period experienced the highest gain in agricultural areas (+0.84%). Once again, arable land transitioned mainly to pastures (0.7%), and pastures transitioned mainly to arable land (0.5%), while heterogeneous areas transitioned to arable land (0.7%) and to pastures (0.7%). Indeed, the most substantial transition recorded was from arable land to permanent crops (1.8%), followed by the transition from forest to pastures (0.7%), and from arable land to pastures (0.6%). Unlike the past two periods, we verified for the first time more losses than gains for the forest class.

3.3. LULC Change Detection

The in-depth analysis of the overall LULC change detection shows that, for the overall period (1995–2018), about 72.1% of the LULC persisted, while 27.9% exhibited transitions. Briefly, the highest LULC transitions were experienced from arable land to permanent crops (4.96%), followed by the transition from pastures to forest (3.48%) (

Table 4. LULC transition matrix 1995–2018 (in %) *.

		2018
	LULC Class	Artificial Surfaces	Arable Land	Permanent Crops	Pastures	Heterog. Areas	Forest	Water Bodies
1995	Artificial surfaces	0.76	0.00	0.01	0.01	0.00	0.00	0.00
	Arable land	0.16	16.47	4.96	2.92	0.17	1.78	0.18
	Permanent crops	0.03	0.53	6.14	0.12	0.42	0.17	0.03
	Pastures	0.13	2.70	0.92	6.52	0.26	3.48	0.12
	Heterog. areas	0.10	0.64	0.40	1.18	16.23	2.53	0.16
	Forest	0.05	0.37	0.11	0.76	2.13	24.66	0.27
	Water bodies	0.01	0.00	0.01	0.01	0.00	0.04	1.36

* The highest LULC transitions are highlighted in bold.

). In addition, the overall period transitions highlighted that forest experienced high gains from agricultural classes, while both artificial surfaces and water bodies showed a minor transition to other LULC classes.

From

Table 5. Summary of LULC changes (1995–2018) (in %).

LULC Class	${\mathit{P}}_{\mathit{j}}$	${\mathit{G}}_{\mathit{j}}$	${\mathit{L}}_{\mathit{j}}$	${\mathit{C}}_{\mathit{j}}$	${\mathit{S}}_{\mathit{j}}$	${\mathit{D}}_{\mathit{j}}$	G/L	${\mathit{L}}_{\mathit{p}}$	${\mathit{G}}_{\mathit{p}}$
Artificial surfaces	0.76	0.48	0.03	0.50	0.05	0.45	18.32	0.03	0.63
Arable land	16.47	4.23	10.16	14.39	8.46	5.93	0.42	0.62	0.26
Permanent crops	6.14	6.42	1.30	7.71	2.59	5.12	4.96	0.21	1.05
Pastures	6.52	4.99	7.61	12.60	9.97	2.62	0.66	1.17	0.77
Heterogeneous areas	16.23	2.97	5.00	7.97	5.94	2.03	0.59	0.31	0.18
Forest	24.66	8.00	3.69	11.69	7.38	4.31	2.17	0.15	0.32
Water bodies	1.36	0.78	0.07	0.85	0.15	0.70	10.47	0.05	0.57
Total	72.14	27.86	27.86	27.86	17.27	10.59

Note: $P_j$ = persistence, $G_j$ = gain, $L_j$ = loss, $C_j$ = total change, $S_j$ = Swap, $D_j$ = net change, G/L = gain-to-loss ratio, $L_p$ = loss-to-persistence ratio, $G_p$ = gain-to-persistence ratio.

, it can be verified that, during the overall period, both the artificial surfaces and water bodies classes experienced a lower percentage of persistence (0.76% and 1.36%, respectively). These results are mostly due to their lower percentage of territory occupancy within the period analysis. Thus, artificial surfaces occupied about 1.2% of the district area, meaning almost 0.5% of gains and only 0.03% of losses. Therefore, the changes in the artificial surfaces class is mostly net change (about 0.45%). Similarly, the water bodies class experienced almost 0.8% of gains against 0.1% of losses between 1995 and 2018. Thus, most of the change in the water bodies class is a net change (0.7%). Hence, artificial surfaces and water bodies experienced a gain–loss ratio of 18.3 and 10.5, respectively; however, according to the loss-to-persistence and gain-to-persistence ratios, both classes presented a greater tendency to persist than to gain or lose. Figure 4 presents the net changes for both the artificial surfaces and water bodies classes. While artificial surfaces gains are highly dispersed across the district and usually founded near spots with the same LULC, water bodies gains are not clustered predominantly near to each other. Due to the minimal losses during 1995–2018, it is difficult to identify the dispersion of losses across the district.

In contrast, forest, arable land, and heterogeneous areas presented a high percentage of persistence (24.66%, 16.47%, and 16.23%, respectively), while pastures and permanent crops presented a relatively low percentage of persistence (6.52% and 6.14%, respectively). These results are mostly due to each LULC class percentage of territory occupancy during all the years under study. Nevertheless, permanent crops and forest classes presented more gains than losses during the overall period. The gain–loss ratio indicates that permanent crops presented 4.9 times more gain than loss, while forest presented a ratio of 2.2 times more. However, the influence of net changes was found to be the main driver in permanent crops (66% of the total change for permanent crops), while swapping was the main driver in the forest class (63% of the total change for forest).

However, according to the loss-to-persistence and gain-to-persistence ratios, the forest class experienced a greater tendency to persist than to gain or lose, while the permanent crops class presented a greater tendency to persist than to lose and a greater tendency to gain than to persist. Figure 5 presents the net changes for both the permanent crops and forest classes. The gains of forestland are dispersedly distributed across the Beja district and are usually found near spots with the same LULC. Comparatively, permanent crops gains are primarily distributed in the north and east. Gain spots are grouped mainly near to each other, while losses are minimal and show no distribution pattern.

The remaining agricultural classes presented more losses than gains, with a loss–gain ratio of 2.4, 1.5, and 1.7 for arable land, pastures, and heterogeneous areas, respectively. These three agricultural classes indicate both swap and net changes between 1995 and 2018. However, the influence of swapping was found to be the main driver of the total change, which was a change by 59%, 79%, and 75% for arable land, pastures, and heterogeneous areas, respectively. While both arable land and heterogeneous areas presented a greater tendency to persist than to gain or lose, the pastures class presented a greater tendency to lose than to persist and a greater tendency to persist than to gain.

Figure 6 presents the net changes for arable land, pastures, and heterogeneous areas. The three classes present a very distinct change pattern. Arable land losses are mainly consolidated in the north and east of Beja, while the gains are dispersed in the south of the district, predominantly near spots with the same LULC. As for pastures, the gains are spread over the entire territory, while the losses can be found in the south and close to the east. Lastly, the heterogeneous areas show large gains in the center of the territory, while losses are distributed mainly close to where there is forestland.

3.4. Detection of Random and Systematic LULC Transitions

The random and systematic LULC transitions analysis showed that not all the inter-class transitions during the period under analysis were systematic (dominant) processes within a pattern of LULC changes. Based on

Table 6. Inter-class LULC transitions gains.

	2018
1995	Artificial Surfaces	Arable Land	Permanent Crops	Pastures	Heterog. Areas	Forest	Water Bodies
	(i) Expected LULC if the process of change in each class were to occur randomly (%)
Artificial surfaces	0.76	0.05	0.05	0.05	0.03	0.09	0.01
Arable land	0.13	16.47	1.85	1.55	1.00	2.98	0.21
Permanent crops	0.04	0.43	6.14	0.43	0.28	0.83	0.06
Pastures	0.07	0.81	0.98	6.52	0.53	1.58	0.11
Heterog. areas	0.10	1.22	1.47	1.23	16.23	2.37	0.17
Forest	0.14	1.63	1.97	1.65	1.07	24.66	0.22
Water bodies	0.01	0.08	0.10	0.08	0.05	0.16	1.36
	(ii) Observed LULC minus the proportion expected under the random process (%) *
Artificial surfaces	0.00	−0.04	−0.05	−0.04	−0.03	−0.09	0.00
Arable land	0.03	0.00	3.12	1.37	−0.84	−1.20	−0.03
Permanent crops	−0.01	0.10	0.00	−0.31	0.14	−0.66	−0.03
Pastures	0.06	1.88	−0.06	0.00	−0.28	1.91	0.01
Heterog. areas	−0.01	−0.59	−1.07	−0.06	0.00	0.16	−0.01
Forest	−0.08	−1.27	−1.85	−0.89	1.06	0.00	0.05
Water bodies	0.00	−0.08	−0.09	−0.08	−0.05	−0.12	0.00

* The magnitude of the value indicates the percent of LULC. The most systematic transitions are highlighted in bold.

and

Table 7. Inter-class LULC transitions losses.

	2018
1995	Artificial Surfaces	Arable Land	Permanent Crops	Pastures	Heterog. Areas	Forest	Water Bodies
	(i) Expected LULC if the process of change in each class were to occur randomly (%)
Artificial surfaces	0.76	0.01	0.00	0.00	0.01	0.01	0.00
Arable land	0.16	16.47	1.61	1.47	2.46	4.19	0.27
Permanent crops	0.02	0.31	6.14	0.17	0.28	0.48	0.03
Pastures	0.11	1.78	1.08	6.52	1.65	2.81	0.18
Heterog. areas	0.08	1.28	0.78	0.71	16.23	2.02	0.13
Forest	0.07	1.13	0.69	0.63	1.05	24.66	0.12
Water bodies	0.00	0.02	0.01	0.01	0.01	0.02	1.36
	(ii) Observed LULC minus the proportion expected under the random process (%) *
Artificial surfaces	0.00	0.00	0.00	0.01	0.00	−0.01	0.00
Arable land	0.00	0.00	3.35	1.44	−2.29	−2.41	−0.09
Permanent crops	0.01	0.22	0.00	−0.05	0.13	−0.31	0.00
Pastures	0.02	0.92	−0.16	0.00	−1.39	0.68	−0.06
Heterog. areas	0.02	−0.64	−0.37	0.46	0.00	0.51	0.03
Forest	−0.02	−0.77	−0.57	0.13	1.08	0.00	0.16
Water bodies	0.01	−0.01	0.00	0.00	−0.01	0.02	0.00

* The magnitude of the value indicates the percent of LULC. The most systematic transitions are highlighted in bold.

, we are able to determine what the dominant transitions in each LULC class would be if the gain or loss in each class were to occur randomly. Therefore, the arable land and permanent crops classes experienced the strongest dominant signal of change. While arable land systematically transitioned to permanent crops, permanent crops were systematically gaining from arable land simultaneously (3.12%), and arable land was also systematically losing to permanent crops (3.35%). However, permanent crops systematically avoided losses to arable land.

Likewise, while arable land systematically transitioned to pastures, the pastures class was systematically gaining from arable land simultaneously (1.37%), and arable land was also systematically losing to pastures (1.44%). By contrast, while pastures were systematically transitioning to arable land, arable land was systematically gaining from pastures simultaneously (1.88%), and pastures were also systematically losing to arable land (0.92%). Thus, the arable land and pastures classes systematically transitioned between each other. Analyzing the situation more closely, based on Table 6, it should also be noted that permanent crops and arable land systematically avoided gaining from forest (−1.85% and −1.27%, respectively), and permanent crops or arable land avoided gaining from heterogeneous areas (−1.07% and −0.59%, respectively). In Table 7, it can be seen that arable land systematically avoided losing to heterogeneous areas (−2.29%) or to forest (−2.41%), and pastures avoided losing to heterogeneous areas (−1.39%).

Even though forest systematically gained from pastures (1.91%), there was a weak signal of pastures losing to forest (0.68). The analysis of the results in Table 6 also recognized that pastures systematically avoided gaining from forest (−0.89%), while Table 7 indicated that pastures avoided losing to heterogeneous areas (−1.39%). On the one hand, heterogeneous areas systematically gained from forests (1.06%), while heterogeneous areas systematically avoid losing to forest (0.51%). However, forest was systematically losing to heterogeneous areas (1.08%). From Table 6, it also appeared that forest and heterogeneous area classes systematically avoided gaining from arable land (−1.20% and −0.84%, respectively), and forest avoided gaining from permanent crops (−0.66%).

3.5. States Transition Probability of the Markov Chain

Table 8. Markov chain probability transition matrix (1995–2018) (in %).

		2018
	LULC Class	Artificial Surfaces	Arable Land	Permanent Crops	Pastures	Heterog. Areas	Forest	Water Bodies
1995	Artificial surfaces	0.97	0.00	0.01	0.01	0.00	0.00	0.00
	Arable land	0.01	0.62	0.19	0.11	0.01	0.07	0.01
	Permanent crops	0.00	0.07	0.83	0.02	0.06	0.02	0.00
	Pastures	0.01	0.19	0.07	0.46	0.02	0.25	0.01
	Heterog. areas	0.00	0.03	0.02	0.06	0.76	0.12	0.01
	Forest	0.00	0.01	0.00	0.03	0.08	0.87	0.01
	Water bodies	0.01	0.00	0.01	0.00	0.00	0.03	0.95

presents the temporal transition probability matrix for the overall period (1995–2018), showing the probability that each LULC type transitions to another (the off-diagonal entries indicate the transitions probability from one LULC class to another). As expected, there is a small probability that artificial surfaces will transition to another use (3%), and the same is observed with the water bodies (5%). For forest areas, there is a probability of an approximately 13% change in use, specifically, approximately 8% towards heterogeneous areas.

With regard to croplands, there is a clear tendency for each type of agricultural use to transition to another agricultural use or to forest. For example, pastures and arable land are the classes that are most likely to transition. Pastures have a 25% probability of transition to forest and a 19% probability of transition to arable land. Arable land presents a 19% probability of transition to permanent crops and an 11% probability of transition to pastures. Heterogeneous areas have a 12% probability of transition to forest. From all croplands, the low probability of permanent crops to transition (about 17%) and the low probability of heterogeneous areas to transition (about 24%) to other LULC classes are noteworthy.

3.6. Classification States

The analysis in Section 3.5 showed that each state has a positive probability to achieve any other state. Therefore, all states of the DTMC communicate, and the matrix can be cataloged as irreducible. As there are no transition or absorbing states, the matrix is also recurrent; i.e., all states are recurrent. Finally, as we always obtain $k=1$, the recurrent states are all aperiodic. Together, since all states are recurrent and aperiodic, communicating the matrix can be cataloged as ergodic. Thus, having a transition probability matrix, the mean first passage time was evaluated.

Table 9. Mean first passage time.

	2018
1995	Artificial Surfaces	Arable Land	Permanent Crops	Pastures	Heterog. Areas	Forest	Water Bodies
Artificial surfaces	0.0	61.3	62.0	53.2	55.4	46.1	145.1
Arable land	232.3	0.0	23.5	26.8	27.0	17.2	134.2
Permanent crops	233.3	27.9	0.0	33.1	24.2	19.8	135.7
Pastures	232.1	26.7	34.9	0.0	25.5	12.1	132.6
Heterog. areas	233.8	36.0	42.1	29.8	0.0	12.7	132.4
Forest	236.6	40.0	46.6	33.2	20.5	0.0	130.7
Water bodies	218.5	58.4	58.3	52.7	46.5	31.6	0.0

gives us the mean first passage time between 1995 and 2018. It is clear that a high number of steps are expected to be necessary for transition to artificial surfaces and water bodies, suggesting that both of these LULC classes are almost stable. Faster transitions are, in all cases, towards a forest state.

Having the recurrent states, we calculated the mean recurrence time, i.e., the expected number of steps to return to the initial state. This measure is narrowly related to the mean first passage and states. Even in our case, a returning probability of 1 does not signify that the predictable return time is fixed. The results showed that all the LULC classes represent positive recurrent states (

Table 10. Mean recurrence time.

LULC Class	1995–2007	2007–2010	2010–2015	2015–2018	1995–2018
Artificial surfaces	22.8	122.0	11.2	76.8	8.62
Arable land	37.7	33.8	5.1	25.1	12.80
Permanent crops	56.4	9.2	6.2	1.5	7.16
Pastures	83.8	99.1	8.7	12.7	17.98
Heterog. areas	21.4	16.5	8.0	12.8	6.74
Forest	7.5	7.8	4.0	9.8	2.97
Water bodies	1.4	1.5	16.1	46.9	7.90

). Thus, if the system initiates in state $s_i$, in determining the number of steps necessary to return to $s_i$ for the first time, it is obvious that it will eventually return to the initial sate. This is because, in the first step, the system either stays at $s_i$ or changes to other state $s_j$, and from that state $s_j$, one will ultimately reach $s_i$ since the chain is ergodic.

From Table 10, it is noticeable that, in the first period (1995–2005), the LULC classes with the lowest mean recurrence time are forest (7.5) and water bodies (1.4), which denotes the resilience of these two LULC classes. The next period (2005–2010) is similar to the previous one, but we start to see the emergence of permanent crops (9.2). The last two periods emphasize the importance of forest and permanent crops in the study area. However, the roles invert themselves, since in 2015–2018, permanent crops achieved a major role in the landscape of Beja, and this may soon become an absorbing state. As for the overall period from 1995 to 2018, we can confirm the forest resilience with a lower value (2.97), whereas the different changing cycles slightly reduce the importance of permanent crops, especially due to the dynamics between 1995 and 2007.

Table 11. Steady sates (%).

LULC Class	1995–2007	2007–2010	2010–2015	2015–2018	1995–2018
Artificial surfaces	0.04	0.01	0.09	0.01	0.11
Arable land	0.03	0.03	0.20	0.04	0.07
Permanent crops	0.02	0.11	0.16	0.67	0.13
Pastures	0.01	0.01	0.11	0.08	0.05
Heterog. areas	0.05	0.06	0.12	0.08	0.14
Forest	0.13	0.13	0.25	0.10	0.33
Water bodies	0.72	0.65	0.06	0.02	0.12

shows that, on average, and after some time, following the 2015–2018 trend, 67% of the area will be permanent crops. Likewise, considering the 2010–2015 trend, about 20% of the territory will be arable land and about 25% forestland. In the former periods of 1995–2005 and 2005–2010, it could be expected that 72% and 65% of the Beja district area will be water bodies, respectively. Lastly, the analysis of the overall period (1995–2018) allows us to embrace the different dynamics observed, avoiding trends too ephemeral to be considered in the long term and reconfirming forestland as the most likely class to occur.

3.7. Simulating Future States

As stated before, in our recurrent positive (i.e., having a stationary distribution) and aperiodic matrix, regardless of the initial probabilities, as the time step increases towards infinity, the chain probability distribution converges. In this case, the chain has a limiting distribution, i.e., a stationary distribution. In Figure 7, we inspect a random simulation of $p_{ij}^{\left(n\right)}$ through time and show that, taking into account the 1995–2015 probability vector (Figure 7a), states converge to two vectors. The probability of being in States 1 (artificial surfaces) and 7 (forest) converges to higher values than the others LULC classes. Over the long term (e.g., after 48 time steps), the increasing probability of becoming forest overcomes the decreasing probability of becoming artificial surfaces.

As for Figure 7b, where we focus on the most recent trend (2015–2018), the situation is somewhat similar. While the system evolves to infinity, the probability vector converges to a vector assigning most of the probability to State 3 (permanent crops). In this case, the probability of becoming permanent crops overpasses that of becoming artificial surfaces only after 28 time steps. In addition, the decreasing rate of the probability of becoming artificial surfaces is much higher in b) than in a). Lastly, it was confirmed that the empirical probabilities behave in the same way (Figure 8).

Therefore, the future state simulation results showed that, in 2040, about 25% of the territory can transition use, and about 21% of that territory corresponds to agricultural land. Accordingly, the major changes concern cropland areas, with a predicted decrease of about −3.67% in 2040 (total 60.3%) (Figure 9). However, this situation hides very different behaviors for each cropland class, since arable land (−3.81%), pastures (−1.35%), and heterogeneous areas (−1.02%) diminish. In contrast, and following the trend of past years, permanent crops will see its importance continue to rise (+2.93%). The remaining LULC classes will continue expanding until 2040, where forest class will increase by about 2.60%, artificial surfaces will increase by about 0.40%, and water bodies will increase by about 0.65%. Figure 9 shows the spatial distribution of LULC by 2040.

In Figure 10a, we have the uncertainty map for the 2040 LULC simulation, with values ranging from 0.03 to 0.53. Crossing these values with LULC classes shows us where the model performs better. Thus, we have two classes that are very resilient and are therefore modeled with a low degree of uncertainty: artificial surfaces and water bodies, both with an average uncertainty of 0.03. Next, also with low values, are the two most dynamic and important classes of the study area, namely, forest (0.13) and permanent crops (0.17). Finally, three classes remain, which, due to their complex structure and overlapping, achieve the worst values: heterogeneous areas (0.23), arable land (0.38), and pastures (0.53). In Figure 10b, we have a spatial notion of the entropy distribution with values oscillating between 0.28 and 1.95. Spatially, both calculations, uncertainty and entropy, show a high rate of agreement.

4. Discussion and Conclusions

In recent decades, the abandonment of agricultural land has been an international trend. In Portugal, this has occurred, confirming marked changes in the agricultural landscape [30,31]. Therefore, in this study, we analyzed the spatiotemporal changes that occurred in the largest district of Portugal, Beja, which has an agricultural economic background, and provided a trend projection for the year of 2040. Accordingly, the monitoring of the dynamics of LULC at the regional level revealed noteworthy LULC dynamics intrinsically associated with distinct agricultural practices of food production and distribution.

4.1. About the LULC Changes Detection

The degree of transition between the LULC classes in the Beja district was significantly different during each period under analysis. The following findings are worth noting:

(i)

In 23 years, about 28% of the Beja district experienced dozens of inter-class transitions, but the most significant changes were in the croplands. Nevertheless, croplands remain the predominant LULC class in more than half of the Beja district (in 2018 about 64%). Particularly, there was a predominance of arable land that was dispersed in a great part of Beja, despite the significant decrease in this practice between 1995 and 2018 (−5.9%). Despite the decrease, the results confirmed the importance of cereal production in this region, as in the case of wheat [43,44]. Likewise, heterogeneous areas experienced a continuous decline (−2%), but even so this class represents the second largest cropland in the district, which confirms the relevance of the Montado system and the cork oak exploitation in the agroforestry areas of the district.

(ii)

Of all the agricultural classes, permanent crops experienced a growing tendency. Thus, this cropland gained prevalence after a certain period and simultaneously experienced some loss (+5.13% between 1995 and 2018). In spite of the reduction in the territory, the substantial increase in permanent crops suggests an emphasis on intensive farming without fallow, proved by the enlargement of the olive grove plantation in this region and the increase in olive oil production [43,44]. These results are consistent with other studies, despite the methodological approach differences [31,32,45].

(iii)

As for the overall period, forests experienced more gains than losses (+4.31% between 1995 and 2018). There was a clear tendency of agricultural land loss to the detriment of forests, these results being in line with what is happening on a national level observed in other studies [30,46]. Additionally, about 63% of the total changes experienced by forests was attributed to a swap type, which may indicate that the Beja district simultaneously experienced forms of reforestation and deforestation between 1995 and 2018.

(iv)

Lastly, artificial surfaces and water bodies experienced more gains than losses and showed a tendency to persist. A main cause for these observed gains is spatial conversions, such as the expansion of the Beja airport area and the construction of new roads in the case of artificial surfaces [47], and the Alqueva dam construction in 2002 in the case of the water bodies [31]. Nonetheless, artificial surfaces only increased about 0.45%, between 1995 and 2018, which indicates that the increase in urban areas was not remarkable in the Beja district, thereby counteracting the Portuguese trend towards an urban area increase [46].

Accordingly, the behavior of the inter-class transitions was significantly different between periods, and explicitly revealed that arable land, pastures, and forest were the most dynamic LULC classes, having experienced both damage and recovery over the years. Hence, the drivers that promote the LULC gains are likely to differ from the drivers that lead to LULC losses [12]. On the one hand, the main cause for cropland changes could be the competition between LULC classes, since such competition may restrict the use of available land even when this land has suitable determinants for a specific LULC class, which would explain higher observed losses compared with gains [12]. In addition, the political measures of incentives for the increase in technology, for example, the almost total abandonment of traditional agriculture or the high level of subsidies to agriculture, could be explanatory factors for the agricultural LULC changes observed [48]. Accordingly, the LULC periodical fluctuations were depolarized by certain drivers, such as political decisions, economic development planning, or even weather conditions, showing the importance of the analysis of the possible drivers of changes in future work [46].

4.2. About the Dominant LULC Transitions

The analysis of the differences between the observed LULC and the expected LULC gains and losses under a random process of changes showed that, of the several observed inter-class transitions, few experienced dominant signals of change between 1995 and 2018 (Figure 11), namely:

(i)

transitions of arable land to permanent crops and to pastures, which accounted for about 5% and 2.9%, respectively;

(ii)

transitions of pastures to forest and arable land, which accounted for about 3.5% and 2.7%, respectively;

(iii)

transitions of forest to heterogeneous areas, which accounted for about 2.1%.

Accordingly, the arable land tendency to lose, especially for the permanent crops class, causes some concerns. The Beja district located in the Mediterranean region is a suitable region to produce olive oil mainly due to favorable weather conditions. Moreover, the positive market environment prevailed with the gradual increase in global olive oil consumption in the past decade, and the super-intensive olive grove plantation also increased to cover the demand [45]. Thus, traditional olive grove plantations are becoming unusual, mainly because the super-intensive cropland production profitability is higher and because of Common Agricultural Policy (CAP) incentives [49]. While this type of agricultural practice will start producing in the 3rd year after planting, the traditional practice can occur after one decade. In addition, mechanization in the super-intensive olive grove is less labor-dependent, meaning a greater capacity for profit. Hence, this all points to one conclusion: the high increase in super-intensive olive crops to produce olive oil is highly associated with the expansion of permanent crops, which is at the expense of arable land.

Thus, the Beja district is experiencing a new landscape dynamic, caused not only by the systematic increase in permanent crops areas but also by the decrease of arable land that produces a historical agricultural product: wheat (the grain of civilization). Indeed, over the last century, Beja consistently showed remarkable importance in wheat production, mainly explained by: (i) the type of soil in Beja, “Barros de Beja” (Beja clays), which is a dense and deep soil that is particularly suitable for this cereal production; (ii) the type of land exploitation (latifundium) [50]; and (iii) the weather conditions (i.e., a temperate climate with Mediterranean characteristics), as the wheat crop has a relative tolerance to water deficiency [51].

Thus, these results indicate that, although arable land predominated in the Beja district between 1995 and 2018, this cropland experienced a significant decrease due to the competition between LULC classes. Accordingly, the transition from arable land to permanent crops was relatively constant over time, and it was more pronounced in the last period (2015–2018), presumably triggered by the substantial increase of foreign investment (especially from Spanish investors) and the strong exploitation of water resources, after the construction of the Alqueva dam that filled its reservoir only 10 years after construction began (in 2012) [45]. Thus, it seems essential to assess the impact of this systematic transition, since everything points to an increasing super-intensive olive grove cultivation system that will result in a negative impact on this land in the future [52].

In addition, arable land also experienced systematic transitions to pastures. This change was especially strong in the first period (1995–2007), possibly due to changes in the CAP, which led to a decrease in the areas under arable crops, in particular cereal production, and an increase in fallow land maintenance and natural pasture. These policies arose out of a concern related to land preservation in the context of rural planning, resulting, among other actions, in reduced dryland cereal cultivation areas, due to the recognition of the negative impact on soil depletion and erosion resulting from previous agricultural campaigns [53] and an international acknowledgement of the role of such areas in mitigating climate change [54,55].

On the contrary, Beja also experienced a systematic transition from pastures to arable land. This dynamic LULC change was substantially stronger in the first period (1995–2007) and partly contradicts the prior interpretation of the possible factors responsible for such changes. Therefore, following the abovementioned explanations, we believe this transition can also be linked with the implementation of the CAP, which favors fallow land maintenance. Indeed, the arable land class mainly comprises a land rotation system between annually harvested crops (rain-fed or irrigated) and fallow plots, which in our study area mainly represents dry-fed and rain-fed cereal production. However, according to the COS map nomenclature, this class can also comprise fallow land for a maximum of five years. Additionally, agricultural pasture systems have limitations regarding their production irregularity over a given year due to the irregularity of rainfall [56,57] and may not be economically viable for farmers if not subsidized [50,58], which may partly justify these changes in the landscape.

In the Beja district, a systematic transition from pastures to forestland was also found, and this was most significant in the first period (1995–2007) and was probably associated with afforestation policy (the 1992 CAP reform) and EU subsidies (Agenda 2000), which promoted the conversion of pastures in Eucalyptus and Pinus forests [50] for conservation of the land’s structure and functions. Lastly, the systematic transition from forest to heterogeneous areas, which occurred significantly between the third period (2010–2015), is probably related to forest fires, as reforestation took place after this event in the form of Montado to combat desertification [59].

For instance, the examination of the systematic and random transition showed that some of the major LULC transitions, such as heterogeneous areas and arable land transitioning to forest (2.53% and 1.78%, respectively) and the heterogeneous areas transitioning to pastures (1.18%), were random or partially systematic. In addition, it is interesting that heterogeneous areas did not experience systematic losses, which indicates a strong signal of persistence in this LULC class. Thus, the Beja district still retains important traditional agricultural areas, such as the Montado, an agrosilvopastoral system of global importance.

4.3. About the LULC Future Changes

The transition probability simulation between uses for the year 2040 allowed us to understand that imbalances can arise in the future, highlighting that about 25% of the territory is predicted to transition, and about 21% of that territory corresponds to agricultural land. These results are a warning that the Beja district will experience strong landscape dynamics in agricultural land, emphasizing the importance of measures in mitigating the negative effects of such strong signals of LULC transitions [12].

Additionally, the DTMC model showed a high agreement value, proving to be a useful mechanism for monitoring and evaluating LULC changes, with a trend projection with high applicability and flexibility. Therefore, the main advantage of using Markov chains is their mathematical logic that enables swift applications, as well as their possible interaction with other techniques such as remote sensing and GIS, both of which commonly work with matrices. However, a disadvantage of this model is that it does not explain the facts, because the process only takes into account the time variable. Thus, predictions made through this model mostly indicate changes in the state of variables in time and not in space. Nevertheless, the stationary probabilities are useful as indicators of the landscape shares of each LULC class and provide highly valuable information in decision-making processes.

4.4. Main Conclusions and Future Studies

The purpose of this study was to analyze the spatiotemporal evolution of LULC in Beja from 1995 to 2018 and simulate the future development for the year 2040. This analysis indicates that in 23 years there were significant changes that deeply modify the Beja district landscape, mainly due to two opposite phenomena: the intensification of agricultural practices and the afforestation of agricultural lands. The observed LULC patterns reflect the last decade’s structural changes in the Portuguese economy and society and confirm the impact of political decisions and different agricultural activity management, which led to the loss of a significant amount of agricultural heritage. In fact, this regional landscape and the traditional features of agricultural land are facing an uneven transformation dynamic that does not take into account either the balance or the environmental enhancement of the local diversity agricultural systems [49].

The acquired spatiotemporal knowledge regarding the agricultural LULC changes allowed us to obtain a general picture of the main transitions of a region of Portugal that has more importance in the regional agricultural production sector, particularly in the production of cereals, such as wheat, barley, oats, and other products such as olive oil and cork. An understanding of the systematic transitions and the locations thereof provides useful insight into ecosystem service management and contributes to future sustainability strategies. Nonetheless, this statistical approach is exclusively based on an analysis of the LULC transition matrix, and the interpretation of the LULC transitions is based on a sequence of historical determinism documented in prior investigations (e.g., [49,57,60]), so any factor mentioned to rationalize the systematic transition is based only on assumptions of historical fact. Therefore, additional research is necessary to assess the many processes that create these LULC changes, which have been linked to dozens of possible drivers and their different levels of feedback. Still, providing an understanding of the inter-class transitions that have occurred as systematic changes has proven to be a straightforward, effective approach that allows us to differentiate the predominating quantity transitions from those that experienced a single abrupt episode process of change influenced by factors that occur suddenly.

The information obtained contributed to sciences related to the earth and human society on a variety of levels, providing effective information for economic and environmental monitoring and evaluation, resulting in informed policy decisions and in the design of anticipatory measures response to change. Furthermore, in this exploratory study, we used COS maps, official Portuguese LULC data produced by a governmental institution (DGT), and these data were advantageous because of their higher cartographic detail and were characterized by temporal continuity. This highlights the importance of having an entity responsible for producing LULC thematic cartography with regularity and reliability and making it available free of charge, as it can be used for multiple purposes, such as territorial planning and management, ecosystem service mapping, and the identification of areas vulnerable to rural fires (e.g., [34,61,62,63]).