Simulation of land use/land cover changes and urban expansion in Estonia by a hybrid ANN-CA-MCA model and utilizing spectral-textural indices

Over the recent two decades, land use/land cover (LULC) drastically changed in Estonia. Even though the population decreased by 11%, noticeable agricultural and forest land areas were turned into urban land. In this work, we analyzed those LULC changes by mapping the spatial characteristics of LULC and urban expansion in the years 2000–2019 in Estonia. Moreover, using the revealed spatiotemporal transitions of LULC, we simulated LULC and urban expansion for 2030. Landsat 5 and 8 data were used to estimate 147 spectral-textural indices in the Google Earth Engine cloud computing platform. After that, 19 selected indices were used to model LULC changes by applying the hybrid artificial neural network, cellular automata, and Markov chain analysis (ANN-CA-MCA). While determining spectral-textural indices is quite common for LULC classifications, utilization of these continues indices in LULC change detection and examining these indices at the landscape scale is still in infancy. This country-wide modeling approach provided the first comprehensive projection of future LULC utilizing spectral-textural indices. In this work, we utilized the hybrid ANN-CA-MCA model for predicting LULC in Estonia for 2030; we revealed that the predicted changes in LULC from 2019 to 2030 were similar to the observed changes from 2011 to 2019. The predicted change in the area of artificial surfaces was an increased rate of 1.33% to reach 787.04 km2 in total by 2030. Between 2019 and 2030, the other significant changes were the decrease of 34.57 km2 of forest lands and the increase of agricultural lands by 14.90 km2 and wetlands by 9.31 km2. These findings can develop a proper course of action for long-term spatial planning in Estonia. Therefore, a key policy priority should be to plan for the stable care of forest lands to maintain biodiversity.

Spatiotemporal LULC changes and urban expansion are widely modeled to detect the consequences of longterm interactions between humans and the environment (Liu et al., 2017a, b), contributing to sustainable development (Abbas et al., 2021). For this purpose, location-based models are commonly utilized (Tong & Feng, 2020) to simulate the LULC dynamics. Cellular automata (CA) is the most widely used location-based model, especially for simulating the spatiotemporal evolution of LULC and urban expansion . The simplicity and explicit representation of LULC changes (Grinblat et al., 2016) make the CA a standard model for simulating such dynamics. While the CA model assumption is based on the effects of past changes on future transitions (Santé et al., 2010), these background conventions limit the model's ability to simulate the complex nature of LULC over time realistically. Besides, the transition rule, which extracts the state of a cell over time  and is an essential function in the CA model, varies in terms of geographical regions and neighborhood interactions. Integration of CA with other statistical and geospatial models can enhance the CA model's predictability (Aburas et al., 2017;Santé et al., 2010). Commonly used reinforcement methods are the logistic regression model (Guan et al., 2019;Jokar Arsanjani et al., 2013;Siddiqui et al., 2018), Markov chain analysis, MCA (Deep & Saklani, 2014;He et al., 2018;Saadani et al., 2020), agent-based model (Liu et al., 2020a, b;Mozaffaree Pour & Oja, 2021;Zeng et al., 2018), and artificial neural network, ANN (Abbas et al., 2021;Xu et al., 2021).
MCA generates the area and probability of transition (Mondal et al., 2016;Xu et al., 2021) for multiple LULC types by applying the past, present, and future states of the cells; chain here refers to the stochastic process of distribution of transition probability on the future time (Mondal et al., 2016). While MCA lacks representation of the effects of LULC drivers  and spatial dependency (Ullah et al., 2019), the ANN algorithm is an enhanced method compatible with CA-MCA to address the spatial probability of changes. ANN that imitates the human brain's neurons is a machine learning algorithm (Wang, 1994) and can be trained to estimate the probability of occurrence from non-linear functions dependent on many input layers (Liu et al., 2017a, b). The ANN models can learn how neurons' weights change and calibrate during the training process to achieve the desired output (Gharaibeh et al., 2020;Saha et al., 2021). Hence, constructing a hybrid ANN-CA-MCA is beneficial to simulate a more realistic projection of LULC.
Consequently, LULC change modeling requires a deep understanding of the driving factors , historical changes , and environmental predictors (Msofe et al., 2019). Extensive research has modeled the effects of physical factors, including elevation, slope, aspect (Bonilla-Bedoya et al., 2020;Feng & Qi, 2018;Girma et al., 2022;Xu et al., 2021), and soil texture (Han et al., 2015;Tian et al., 2011) on LULC changes. Several attempts have explored the importance of proximity factors on the LULC changes. Among others, distance to roads (Rahnama, 2021;Ullah et al., 2019), distance to water bodies and rivers (Paterson et al., 2015;Sohl et al., 2012), and distance to city (Abbas et al., 2021;Gharaibeh et al., 2020;Patra et al., 2018) were the most important factors. A growing body of literature has investigated the environmental impacts such as temperature on LULC transitions (Fattah et al., 2021;Li et al., 2016;Rahaman et al., 2020;Ullah et al., 2019). Therefore, understanding the changes and dynamics of LULC is challenging and requires modeling techniques utilized with spatiotemporal data (Jokar Arsanjani, 2018;Singh et al., 2015). Likewise, selecting proper predictors which detect the nature structure of LULC and determine the pattern of changes are essential in LULC changes modeling.
The spectral-textural data from satellite sources provide reliable remote sensing information for monitoring and mapping the LULC (Zhang et al., 2014) and landscape morphology over time. The Gray Level Co-Occurrence Matrix (GLCM) indices provide spatial interactions between each pixel and its neighbor (Löfstedt et al., 2019). Haralick et al. (1973) introduced the GLCM as a quantifying method for detecting the neighboring pixel relationship in an image. Previous studies investigated the importance of spectral-textural data for LULC classification. Kupidura (2019) studied different techniques and satellite imagery tools analyzing the texture features for LULC classification and indicated that spectral and textural GLCM data increase the overall accuracy of the results. Hurskainen et al. ( 2019) applied the spectral-textural features to improve the LULC classification accuracy in a heterogeneous landscape and demonstrated the potential of these indices for LULC mapping. Others also examined spectral and textural features as a beneficial tool for classifying LULC types like agricultural lands (Okubo et al., 2010;Wikantika et al., 2004) and urban areas (Chaturvedi et al., 2020;Mhangara & Odindi, 2013).
While determining spectral-textural indices is quite common for LULC classifications, utilization of these continues indices in LULC change detection and examining these indices at the landscape scale is still in infancy. Zhang et al. (2019) have established a combined textural and spectral framework to obtain a more accurate LULC classification map applicable for the CA-MCA model to interpret the past trends of a coastal river basin in China and simulate its future. Halmy et al. (2015) predicted the LULC changes in Egypt using the CA-MCA model utilized with the textural indices to improve the LULC classification accuracy. Ghasemi et al. (2021) applied the GLCM of texture features to increase the accuracy of the classification of LULC in the city of Tabriz, Iran, and found that using the eight GLCM texture feature improved the accuracy of the results by approximately 9%. Others also examined spectral and textural features as a beneficial tool for classifying LULC types like agricultural lands (Okubo et al., 2010;Wikantika et al., 2004) and urban areas (Chaturvedi et al., 2020;Mhangara & Odindi, 2013).
A considerable amount of literature has been published on the importance of patch-based landscape metrics in modeling the spatial characteristics of LULC and urban expansion from time-series data (Das & Angadi, 2021;Hesselbarth et al., 2019;Jiao et al., 2015;Mirsanjari et al., 2021;Yang et al., 2016a, b;Zhou et al., 2014). However, an essential aspect of spectral-textural indices is their pixel-based capabilities and potentials in detecting the discrete cells of multiple LULC transitions over time and very little is currently known about applying the pixelbased predictors by a cell-based model to predict LULC changes.
In this work, we hypothesized that spectral and textural properties of landscapes, which can be captured with optical imagery, provide an adequate proxy for LULC transitions. The primary objectives of this study were (i) to identify the importance of spectraltextural predictors on LULC change detection and urban expansion between 2000 and 2019; (ii) to quantify the future LULC demand quantities; (iii) to estimate the spatial probability of occurrence of multiple LULC types; (iv) to predict LULC changes by 2030 in Estonia; and (v) to monitor the urban expansion in the main cities of Tallinn, Tartu, and Pärnu. We performed multiple regression analysis to investigate the importance of spectral-textural predictors on LULC change detection and urban expansion and their interrelation. Then, we utilized spectral-textural indices in the hybrid ANN-CA-MCA simulation model to explore the spectral-textural indices' capabilities in representing the past trends of LULC changes and modeling the future LULC in Estonia. Hence, by applying remote sensing indices, the hybrid ANN-CA-MCA model contributed to addressing the heterogeneity in space, determining the spatiotemporal transitions of LULC, improving the model accuracy, and enhancing the simulation of LULC and urban expansion. It is worth pointing out that the present work provided the first comprehensive simulation for the future changes of LULC in Estonia.

Study area
The study area is Estonia (Fig. 1), covering 45,338 km 2 of northeastern European lands. This study area was chosen for its apparent LULC dynamics and population changes over two decades. Estonia experienced reforestation at the beginning of the twenty-first century (Poska et al., 2018) and loss of agricultural lands (Kaasik et al., 2008). At the same time, based on the statistics of Estonia (Statistical database, 2022), the country's population has decreased by 11% from 2000 to 2019. The population minimum occurred in 2015; since then, slight growth has been noticeable. In contrast to the overall country's population decrease, the internal migration increased the population in the immediate vicinity of major cities, including Tallinn, Tartu, and Pärnu, and respectively, new built-up areas shaped a scattering pattern. In 2019, Tallinn had almost 9396 internal migrants, more than 2000 (increased by 86%). In Tartu, 3117 people immigrated from the other territories of Estonia with an internal migration rate of 43.70% between 2000 and 2019. Also, 1469 migrants were settled in Pärnu in 2019 (Statistical database, 2022). Conversely, the rest of the regions (majority of the territory) have decreased population. As a result, the meaning of the rural lifestyle has changed in Estonia within the last decades. Approximately two-thirds of the population in Estonia is settled in municipalities considered administratively urban, while lifestyle of most formally "rural" municipalities around Tallinn, Tartu, and Pärnu is according to an urban or suburban (Oja, 2020). Dataset and data pre-processing LULC data LULC data were extracted from the data provided by Parente et al. (2021) for continental Europe. This dataset was based on the ensemble machine learning method using the samples from primary European datasets, including the LUCAS and CORINE land cover. Then we extracted the data to Estonia's extent and used the level 1 class consisting of the five main categories of LULC. The data was reprojected using the Lambert Conformal Conic projection (Estonian national grid of 1997, EPSG 3301).

Optical remotely sensed data
We applied the cloud computing platform Google Earth Engine provided by Google LLC (Gorelick et al., 2017) as it combines freely available remote sensing datasets with rich geospatial processing capabilities. Google Earth Engine was employed to derive a set of remote sensing-based indicators of landscape physiognomy from cloudless Landsat mosaics. We focused on three anchor years to study LULC changes in Estonia: 2000 and 2011 with Landsat-5 data, and 2019 with Landsat-8 data (Supplementary Tables 1-3). Median surface reflectance was used for the summertime period of these years, and cloudy pixels were filtered out using the available quality assessment band. Due to the cloudy atmospheric conditions in 2000, the respective cloudless Landsat-5 mosaic was compiled for 1999-2000 imagery. Textural metrics based on the GLCM (Conners et al., 1984;Haralick et al., 1973), focal statistics, and local indicators of spatial autocorrelation (Anselin, 2010) performed well in previous studies on the physiognomy of landscapes in Estonia (Karasov et al., 2022) and urban classifications (Sapena et al., 2021;Wurm et al., 2017). A 21 × 21 kernel size was previously shown as optimal for GLCM-based urban classification (Kuffer et al., 2017); a smaller size captures local spectral variations, while the bigger size is not computationally adequate and starts overlooking landscape diversity. The following groups of indicators were chosen, as shown in Table 1.

Data analysis
In this study, a CA model was performed on the dynamics of LULC categories which affect each other and interact simultaneously to simulate a more realistic future LULC. MCA and ANN were applied in this research to enhance the mining methods of transition rules. To construct this hybrid ANN-CA-MCA model, the data was processed and analyzed as shown in the schematic framework in Fig. 2 and described in the following subsections:

Spatiotemporal analysis of LULC and urban expansion
The LULC data was extracted for 2000, 2011, and 2019. Level 1 class (five main LULC types) was used: artificial surfaces (including the urban fabric, construction and industrial sites, mineral extraction sites, roads, and dump sites), agricultural lands (consisting of non-irrigated arable land and pastures), forest lands (including the broad-leaved forest, coniferous forest, natural grasslands, transitional woodland-shrub, and sparsely vegetated areas), wetlands (covering the inland wetlands and maritime wetlands), and water areas (consisting of the watercourses and water bodies) (Parente et al., 2021). Then the LULC changes were calculated and mapped from 2000 to 2011 and 2011 to 2019 for further analysis. The whole surrounding water areas consisting of the Gulf of Finland, the Baltic Sea, and Lake Peipsi were excluded from the calculation and simulation.
It is essential to note that urban fabric, industrial, and mineral extraction sites were classified as a single main class of artificial surfaces on a country-wide scale. Thus, to capture and analyze the urban expansion in the three main cities of Tallinn, Tartu, and Pärnu, which are far from the vast mining and industrial sites in the northeast of Estonia, urban areas were identified as artificial surfaces.

Multiple regression analysis
In this study, we performed the multiple regression analysis using the MULTIREG tool in IDRISI software, which adopts a least-squares approach to multiple regressions (Eastman, 2016) and provides information about the interrelations between the dependent variable and the independent variables. We applied multiple regression analysis twice: first, to explore the relative influence of the spectraltextural indices (independent variables) on LULC changes (dependent variable) for the 2000-2011 and 2011-2019 periods; second, to analyze the relationships between artificial surfaces' probability of occurrence (dependent variable) and 19 spectraltextural indices (independent variables) to estimate the importance of spectral-textural indices on artificial surfaces' probability of occurrence in 2011, 2019, and 2030.

Projection of LULC demand and transition probability by Markov Chain Analysis
MCA is a random mathematical tool for finding the transition probability, demand quantity, area (Aburas et al., 2017;Xu et al., 2021), and simulation constraints (Zhao et al., 2020) for each LULC type. MCA took advantage of two timestamps spatial data to learn the transition probability (Gharaibeh et al., 2020) and calculate the future transition probability (Zhou et al., 2020). The final allocations of LULC in the simulation are the results of the transition probability of two where S t+1 and S t represent the LULC at time t and t + 1 , and P represents the state of the transition matrix (Zhou et al., 2020). The transition matrix reflects the probabilities of changes by each LULC category will change to the other categories and the area represents the number of pixels predicted to change for each LULC category.
To enhance the simulation accuracy (Liang et al., 2018a, b;Liu et al., 2017a, b), actual LULC areas equal to the demand's quantity in 2011 and 2019 were employed.

Estimation of probability of occurrence by artificial neural networks
ANN was utilized to compute the non-linearity of spatial probability of occurrence considering the predictors of LULC changes over time by a self-learning algorithm suitable for the CA model simulation. A (1) S t+1 = P × S t weight is allocated to the neurons in the hidden layer based on the received signals and the information of input neurons on the training procedures. The hidden layer neurons' range was recommended to set between at least (2n/3) and (2n + 1) neurons (n represents the number of input neurons) (Wang, 1994;Zhao et al., 2020). To act in a non-linear way, a sigmoid activation function was added to the ANN algorithm (Liu et al., 2017a, b). The probability of occurrence was the result of the ANN process, and the final value of each cell indicated the probability of occurrence of each LULC type.
Here, 19 spectral-textural indices were selected as predictors of LULC in each simulation process. The ANN input variables were considered neurons of spectral-textural indices and were set to 19 neurons which were normalized prior to training. Hidden layer neurons were set to 19, and the output neurons correspond to the five categories of LULC. A total of 332,348 samples (17,492 samples for each input neuron) were proportionally selected to train the ANN model. The ANN algorithm was self-adaptive, so the learning rate was not preselected.

LULC and urban expansion simulation by cellular automata
In this study, the CA-based future land-use simulation (GeoSOS-FLUS), free software, was applied to simulate the future projection of LULC in Estonia. As a bottom-up model, the CA model consists of discrete cell space, cell state, neighborhood, and transition rules developed (Mozaffaree Pour & Oja, 2021;Santé et al., 2010;Zhou et al., 2020). State of the cells is defined as the qualitative values representing the LULC categories. While the cell state is affected by its previous state, it will change regarding the implementation of different transition rules (Mozaffaree Pour & Oja, 2021). Neighborhood effects measure the spatial interactions of LULC types and the surrounding neighbors (Omrani et al., 2017;Xu et al., 2021). So, the neighborhood's size and weight are essential for the CA model. The weights of the neighborhood that are variant through time were calculated based on the evolution ratio of each LULC type.
Here, the initial configuration of the CA model was set as follows: a discrete cell space with a spatial resolution of 50 × 50 m 2 described the quantity and degree of development in each LULC type. The Moore neighborhood of 7 × 7 grid units was applied as the neighborhood size implemented in much research ( Huang et al., 2018( Huang et al., , 2021Shafizadeh-Moghadam et al., 2017;Xu et al., 2021;Zhou et al., 2017). Defining the transition rules is a critical step in the CA model and plays a significant role in the accuracy of simulation results (Liu et al., 2017a, b).
Here, multiple LULCs compete with each other, and subsequently, the LULC allocation and quantities of changes are defined by the transition rules (Liang et al., 2018a, b). To this end, the self-adaptive inertia addressed the interactions and competitions among the LULC categories, and the combined probabilities of change and occurrence help the CA model simulate the future projection of LULCs (Li et al., 2017a, b;Yang et al., 2020). It is notable to mention that spectral-textural indices provide a valuable local understanding of the LULC dynamics in detail scale and trends of LULC changes.

Simulation validation
The kappa coefficient, overall accuracy, user's accuracy, and producer's accuracy were calculated to identify the simulation accuracy of LULC results and model performance for 2011 and 2019. The kappa coefficient calculates the agreement between the observed and expected correct in the simulation and LULC maps (Cohen, 1960) and measures the pattern's consistency between the simulation and LULC maps (Dan-Jumbo et al., 2018;Liu et al., 2017a, b;Xu et al., 2021). To quantify the agreement between the simulated map and reference map, the kappa coefficient is calculated as follows (Eq. 2) (Cohen, 1960): where P 0 indicates the proportion of units (measurement of area, pixel, or polygon counts (Stehman & Czaplewski, 1998)) in which the maps agreed and P e is the proportion of units for which agreement is expected by chance. The kappa coefficient = 1 indicates the complete consistency with the reference map, the kappa coefficient = 0 expresses no agreement, and the upper limit occurs when there is a perfect agreement between maps (Cohen, 1960).
The overall accuracy assesses the overall proportion of correctly simulated area (Stehman & Czaplewski, (2) kappa = P 0 − P e 1 − P e 1998) or percentage of corrected simulated units by LULC unit samples (Abdelkareem et al., 2018;Rwanga & Ndambuki, 2017) and derived from the samples in the confusion matrix (Olofsson et al., 2014), calculated based on Eq.
(3) (Stehman & Czaplewski, 1998): where N represents classes in the confusion matrix, A jj indicate the items in the major diagonal in the confusion matrix and represent the number of samples correctly identified, and n is the total number of collected samples (Conese & Maselli, 1992;Congedo, 2021). Where the area of a specific LULC category is small, the overall accuracy is less significant (Santé et al., 2010). The user's and producer's accuracy are the statistics of the quality of the simulated maps in each LULC class and were calculated based on the Eqs. (4) and (5) (Olofsson et al., 2014;Stehman & Czaplewski, 1998): where i refers to the class and p is the proportion of the area mapped as class i which is referenced as class i.
Producer's accuracy for class j calculated from the area proportion of class j, which was mapped as class j. Figure 3 and Table 3 represent Estonia's LULC dynamics for 2000, 2011, and 2019. Between 2000 and 2019, about 24.16 km 2 was added to the artificial surfaces, placed on previous agricultural lands. During the later years, the addition of newly built areas has slowed down, so from 2000 to 2011, the rate of increase in artificial surfaces was 1.64%, reaching 765.00 km 2 , while from 2011 to 2019, this rate has decreased by 1.49% (776.59 km 2 ).

LULC dynamics and urban expansion from 2000 to 2019
(3) Despite the increase of artificial surface areas, the decrease in agricultural lands with strong human impact was noticeable between 2000 and 2011. About 12.58 km 2 of agricultural areas were converted to artificial surfaces, and agricultural lands with a decline rate of − 5.23% reached 12,852.92 km 2 in 2011. However, from 2011 to 2019, these lands experienced a slight increase (+ 0.12%), reaching 12,867.72 km 2 due to forest land conversion. From 2000 to 2019, forest lands experienced an increased rate of 1.81%. While in 2000, forests covered 62.22% of Estonia, it extended to 27,206.92 km 2 in 2019, covering 63.36% of the total land in Estonia. Furthermore, the 6.69% increase in wetlands was notable during the study period. While the wetlands area covered 1752.26 km 2 in 2000 (4.01% of total lands), it expanded by 125.63 km 2 in 2019, resulting from agricultural and forest land transitions. Overall, in 2019, wetlands formed 4.29% of land in Estonia.

Spectral-textural indices as predictors of LULC change
The relationships between the LULC changes and 19 spectral-textural indices are presented in Table 4 and reveal to what extent the spectral-textural indices represent the LULC changes. The results explored that the uniformity of distribution of color saturation and entropy on the pixel difference of color value (sat_asm_21 and val_dent_21, respectively) were the best positive indicators of LULC changes between 2000-2011 and 2011-2019. Additionally, entropy on the difference of color value and saturation of pixel pairs, spatial clusters of high normalized difference built-up index values, and normalized difference water index (val_dent_21, sat_dent_21, ndbi_gearys_21, ndwi) had positive relationships with LULC shifts in the period between 2000 and 2011. And normalized difference water index and correlation among normalized difference built-up index pixel pairs  Tables 3 and 4 illustrates that from 2000 to 2011, the most influential index in determining the LULC shifts was the angular second moment (an indicator of orderliness) for color saturation (sat_asm_21), which indicates the appearance of wetlands in landscapes. Interestingly, the highest shifts of wetlands (6.06% increase) happened in the same period, and this indicator correctly detects this transformation over time. Furthermore, the artificial surfaces increased by 12.58 km 2 from 2000 to 2011, experiencing a 1.64% rise. Likewise, detecting changes in these areas by spectral-textural indices revealed the clusters of built-up areas, roads, cultural features, and bare soil (indicated by variables ndbi_gearys_21, sat_dent_21, and val_dent_21, respectively), and had strong positive influences on LULC shifts' detection in this period. Besides, from 2011 to 2019, urban areas' determinants (val_dent_21 and ndbi_corr_21) had the highest positive effects on representing the LULC changes. It is in line with the 11.59 km 2 increase (+ 1.49%) of artificial surfaces in this period.
The decrease in agricultural lands was also a prominent direction of LULC changes between 2000 and 2011. Likewise, inverse difference moment of normalized difference built-up index and entropy of the blue band (ndbi_idm_21, blue_entropy_21) as determinants of vegetation clusters and diversity of landscape composition regarding vegetation edges represented this transition. Furthermore, the increase in water areas from 2000 to 2019 was also determined by the high importance of the normalized difference water index, representing this shift.
Projection of LULC demand and transition probability in 2030 using MCA The demand quantity was another primary input of the simulation model calculated by MCA and is presented in Table 5. The estimation of demands was based on the transition probability and area of the changes in LULC from 2011 to 2019. The results show that it is 78.30% probable that artificial surfaces remain the same, and its probability of transitions was mainly affected by agricultural lands by 13.03%. While there is a 94.87% likelihood of agricultural lands to remain unchanged, forests impact 4.15% on its transition probability. Also, the probability of the forest lands remaining constant was estimated at 96.90%, with 1.97% effects on agricultural lands. Furthermore, there is an 86.30% probability of wetlands remaining intact, while the likelihood of forest being turned into wetlands is estimated at 11.72%, and wetlands' transition probability of being water is estimated a 2.09%, which result in a 97.17% likelihood of water areas to remain unchanged in 2030.

ANN probability of occurrence results
This study explored the LULC changes using the hybrid ANN-CA-MCA model in the GIS environment. Taking advantage of a back propagation-ANN algorithm, the probability of occurrence was estimated to be one primary input in the simulation model to achieve more realistic results. Figure 4 depicts the ANN probability of occurrence (spatial location) in each LULC category (rows) as the result of competition with other LULC types and its final likelihood of changes applied for simulation of 2011, 2019, and 2030 (columns). The probability of occurrence, which is the result of competition between LULC categories, ranges between 0 and 1; the higher value (dark gray) represents the higher combined probability of occurrence for each LULC category, and the lower values (light white) show less probability of occurrence and conversions in the state of cells over time. Table 6 illustrates the relationships between artificial surfaces' probability of occurrence (dependent variable) and 19 spectral-textural indices (independent variables). The statistical results of the multiple regression analysis revealed the importance of spectral-textural indices on the expansion of artificial surfaces. The adjusted R 2 value was 0.94 for the probability of occurrence in 2011 and 0.93 in 2019 and 2030, which indicated that artificial surfaces' probability of occurrence is explained extremely high by spectral-textural indices. We explored that the uniformity of distribution of color saturation, entropy on the pixel difference of color value, saturation value, and average on the sum of color saturation pixel pairs (sat_asm_21, val_dent_21, sat_ dent_21, and sat_savg_21, respectively) were the best positive indicators of the probability of occurrence in 2011 and 2019. Also, correlation among normalized difference built-up index pixel pairs and NIR band mean (ndbi_corr_21 and nir_mean_21) most influenced the artificial surfaces' probability of occurrence in 2030. In contrast, homogeneity of normalized difference built-up index pixel pairs (ndbi_idm_21) was the most negative influential factor over these periods.

CA Simulation for LULC and urban expansion in 2011 and 2019
The CA model was performed to simulate the future projection of LULC. The simulation was run for LULC in 2011 and 2019. Then, the results were validated against the actual maps of those years to explore the accuracy of the predicted indices capable of LULC simulation in 2030. Figure 5 represents the implementation of the hybrid ANN-CA-MCA model in the simulation of LULC in 2011 and 2019 and the actual LULC. Also, the overlay was performed to visualize the spatial changes between simulated and LULC maps .
The largest spatial changes between the simulation and the actual LULC maps were found in the mining sites in the northeastern Estonia and the surroundings of Tallinn and Tartu; these changes may be related to the spectral responses of construction covers.

Accuracy assessment for the LULC Simulation in 2011 and 2019
Accuracy assessments utilizing the kappa coefficient and overall accuracy were calculated for the simulation of LULC in 2011 and 2019 and are  presented in Table 7. The kappa coefficient was 0.87, with an overall accuracy of 93.46% for the simulation of LULC in 2011, indicating the excellent simulation results and consistency in the simulation of the patterns. The accuracy assessment was performed for the LULC 2019 simulation to ensure the accuracy of the model simulation and parameter selection. Respectively, the kappa coefficient of 0.91 and overall accuracy of 95.38% significantly proved the applied model's accuracy. Besides, the kappa coefficient and overall accuracy for the simulation of 2019 were higher than for 2011. Producers' and users' accuracy for both simulations were higher than 90% in water areas, agriculture, and forest lands. At the same time, wetlands' and artificial surfaces' accuracy was over 70%, indicating the overall model parameters' high accuracy.

Simulation of LULC in 2030
Based on the overall conditions set for the hybrid ANN-CA-MCA model, we predicted the LULC changes in Estonia for 2030. Figure 6 and Table 8 represent the proposed model LULC simulation results. Using spectral-textural remote sensing indices, the prediction of changes in artificial surfaces was estimated to increase at a rate of 1.33% and expected to reach 787.04 km 2 in total, with a similar pattern to the previous periods. The most significant change was modeled for forest lands with a decrease of 34.57 km 2 areas. In comparison, the agricultural lands were expected to experience an increase by 14.90 km 2 and wetlands by 9.31 km 2 estimated in 2030. A closer inspection of Table 8 illustrates that from 2000 to 2030, a 34.62 km 2 growth of artificial surfaces (+ 4.40%) (Supplementary Fig. 6) and an increase of 134.94 km 2 of wetlands (+ 7.15%) are expected to be the most significant LULC changes in Estonia. Besides, the water class was the most stable class with a 6.89 km 2 increase in total, and the agricultural lands were the most dynamic class with conversion areas of 642.31 km 2 , respectively. Urban expansion in Tallinn, Tartu, and Pärnu in 2030 Here, the spatial details of urban expansion in 2030 can be examined at the local scale. In this regard, the artificial surfaces indicating the urban areas of 2019 and 2030 were overlayed in three major cities in Estonia, including Tallinn, Tartu, and Pärnu, represented in Fig. 7. While these cities will expand from the existing urban areas to exceed the cities' boundaries, the new development areas will primarily be distributed around the existing urban areas, indicating a prominence of infilling expansion. In Tallinn, where the water restricted expansion to the north and south, the expansion will be situated in the west and east directions. Expansion of artificial areas in the city of Tartu will be located mainly to the south and south-west. Indeed, scattering patterns of expansion in Pärnu as a summer capital will be placed in the northeast. A visual interpretation of urban expansion in these three cities also reveals that the vast majority of expansion was predicted in Tallinn, the capital city, with more population. As the second major city in Estonia, Tartu predicted more expansion in comparison to Pärnu.

Spectral-textural indices application to predict LULC changes
In the current study, we revisited the usage of spectral and GLCM-based textural indices, common for LULC classifications, in the context of change detection with the hybrid ANN-CA-MCA model. While utilization of the spectral-textural indices in LULC change detection at the landscape scale is still in infancy, other studies employed spectral-textural indices for LULC classifications based on satellite bands (Abubakar et al., 2021;Chaturvedi et al., 2020;Jin et al., 2018;Lu et al., 2014;Rodriguez-Galiano et al., 2012;Tassi et al., 2021), principal component analysis data (Ghasemi et al., 2021;Park & Guldmann, 2020), and selected spectral indices (Halmy et al., 2015;Naboureh et al., 2017). One interesting finding was that textural indices representing homogeneity and similarity (such as angular second moment for color bands) represent (semi) natural landscape areas, such as wetlands and forests. In contrast, diversity-related indices (e.g., difference entropy) indicate artificial surfaces and cities. This finding supports the idea that textural indices calculated based on the HSV space bands extend our opportunities of tracking landscape transitions by accounting for meaningful information about coloristic landscape attributes.
Another important finding was that the spectraltextural indices provided the structural landscape information at pixel-level, addressing the heterogeneity in space for our hybrid cell-based model. Our study is one of the first attempts to thoroughly examine the effectiveness of non-uniform spatial information derived from spectral-textural indices in pixel-level detecting changes in multiple LULC categories. In accordance with the present results, previous studies have demonstrated that spectral and textural indices are useful indicators for monitoring and understanding the heterogeneous LULC. Yatoo et al. (2022) used various pixel-based spectral indices to monitor built-up expansion with a CA-ANN model and found that spectral indices are appropriate in understanding the spatial extent of urban classes. He et al. (2011) proved that GLCM-textural features are effective in detecting LULC changes in rural-urban fringe areas, which are often heterogeneous and fragmental in nature. Other studies also developed different approaches to answer the heterogeneity in the CA model cell space, including the patch-based approach of landscape metrics (Fenta et al., 2017;Lin et al., 2020;Liu et al., 2020a, b;Yang et al., 2016a, b), applying different cell sizes and partitioned grid cells Xu et al., 2018Xu et al., , 2021, and specifying the different neighborhood characteristics (Feng & Tong, 2019;Mozaffaree Pour & Oja, 2021). Since the accuracy of our implemented model using the pixel-based spectral-textural data was similar to those studies mentioned, this finding suggests the spectral-textural data usefulness for other modeling approaches.
The hybrid ANN-CA-MCA model application to predict LULC changes and urban expansion We hypothesized that the spectral-textural properties of landscapes provide a sufficient proxy for LULC transitions. Applying the spectral-textural indices by the hybrid ANN-CA-MCA model determined the accuracy of predictions over 90%, indicating the   Huang et al. (2018); Li et al., (2017a, b); Liang et al., (2018a, b); Liu et al., (2017a, b); and Zhao et al. (2020). Liu et al., (2017a, b) indicated the excellent accuracy of the CA model implemented with GeoSOS-FLUS compared to the other platforms. While they used the soil, climate, and socioeconomic data as drivers of LULC changes, the overall accuracy of their simulation was 75%, and the Kappa coefficient reached 0.67. Furthermore, Liang et al., (2018a, b) investigated the planning driving effects on the LULC dynamics and explored the higher simulation accuracy (approximately 5%) that could be achieved when the reality of policy implementation came to the model. Our study confirms that the ANN probability of occurrence algorithm is an effective pixel-based tool to discover and mine the transition rules based on the spectral-textural predictors for multiple LULC types. Besides, adding ANN to the CA model allows the dynamic land transformation and competition in the cell level using the spectral-textural indices information. While the probability of occurrence is distributed between multiple LULC types, its' self-adaptive inertia addresses the randomness and uncertainty in the simulation of LULC changes (Liang et al., 2018a, b). The high value of adjusted R-square also detected the effectiveness of this algorithm for estimation of the probability of occurrence; indeed, it reflected the diversity and complexity of LULC transitions over time. This finding corroborates the ideas of Yang et al., (2016a, b), who suggested that utilizing indicators characterizing the LULC patterns and composition, such as landscape metrics and physiognomy, with the ANN model can generate more accurate local transition rules and increase the quantity of simulation. While the scale dependency and calculation methods of landscape metrics are challenging (Yang et al., 2016a, b), the capability of ANN in capturing the transition rules from heterogeneous LULC is acceptable. In contrast to the mechanism of logistic regression, which is a wellknown linear-fitting mining algorithm    Mustafa et al. (2017), the ANN algorithm can discover the spatial probability of changes based on a competition mechanism and proper for the non-linear systems like LULC changes. Other studies applied physical, socioeconomic, and proximity drivers for LULC changes simulation based on integrated models, including the ANN-CA-MCA (Girma et al., 2022;Sankarrao et al., 2021), ANN-CA (Gharaibeh et al., 2020;Hu et al., 2018;Xu et al., 2021;Yang et al., 2016a, b), and CA-MCA (Guan et al., 2019;Mondal et al., 2016;Rimal et al., 2017;Singh et al., 2015). While covering all aspects and predictors of the LULC transitions requires sufficient knowledge about the study area and trends of changes, adopting spectral-textural indices helps to overcome the selection limitations and reliability in simulation results. Therefore, this research was based on observable remote sensing landscape physiognomy indices. While the simulation accuracy was similar to those studies, utilizing spectral-textural indices straightly into a hybrid model has not previously been explored.

LULC and urban expansion simulation
While construction of new settlements in scattering form is a characteristic attribute of urban expansion in main cities of Estonia after 2000, this study confirmed that the scattering patterns of new constructions are expected to continue as the infilling development form fed by the existing infrastructures from 2019 to 2030. The same pattern was explored by Mozaffaree Pour and Oja (2021) in the simulation of Tallinn and its suburbs, indicating the infilling patterns of new developments in 2030. Besides, new constructions will be located near the existing expansions of Tallinn, Tartu, and Pärnu. This finding is in line with the spatial plan of Estonia 2030 + , which supports the living and economic environments of the existing settlement to prevent the scattering.
Estonia is rich in natural forest lands. The share of forest areas from the total LULC was 62.22% in 2000. The results showed a decrease in forest lands by − 0.16% from 2011 to 2019, and the hybrid ANN-CA-MCA model results confirmed a continued decline (− 0.12%) from 2019 to 2030. This result reflects those of Oja (2020), who also found that the structure of the forests in Estonia is changing due to the overgrown agricultural lands, increase in mixed forests, and decrease in coniferous forests. While forest shifts significantly impact animal habitats (Lee & Chang, 2011), climate (Cegielska et al., 2018), and soil attributes , effective longterm reforestation regulations must be established to preserve the forest areas. Correspondingly, as proposed in the Estonian Environmental Strategy 2030, the planning should provide a balance satisfaction of socioeconomic, cultural, and ecological needs relating to the utilization of forests (Estonian Ministry of the Environment, 2007).
Agricultural lands are the second-largest type of LULC in Estonia (30.93% of LULC in 2000). Estonia experienced a significant loss of agricultural lands (− 5.23%) from 2000 to 2011. Poska et al. (2014) indicated that the decrease in arable farming lands resulted from land ownership and agricultural policy changes after 1991 in Estonia. Conversely, a slight increase of 14.80 km 2 was observed from 2011 to 2019, and our hybrid model predicted adding about 14.90 km 2 (+ 0.12%) to agricultural lands from 2019 to 2030. This slight increase is essential, especially in Estonia, where cultivated land is the primary source of feeding the population, quite consistently forming 2/3 of the agricultural land (Oja, 2020). However, considering this slight rise, the increase in the future share of cultivated lands in active use from 2019 to 2030 is an optimistic perspective toward sustainability and a primary goal of sustainable Estonia (Republic of Estonia; Government office, 2021).
Wetlands covered 4.01% of land in Estonia in 2000 and experienced a drastic increase (+ 6.69%) in 2019. From 2000 to 2019, the changes in wetlands resulted from 1.27% loss of forest lands, − 2.47% of water areas, and 0.39% conversion of agricultural lands. During the past decades, wetlands, as common features in the Estonian landscape (Kimmel et al., 2010), were wisely protected. It was predicted that 9.31 km 2 will be added to wetlands by 2030. This pattern will continue with the gain of wetlands in 2030. Moreover, international wetland conservation projects such as Ramsar (Blumenfeld et al., 2009) and Natura 2000 Network (Kimmel et al., 2010) motivated preserving wetlands. Still, the probable loss of living biomass and organic matter needs proper mitigations.
The results of the hybrid ANN-CA-MCA model also predicted a 0.02% loss of water areas from 2019 to 2030. While water experienced fewer changes than the other LULC types, from 2000 to 2019, the inland water increased by 6.97 km 2 . A minor loss is expected regarding the transitions that ANN and MCA estimated in the model. Interestingly, our model implementation results support a realistic explanation of the changes in water areas in Estonia. During the past few years, respectively, Estonia started the process of dam removal, restoring the fish pathways (Tamm, 2018) and decreasing the artificial lakes, which the significant result would be a decrease in the areas of water.

Conclusions
This study investigated the Estonian LULC changes over the past two decades and, for the first time, predicted Estonian LULC in 2030. Here, MCA and ANN were applied to enhance the mining of transition rules by the CA model utilizing the spectral-textural indices. The results of this study indicate the following: 1. Spectral-textural indices powerfully explained the LULC changes and their probabilities of occurrence. In our work, 19 spectral-textural indices explained over 70% of the variability in the LULCs and more than 90% of the probability of occurrence in artificial surfaces. 2. The changes in LULC from 2019 to 2030 are similar to the changes from 2011 to 2019. The predicted change in the area of artificial surfaces was an increased rate of 1.33% to reach 787.04 km 2 in total. The most significant change was modeled to be a decrease of 34.57 km 2 of forest lands in 2030. Agricultural lands and wetlands were modeled to increase by 14.90 km 2 and 9.31 km 2 , correspondingly. 3. Modeling results showed that Tallinn, Tartu, and Pärnu cities can expand from the existing urban areas to exceed the city's boundaries in 2030; the new development areas will mainly be located near the existing expansion indicating a prominence of infilling expansion.
These findings suggest several courses of action for long-term spatial planning in Estonia. Some initial procedures include maintaining the importance of living and economic environments of the existing settlements to prevent the scattering of new ones, preservation of Estonia's natural landscapes, a policy priority for the stable care of forest lands to maintain biodiversity, encouraging environmentally clean production by the motivation of people to be involved in the conventional agricultural sector, and organizing the cultural and nature tourism in the way of learning sustainability.
Further research is required to fully understand the cellular automata model's scale effects and neighborhood interactions by implications of different cell sizes of landscape physiognomy indices. Besides, more research could be conducted to characterize the spatial planning policies to explore the reflection of macro-level effects on the pixel-level morphological predictors of LULC changes.

Funding This research was funded by Estonian Research
Council grant number PRG352. Funding from the Academy of Finland (PEATSPEC, decision no 341963) is acknowledged by Iuliia Burdun. The authors are thankful for the thorough comments of anonymous reviewers on the earlier draft of this paper. The authors also thank Dr. Liang Xun, developer of the GeoSOS-FLUS software, which guided the implementation of the model process.

Data availability
The authors confirm that the data supporting the findings of this study are available within the article and/or its supplementary materials.

Declarations
Competing interests The authors declare no competing interests.