**Research on Construction Land Use Benefit and the Coupling Coordination Relationship Based on a Three-Dimensional Frame Model—A Case Study in the Lanzhou-Xining Urban Agglomeration**

**Weiping Zhang 1, Peiji Shi 1,2,\* and Huali Tong <sup>1</sup>**


**Abstract:** Coordinating the social, economic, and eco-environmental benefits of construction land use has become the key to the high-quality development of Lanzhou-Xining urban agglomerations (LXUA). Therefore, based on the coupling coordination connotation and interaction mechanism of construction land use benefit (CLUB), we measured the CLUB level and the coupling coordination degree (CCD) between its principal elements in LXUA from 2005 to 2018. Results showed that: (1) The construction land development intensity (CLDI) in the LXUA is generally low, and spatially presents a dual-core structure with Lanzhou and Xining urban areas as the core. (2) The comprehensive construction land use benefit has increased over time, but the overall level is not high. The spatial differentiation is obvious, and the core cities (Lanzhou and Xining) are significantly higher than other cities. (3) The regional differences in the subsystem benefit of construction land use are obvious. The social benefit and economic benefit showed a "convex" shape distribution pattern of "high in the middle and low in the east and west wings", and regional differences of economic benefit vary greatly. The eco-environmental benefit was relatively high, showed a "concave" shape evolution in the east–west direction. (4) In addition, the CCD of the CLUB were still at a medium–low level. The higher the administrative level of the city, the better the economic foundation, and the higher or better the CCD of the social, economic, and eco-environmental benefits. (5) The CCD is inseparable from the influence of the three benefits of construction land use. Therefore, different regions should form their own targeted development paths to promote the coordinated and orderly development of LXUA.

**Keywords:** construction land development intensity; construction land use benefit; coupling and coordination relationship; spatiotemporal evolution; Lanzhou-Xining urban agglomeration

#### **1. Introduction**

Land is the basic carrier of all social and economic activities, and can provide space and resources for human production and life [1]. As an important part of land resources, construction land mainly includes urban construction land, rural residential land, transportation land, water conservancy land, and other construction land. Its utilization is not only related to the development of the secondary and tertiary industries of the national economy, but also has an important impact on the coordinated development of urban and rural areas [2,3]. Land use benefit refers to the comprehensive output of social, economic, and eco-environmental benefits obtained by human capital, labor, and technical input, which is related to the sustainable development of a region [4].

At present, due to rapid global urbanization, most developed countries and regions are striving to complete the simultaneous optimization of land use in society, economy, and

**Citation:** Zhang, W.; Shi, P.; Tong, H. Research on Construction Land Use Benefit and the Coupling Coordination Relationship Based on a Three-Dimensional Frame Model—A Case Study in the Lanzhou-Xining Urban Agglomeration. *Land* **2022**, *11*, 460. https://doi.org/10.3390/ land11040460

Academic Editors: Víctor Hugo González-Jaramillo and Antonio Novelli

Received: 18 February 2022 Accepted: 22 March 2022 Published: 24 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

ecology, and have made some achievements [5–7]. As the main target of urban expansion in the coming decades [8], driven by industrialization, new urbanization and modernization, the urban population of developing counties has increased rapidly, and the scale of urban construction land development has been expanding continuously. Especially in China, the urbanization rate has increased from 17.90% in 1978 to 59.58% in 2018 [9]. However, with the rapid increase in urbanization level, the contradiction between man and land has become increasingly acute, and problems such as extensive land development and utilization, disorderly expansion, and land pollution have frequently occurred, which has led to an imbalance between social, economic, and eco-environmental benefits of land use, especially in urban agglomeration areas where social and economic activities are more prevalent [10,11]. In addition, urban sprawl will inevitably affect the implementation of cultivated land protection policies and even national food security [12]. Therefore, it is of great significance to study the development intensity and multi-dimensional use benefits of construction land for the rational and optimal allocation of land resources and the sustainable development of urban and rural areas.

Scholars have done a lot of research on land use benefit, mainly focusing on the evaluation of land use benefits, the diversification and applicability of research methods and models, and the excavation of influence mechanisms. In the early days, people only paid attention to the economic benefits of land. For example, the law of land rent in western economics laid a theoretical foundation for the study of the economic benefits of land use [13]. Fulton et al. also verified that urban land expansion was closely related to urban land economic benefits [14]. However, urban sprawl also leads to the increase of social and ecological costs such as prolonged commute time, waste of resources, and destruction of ecological environment pollution [15]. Therefore, while discussing the principle of maximum return brought by non-renewability and restriction of land in principles of land economics, scholars emphasize that land use should meet social goals such as wealth production, fair distribution, and ecological protection; and improve the overall benefit of land use by means of political, legal and economic leverage [16]. With the challenge of sustainable development, people are more aware that the goal of land use is the sustainable use of land resources, which must take into consideration the coordinated development of social, economic and eco-environmental benefits [17,18]. According to previous research, the research object and perspective of land use benefit is not only to evaluate the land use benefit, but also to change the coordination relationship between land use benefit and urbanization and other factors. For example, Zhang et al. analyzed the coupling and coordination between urban land use efficiency and urbanization in 34 prefecturelevel cities in the three northeastern provinces, and found that although there is a mutual response relationship between them, urbanization pays too much attention to development speed and despises development quality, resulting in low overall development level and low coupling and coordination degree [19]. He et al. constructed the theoretical framework of land use benefit and industrial structure evolution, and found that industrial structure evolution has an obvious single effect on land use benefit [20].

At the same time, there are regional differences in the study of land use benefits. In developed regions such as America and Europe studies of combined land use, urban expansion, and ecological environment management to try to analyze and evaluate ecosystem services and economic benefits through land use modeling, land protection, and planning behavior [21–23]. While in developing countries, the research is more related to urbanization and urban sprawl, focusing on the rational use of urban land, coordination of urban and rural land, and management strategies for the process of rapid urbanization [24,25]. Research methods and models of land use benefits are also hot topics that scholars pay attention to. The existing research is mainly based on the construction of a land use benefits index system, with the data envelopment analysis model, neural network model, multiobjective linear programming model, SWAT model, coupling degree model, and other methods [26–30]. With the development of science and technology, GIS, remote sensing and spatial measurement methods have been introduced into the study of land use and

land expansion, and the spatial monitoring and analysis functions of these methods have been applied to realize the visual expression of the research results [31,32]. Of course, due to data availability and the research objects being different, the weighting methods and indicators selected are also different. Methods such as the analytic hierarchy process (AHP), coefficient of variation method, expert scoring method, and entropy weight method are widely used [33–35]. Based on this, Ran et al. used the Friedman chi-square test, Spearman correlation, consistency test, and one-way analysis of variance to compare different subjective and objective weighting methods, they found that the comprehensive index method, rank-sum ratio method, entropy method, and the integrated entropy methods all have significant statistical characteristics, and each has its own advantages and disadvantages. A more scientific, comprehensive approach is required [36]. Research on the influencing factors of land use benefits generally focuses on economic, social, transportation, and political factors [37–40]. Some studies have also found that the determinants are related to urbanization and industrialization, accessibility, and economic transformation [41,42]. These studies provide policy and guidance for improving land use efficiency and promoting sustainable development of cities.

From the current research progress, we also find that most of the land use benefits studies focus on the desirable outputs of limited land use, such as economic benefit, while there are few studies on undesirable outputs such as pollution and industrial emissions from land use. In addition, the study areas are mainly concentrated in areas with a high economic level, while the research on land use benefits of the northwest region China, where the ecological environment is relatively fragile, is relatively weak. Therefore, this paper selects the LXUA which is the important urban agglomeration on the Silk Road Economic Belt as an example. On the basis of elaborating and analyzing the connotation and mechanism of the coupling coordination of multi-functional benefits of construction land use, a threedimensional framework system of CLUB is constructed including society, economy, and eco-environment. Secondly, using an entropy method, a composite index model and coupling model, and the application of CLUB is demonstrated. Finally, according to the evolution law and spatial differentiation characteristics of the coupling coordination relationship of CLUB in different regions, an optimization strategy for sustainable land use development is proposed according to local conditions. The overall workflow of the study is shown in Figure 1.

**Figure 1.** The coupling interaction relationship of the CLUB system and workflow of the study.

The structure of this article is as follows. Section 2 explains the conceptual framework of the interactive coupling relationship of the CLUB system. Section 3 introduces the current situation of LXUA, data sources, and the research method. Section 4 analyzes the spatiotemporal differentiation of the CLUB based on the CLDI, and explores the coupling and coordination relationship between the social, economic, and eco-environmental benefits of the construction land use. Section 5 discusses and analyzes the important conclusions. Section 6 summarizes the main conclusions and future research directions.

#### **2. Interaction Coupling Mechanism of CLUB**

The economic, social, and eco-environmental benefits of construction land use are closely related, restricted, and promoted by each other [43]. On one hand, land is an important carrier of all human activities. With the rapid advancement of urbanization, construction land has been continuously developed and utilized, providing space and resources for various activities. Its quantity and quality are closely related to the social and economic benefits. If people develop land and do not exceed the resource's environmental carrying capacity, we can get economic benefit continuously. If the economic benefit is invested in improving people's livelihood and infrastructure construction, good social benefit will be produced. The improvement of social benefit will improve regional production conditions, promote further development of the regional eco-environment. On the other hand, in the process of land use, due to the limitation of technology, capital, and knowledge level, people will have negative economic scale effects, which can lead to environmental pollution, ecological systems destruction, and affect the eco-environmental benefit of land. The destruction of the eco-environment will worsen the local environmental conditions and affect the social benefit of land use. Moreover, it will also reduce the local resources environmental carrying capacity, resulting in the loss of land economic benefit. To sum up, only when the three achieve dynamic coordination and balance, can they promote the effective improvement of the whole system benefits to a greater extent and maximize the benefits (Figure 1).

#### **3. Materials and Methods**

#### *3.1. Study Area*

The LXUA is a 9.75 × 104 km regional urban agglomeration in the upper reaches of the Yellow River in the inland northwest China (Figure 2), including nine cities, Lanzhou, Baiyin, Dingxi, Linxia, Xining, Haidong, Hainan, Haibei, and Huangnan, and a total of 39 counties. The population concentration degree is relatively high, and it is an important radiation center and growth pole in Northwest China. As of 2018, the GDP of the LXUA reached CNY 515.59 billion, accounting for 51.13% of the total of Gansu and Qinghai province. As the hub of the Silk Road Economic Belt and the Yangtze River Economic Belt, LXUA's geographical advantages have become increasingly prominent, especially the construction of the new land–sea passage in the west, which makes LXUA's hub position of the "Sixth Ring in the Middle" more prominent. This region has good resource endowment and belongs to a region with good soil–water combination conditions in Northwest China. However, the economic level within the urban agglomeration varies greatly, and the imbalance of regional development is outstanding. Therefore, under the background of the new round of western development, how to balance the benefit relationship between land use, promote the coordination between development, utilization and protection of land, and realize the promotion of comprehensive benefit of land use in the region, has become the primary task of high-quality development of urban agglomerations.

**Figure 2.** Study area.

#### *3.2. Data Sources*

This research takes 39 counties in LXUA as the basic research unit. The vector boundaries of urban agglomerations and the construction land data are all from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn accessed on 15 November 2021). According to the national countylevel land use classification system, the data of the corresponding year was revised to extract and divide the construction land area data of each county. The socioeconomic data of each research unit was mainly derived from the Statistical Yearbooks of the Chinese County, Statistical Yearbooks of Gansu and Qinghai Province. The eco-environmental data was mainly derived from the environmental protection bureaus of cities and prefectures and field survey.

#### *3.3. Methods*

#### 3.3.1. CLDI Measurement

Development intensity can be used to measure the extent of land development and utilization by human activities in a certain region. As for the development intensity, 20% is generally regarded as the standard line of livability and 30% as the warning line, globally. The National Land Planning Outline (2016–2030) proposes that the development intensity of China's land should not exceed 4.62%. This paper selects the proportion of construction land area to the regional total area proposed by Fan in the main functional zone as a measure of CLDI [44]. The equation is as follows:

$$\text{CLDI} = \frac{\text{CLA}}{\text{TA}} \times 100\% \tag{1}$$

In Equation (1), CLDI represents construction land development intensity, *CLA* is construction land area, and *TA* is the total area of region.

#### 3.3.2. CLUB Index System

The CLUB mainly refers to the comprehensive output of social, economic, and ecoenvironmental benefits caused by the labor, capital, and technology invested by human beings in construction land. Therefore, drawing on previous studies [43,45,46], and combined with the connotation of CLUB and the scientific principle of index selection, and the availability of data, this paper constructed a three-dimensional evaluation index system including social, economic, and eco-environmental benefits (Table 1).

**Table 1.** Three-dimensional framework evaluation index system of CLUB.


3.3.3. CLUB Index Weight Setting and Score Calculation Data Preprocessing

Different indicators have positive and negative effects, in order to ensure the rationality of the evaluation results, it was necessary to standardize the data to ensure the uniformity of the dimensions of each indicator by the following:

$$\mathbf{x}'\_{ij} = \begin{cases} \begin{array}{c} \frac{\mathbf{x}\_{ij} - \mathbf{x}\_{\text{min}}}{\mathbf{x}\_{\text{max}} - \mathbf{x}\_{\text{min}}} \text{Positive indicator} \\ \frac{\mathbf{x}\_{\text{max}} - \mathbf{x}\_{ij}}{\mathbf{x}\_{\text{max}} - \mathbf{x}\_{\text{min}}} \text{Negative indicator} \end{array} \tag{2}$$

In Equation (2), *x ij* represents the standardized value, *xij* is the original data value of the *j*-th indicator, *x*max and *x*min are the maximum and minimum values of the *j*-th indicator.

#### Weight Calculation

Due to the different degrees of dispersion of different indicators, the impact on the comprehensive evaluation is also different. Therefore, this study used the more objective entropy method to calculate the index weight [46,47]. The weight results are shown in Table 1.

Comprehensive Value and Subsystem Score Calculation

According to the above standardized values and weights, the social benefit value (*Usoc*), economic benefit value (*Ueco*), and eco-environmental benefit value (*Uenv*) of construction land use were calculated in combination with the composite index method [46]. Then, the comprehensive value (*Ucom*) of construction land use of the *i*-th sample was obtained by using the weighted evaluation method of the criterion layer.

$$
\Omega I\_{\text{soc}/\text{eco}/\text{encv}} = \sum\_{j=1}^{n} (w\_{\circ} \times \mathbf{x}'\_{\circ j}) \tag{3}
$$

$$
\hbar L\_{\rm com} = \varkappa \mathcal{U} L\_{\rm soc} + \beta \mathcal{U} L\_{\rm cro} + \gamma \mathcal{U} L\_{\rm env} \tag{4}
$$

In Equations (3) and (4), *Usoc*, *Ueco*, and *Uenv* represent the social benefit, economic benefit, and eco-environmental benefit of construction land use, respectively; *x ij* is the standard value of the *j*-th index in three systems; *wj* represents the weight of the corresponding index. *Ucom* represents the comprehensive benefit of construction land use. *α*, *β*, and *γ* are undetermined coefficients, which represent the importance of the social, economic, and eco-environmental benefits of construction land use, respectively. Because the three need to coordinate and advance together and are almost equally important, this paper takes *α* = *β* = *γ* = 1/3.

#### 3.3.4. Coupling Degree (CD) and Coupling Coordination Degree (CCD)

Coupling is usually used to describe the degree to which two or more systems interact and affect each other [48,49]. This paper constructs a CD model among the threedimensional benefits of social, economic, and eco-environmental. Based on the CD, this paper further draws on the CCD to measure the degree of harmony and consistency of the social, economic, and eco-environmental benefits of construction land use in the process of change [50]. The equations are as follows:

$$\mathbf{C} = \left[ \frac{\mathbf{U}\_{\rm cco} \times \mathbf{U}\_{\rm soc} \times \mathbf{U}\_{\rm cnv}}{\left( \frac{\mathbf{U}\_{\rm cco} + \mathbf{U}\_{\rm cco} + \mathbf{U}\_{\rm cun}}{3} \right)^{\frac{1}{3}}} \right]^{\frac{1}{3}} \tag{5}$$

$$D = \sqrt{\mathbb{C} \times \mathcal{U}\_{\text{conv}}} \tag{6}$$

In Equations (5) and (6), *C* denotes the coupling degree of social, economic, and ecoenvironmental benefits of construction land use, and *C* ∈ [0, 1]. *D* denotes the coupling coordination degree between social, economic, and eco-environmental benefits of construction land use, and *D* ∈ [0, 1]. The larger the *D*, the higher the CCD. When *C =* 0, the CD is extremely low, and the systems, or between elements within the system, are in a disordered state. When *C* = 1, the CD is the largest, and benign coordinated coupling is achieved between systems or between internal elements of the system.

On the basis of existing research, we found that many scholars have used the critical value method to classify CD and CCD. However, this paper does not subjectively carry out quantitative division, but chooses the objective quartile method. The CD and CCD are divided into four stages from low to high level: lower level, medium level, higher level, and benign (optimal) level.

#### **4. Results**

#### *4.1. Spatiotemporal Change of CLDI*

According to Equation (1), it was found that the CLDI in LXUA had an upward trend from 2005 to 2018. The average of CLDI in 2005, 2010, and 2018 were 6.53%, 7.31% and 8.54%, respectively. In 2005, CLDI was dominated by low-value and lower-value counties, accounting for 38.46% and 23.08%, respectively. In 2010, the CLDI showed a steady increase, compared with 2005, the proportion of low-value areas decreased, and the number of medium and high-value areas increased. By 2018, the CLDI showed an accelerating trend, the number of counties with a development intensity of 2% and above continued to increase. However, compared with cities in eastern and central China [31], the CLDI in LXUA is relatively low.

As shown in Figure 3, the CLDI differs on different spatial scales. On the provincial scale, the CLDI in Gansu area (8.63%) is higher than that in the Qinghai area (6.95%). On the prefecture-level city scale, Xining, Lanzhou, and Baiyin ranked the top three in terms of development intensity. At the county level, the CLDI in the municipal districts of Xining, Linxia, Lanzhou, and Baiyin all exceeded the national land development intensity of 4.60% in the National Land Planning Outline (2016–2030). In terms of changes in the entire urban agglomeration, the CLDI in the LXUA generally presents a "core-periphery" spatial structure with the Lanzhou and Xining urban areas as the core, and other peripheral areas

decreasing in turn. The closer it is to the city center, the stronger the gathering capacity of construction land and the larger the development scale. On the contrary, the further from the city edge, the CLDI gradually decreases.

**Figure 3.** The spatial pattern of CLDI in LXUA.

#### *4.2. Spatiotemporal Change of CLUB*

As shown in Figure 4a, the social benefit of construction land use showed a trend of "first rise and then decline", and first rose from 0.28 in 2005 to 0.35 in 2010, and then declined to 0.33 in 2018. From Figure 4b, the variation coefficient of social benefit fell from 0.57 to 0.44, a decrease of 22.81%, which indicates that the regional differences of social benefit tend to narrow. From Figure 5a, the social benefit is characterized by a "convex" shaped spatial distribution. The social benefit of Lanzhou and Xining was always higher than that of other regions. Among them, the middle and high value regions are gradually concentrated in the central area connected with Lanzhou and Xining, while the social benefit in the peripheral counties of the urban agglomeration is low.

**Figure 4.** CLUB value and its variation coefficient. (**a**) Construction land use benefit level, (**b**) Variation coefficient of construction land use benefit.

As shown in Figure 4a, the economic benefit of construction land use from 2005 to 2018 showed a "continuous increase", and increased from 0.12 to 0.17, an increase of 41.67%. However, as can be seen in Figure 4b, the variation coefficient of economic benefit decreased from 1.01 to 0.93, with a small reduction of only 7.92%, which means that the regional differences of economic benefit were still large. From Figure 5b, the economic benefit presented a spatial pattern of "high in the west and low in the east, high in the middle and low on the outside". The middle and high value areas were concentrated in Huangzhong, Chengguan, and Ledu, and the low value areas were mainly distributed in the eastern fringe counties of the urban agglomeration.

**Figure 5.** Spatial pattern of the CLUB in LXUA. (**a**) Social benefit, (**b**) Economic benefit, (**c**) Eco-environmental benefit, (**d**) Comprehensive benefit.

From Figure 4a, the eco-environmental benefit of construction land use showed a characteristic of "first decline and then increase". Itfirst declined from 0.61 in 2005 to 0.57 in 2010, and then increased to 0.73 in 2018. From Figure 4b, the regional differences of the eco-environmental benefit were smaller and tended to shrink slightly, with a shrinking rate of 10.71%. As can be seen in Figure 5c, the eco-environmental benefits were characterized by a "concave" shaped spatial distribution. In 2005, the high value areas were mainly concentrated in the western region. By 2018, the eco-environmental benefit of the eastern of the urban agglomeration had risen significantly, and showed a trend of "retreat from west to east".

As shown in Figure 4a, the comprehensive benefit was a continuous upward trend during the study period, and increased from 0.27 to 0.32, an increase of 18.52%. However, we can see from Figure 4b the variation coefficient of comprehensive benefit decreased from 0.38 to 0.37, with a decrease of only 2.63%, which indicates that the regional differences of comprehensive benefit demonstrated only a slight shrinking trend. From Figure 5d, the comprehensive benefit showed a distribution of "high in the middle and low in the outside", and the urban land use benefit with Lanzhou and Xining as the core were always higher than that of other counties. At the same time, the number of low-level counties

decreased significantly, from 20.52% to 7.69%. Overall, the CLUB of LXUA are mainly at low and lower level, and showed a sequential shift from low level area to lower level area.

#### *4.3. Coupling and Coordination Analysis of CLUB*

During the important period of socio-economic transformation and development of LXUA from 2005 to 2018, its land development, utilization structure, and output benefits were affected by population agglomeration, industrial development, and urban policies. As shown in Figure 6a, the coupling degree (CD) value C of CLUB generally fluctuated and rose, and the western counties were higher than the eastern counties. Higher and benign level counties of the CD gradually shifted from a scattered to centralized distribution, mainly in the western and central areas of the urban agglomeration. As the "back garden" of the central city, the demand for industrial development, infrastructure, and residential space in such areas may be in a stage of continuous increase, coupled with the strong promotion of the policy of linking the increase and decrease of construction land, which made the social, economic, and comprehensive benefits of construction land at a good level. Therefore, the CD is also high. Lower level counties of the CD were concentrated in the eastern and northern fringes of LXUA. Such areas were mostly development areas of traditional industry with slow economic development and a single structure, with primary industry accounting for a considerable proportion. Its low level of urbanization and industrialization, limited population agglomeration capacity, and slow expansion of built-up areas, resulted in social, economic, and comprehensive benefits of construction land use being low. Therefore, the CD is also lower.

**Figure 6.** Spatial differentiation of the CD and CCD of social, economic, and eco-environmental benefits of construction land use in LXUA. (**a**) Coupling degree, (**b**) Coupling coordination degree.

To judge whether the input and output of regional construction land development is reasonable and orderly, in addition to considering its interaction and correlation, it is also necessary to explore the coupling coordination degree (CCD) of CLUB. Figure 6b showed that the spatial pattern of the CCD was basically consistent with the CD, but the CCD was generally lower than the CD. In Yongdeng, Gaolan, Jingtai, Yuzhong, Anding, Longxi, Haiyan, and Huangyuan counties, where the CD was at a medium-low level, their coupling characteristics showed a low-level orderly coordination state, and the CLUB in such areas was at a lower level of coordinated evolution, and the industrial development was mostly in the primary stage, so the response sensitivity to the improvement of construction land

use efficiency was slow. Correspondingly, the CCD of higher CD regions of construction land use benefit was also higher, which reflects a high-level synchronous evolution state. From the spatiotemporal evolution, the CCD showed a continuous upward trend, and a regional distribution with optimal level coordination stages did not change significantly from 2005 to 2018, and was mainly concentrated in Lanzhou and Xining urban areas. The medium and lower level counties accounted for 41.03% in 2018. Overall, the CCD between the social, economic, and eco-environmental benefits of construction land use in LXUA is still in the lower level coordination stage, and the CCD needs to be improved urgently.

#### *4.4. Spatial Autocorrelation Analysis*

Spatial autocorrelation can well express the spatial relationship of CCD between the CLUB, so as to further reveal the spatial connections and differences among the research units [50]. Figure 7 shows the Moran's I of the CCD from 2005 to 2018. The Moran's I values are all positive, and *p* < 0.01, which indicates that the CCD had a significant positive spatial autocorrelation. The Moran's I showed a fluctuating trend. From 2005 to 2010, the Moran's I had a decreasing trend, indicating that the spatial agglomeration of the CCD was weakened. From 2010 to 2018, the Moran's I had a slowly increasing trend, indicating that the spatial agglomeration of CCD was slowly increasing.

**Figure 7.** Spatial agglomeration of the CCD of CLUB from 2005 to 2018.

In addition, local Moran's I was used to reveal the spatial association type of CCD, and was visualized in ArcGIS 10.6 (Figure 8). In 2005, 2010, and 2018, the spatial distribution of the high–high cluster was basically the same, mainly in the Xining urban area and some surrounding counties. Low–low clusters were mainly distributed in Longxi, Weiyuan, Tongren, and Xunhua. The spatial distribution of the two categories of spatial outliers that were statistically significant (high–low and low–high outliers) was fragmented.

**Figure 8.** Spatial correlation type diagram of CCD of CLUB in the LXUA.

#### **5. Discussion**

#### *5.1. Reasons for the Change of CLUB and CCD Relationship*

The evaluation of CLUB is conducive to promoting the intensive and efficient use of land resources and the sustainable development of the region, especially for the northwest region where the ecological environment is relatively fragile. As an important ecological barrier support area in western China, LXUA plays a key role in promoting the high-quality development of the Yellow River Basin, promoting the development of the western region, and ensuring the ecological security of Northwest China [51,52]. The research results show that the CLUB of LXUA is generally low, the polarization of difference is obvious, and the regional development imbalance is a prominent problem. This is consistent with the findings of Shi et al. [53]. Especially, there are significant regional differences in the social and economic benefits. The lack of long-term spatial control measures made Lanzhou and Xining the dominant places. The excessive agglomeration of low-level urban functions, and the limited radiation effect of polarized regions restrict the urban agglomerations integration as well as balanced economic development. In addition, the CCD among the social, economic, and eco-environmental benefit of construction land use increased from 0.48 to 0.55, and is still at a medium and low level coordination stage, which is also consistent with previous research [54–56]. However, the northwest region is currently in a stage of accelerated urbanization. The rapid population growth in some areas has increased the development of land resources. At the same time, the economic structure of the counties is similar, the industrial level is relatively low, and there are still many traditional industries with high input, consumption, and emissions. Coupled with the fragile ecological background and increasing resources, environmental pressure caused most counties to be assessed as lower level regions. This unsustainable development mode needs to be alleviated urgently.

#### *5.2. Suggestions for Promoting the CCD of CLUB*

With the promotion of ecological civilization construction and the transformation of new urbanization development, combined with the research results of this paper, it has become an inevitable choice for the future land use of LXUA to promote the coordinated development of the social, economic, and eco-environmental benefits. In this process, different types of regions should form their own targeted development paths. For lower-level coordination areas, they should first transform and upgrade the industrial structure to achieve the mutual assistance of urban functions, the industrial dislocation layout, and industrial chain linkage close development, to enhance the economic benefit of land and improve its social and eco-environmental benefits, promote the harmonious and orderly development of construction land in society, economy, and ecology. For medium coordinated regions they should rely on local characteristic industries and regional advantages to first achieve healthy economic development, thereby increasing investment in the construction of basic, urban public services and other facilities, and to promote the improvement of the social benefit of construction land use, so as to improve the overall CCD of such areas. For higher coordination regions, they should strengthen the construction of resource-based industries to continue industrial transformation, and promote the orderly development of the cities, and the development of green, efficient, clean production, and a circular economy. The optimal coordination areas such as Lanzhou and Xining should take the opportunity of the new round of western revitalization to strengthen the intensive use of construction land and improve the infrastructure construction level to accelerate the transformation and upgrading of traditional advantageous industries, to increase the proportion of the tertiary industry, and to consolidate the industrial foundation and employment support for urban development, so that the urban population capacity can be improved, and realize the transformation of the multi-functional benefits of construction land utilization a at a more highly coordinated level.

#### **6. Conclusions**

This paper uses panel data of 39 county-level cities of the LXUA in China from 2005 to 2018 to construct a three-dimensional system of CLUB, covering society, economy, and eco- environment. We analyzed the temporal and spatial evolution characteristics of CLDI and CLUB and explored the coupling coordination relationship among the benefits. The main conclusions are as follows: (1) The CLDI generally presents a dual-core spatial distribution with Lanzhou and Xining urban areas as the core. In the center of the city, the CLDI is greater. In contrast, in the edge cities, as the distance increases, the CLDI gradually decreases. (2) The social benefit of construction land use showed a "convex" shaped spatial distribution pattern of "high in the middle and low in the east and west wings". The economic benefit was basically the same as the social benefit in terms of spatial distribution, with great regional differences. The eco-environmental benefit was relatively high andregional differences were small. (3) The comprehensive benefit construction land use had significant spatial differences, and generally presented a gradient spatial structure of "core cities > node cities > general county towns", and the comprehensive benefit was mainly lower level and low level. (4) From 2005 to 2018, the interaction between the social, economic, and eco-environmental benefits of construction land use developed in a coordinated direction in terms of time evolution and spatial correlation, but generally the CCD was still at a medium-low level, and the spatial clustering features were significant.

This paper conducts a more comprehensive analysis of the spatiotemporal evolution of CLDI, CLUB, and the coupling coordination relationship of the multi-functional use benefits of construction land in LXUA, which is of great significance to the sustainable development of future land use of different types. However, there are still some deficiencies that need to be pointed out: (1) Based on the availability of county data, the selection of indicator systems still can be improved. (2) With the update of technical means such as big data, the "flow" element is becoming more and more important. However, this paper does not consider the impact of "data flow" and "feature flow" on regions. (3) Compared with other studies [57–59], this study lacks consideration of more refined data and techniques, and impact mechanisms. In the future, on the basis of comparing and learning from the more mature urban agglomeration paths, a more in-depth analysis and research should be carried out in combination with the characteristics of the study area and the shortcomings of this paper.

**Author Contributions:** Conceptualization, W.Z. and P.S.; methodology, W.Z.; software, W.Z.; validation, W.Z., P.S. and H.T.; formal analysis, W.Z.; investigation, W.Z.; resources, P.S.; data curation, W.Z.; writing—original draft preparation, W.Z.; writing—review and editing, W.Z. and H.T.; visualization, W.Z.; supervision, P.S.; project administration, P.S. and H.T.; funding acquisition, P.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 41771130, 42161043).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the first author.

**Acknowledgments:** We sincerely thank the reviewers for their helpful comments and suggestions about our manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


## *Article* **Downscaling Switzerland Land Use/Land Cover Data Using Nearest Neighbors and an Expert System**

**Gregory Giuliani 1,2,\*, Denisa Rodila 1,2, Nathan Külling 1, Ramona Maggini <sup>3</sup> and Anthony Lehmann <sup>1</sup>**


**Abstract:** High spatial and thematic resolution of Land Use/Cover (LU/LC) maps are central for accurate watershed analyses, improved species, and habitat distribution modeling as well as ecosystem services assessment, robust assessments of LU/LC changes, and calculation of indices. Downscaled LU/LC maps for Switzerland were obtained for three time periods by blending two inputs: the Swiss topographic base map at a 1:25,000 scale and the national LU/LC statistics obtained from aerial photointerpretation on a 100 m regular lattice of points. The spatial resolution of the resulting LU/LC map was improved by a factor of 16 to reach a resolution of 25 m, while the thematic resolution was increased from 29 (in the base map) to 62 land use categories. The method combines a simple inverse distance spatial weighting of 36 nearest neighbors' information and an expert system of correspondence between input base map categories and possible output LU/LC types. The developed algorithm, written in Python, reads and writes gridded layers of more than 64 million pixels. Given the size of the analyzed area, a High-Performance Computing (HPC) cluster was used to parallelize the data and the analysis and to obtain results more efficiently. The method presented in this study is a generalizable approach that can be used to downscale different types of geographic information.

**Keywords:** land cover; land use change; downscaling approach; Switzerland; geographic information system; aerial photo interpretation; topographic map; inverse distance weighting; expert system

#### **1. Introduction**

*1.1. Pressures on Land Resources in Switzerland*

The Swiss Federal Department of Environment, Transport, Energy and Communications (DETEC) in its 2016 Strategy stated that by 2030, Switzerland is aiming at becoming a sustainable country while remaining an attractive and competitive business location with a high quality of life [1]. This ambitious objective is presently challenged by several trends such as population growth, increased mobility, energy demand, high consumption of resources, urbanization, loss of biodiversity and associated ecosystem services, and the digitalization of society along with related big data [2]. These trends have an important impact on the environment. Therefore, protecting the environment is a central mission for the Swiss Government, who wants to promote and adopt more sustainable approaches for the exploitation of natural resources [3]. To this end, actions such as protecting natural resources, improving urban planning, reducing emissions of greenhouse gases, preserving water quality, retaining biodiversity and ecosystem services, protecting soils, and preserving countryside are essential [4–6].

All these trends are also placing unprecedented demands on land. Between 1985 and 2009, 15% of the country's surface area changed [7]. Settlement and urban areas have

**Citation:** Giuliani, G.; Rodila, D.; Külling, N.; Maggini, R.; Lehmann, A. Downscaling Switzerland Land Use/Land Cover Data Using Nearest Neighbors and an Expert System. *Land* **2022**, *11*, 615. https://doi.org/ 10.3390/land11050615

Academic Editors: Víctor Hugo González-Jaramillo and Antonio Novelli

Received: 22 February 2022 Accepted: 20 April 2022 Published: 21 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

expanded, agricultural areas have been lost, forested areas have increased, and glaciers have receded [8–10]. During the last 50 years, it is estimated that human activities have affected globally about 83% of the terrestrial land surface and have degraded about 60% of services provided by ecosystems. Land degradation is now at a critical point and will undermine the well-being of 3.2 billion people by 2050 [11]. Consequently, Land Cover (LC) and Land Use (LU) changes are considered as a major tangible indicator of the human footprint [12].

To preserve its potential to deliver goods and services, land should be efficiently and sustainably managed. National policies such as the Green Economy, the Spatial Planning Act, the Spatial Strategy for Switzerland, the Sustainable Development Strategy 2016–2019, or the Strategy on Biodiversity are essential components to support this vision [13]. They generally acknowledge that a given area of land can offer many environmental, social, cultural, and economic benefits at once. However, most ecosystems are being degraded by unsustainable exploitation, fragmentation, urban growth and development of transport, and energy networks. This reduces the spatial and functional coherence of the landscape and consequently, degraded ecosystems are unable to provide the same services as healthy ecosystems [14].

Detailed and accurate knowledge on Land Use and Land Cover Change (LU/LCC) is crucial for many scientific and operational applications, such as watershed analyses [15,16], land use impact on stream ecology [17–19], species and habitat distribution modeling [20], dynamic modeling of species migration [21,22], reserve site selection [23], impact assessment on biodiversity [24], land use planning [25], or monitoring of land use changes [26]. LU/LC affects many aspects of policy and decision-making processes related to climate, water, biodiversity, ecosystems, agriculture, or disasters. Additionally, LU/LCC assessment also contributes to many Multilateral Environmental Agreements (MEAs) and Global Environmental Goals (GEGs) to guide and assess progress toward policy outcomes [27,28]. The importance of sustainable management of land resources is recognized in regional and global policies such as the 2030 Agenda for Sustainable Development, which contains land-related targets and indicators under 14 out of the 17 Sustainable Development Goals (SDGs) [29–31]. Many land organizations and stakeholders are committed to fully implement the SDGs and to monitor the land-related indicators to promote responsible land governance. Land is a significant resource for many sectors; timely and high-resolution LU/LC data therefore constitute critical information for the achievement of the SDGs [32]. Accurate and up-to-date LU/LCC information and related changes are also the base of a sustainable development assessment [33] based on structural (both temporal and spatial) and functional (social, ecological, economical) attributes of the landscape [34]. A supplementary challenge is represented by the spatial scale at which the assessment is performed. For each problem under study, an appropriate scale must be identified, especially when relating ecological processes to landscape patterns [35].

#### *1.2. Land Use/Land Cover Data in Switzerland*

LU/LCC is increasingly acknowledged as both a driver and a consequence of climate and biodiversity changes [36–38]. This important role has been featured by the fact that land cover is considered as an Essential Climate Variable (ECV), a supplementary Essential Water Variable (EWV), and a candidate Essential Biodiversity Variable (EBV) [39–42]. LU/LCC affects the biophysics, biogeochemistry, and biogeography of both the atmosphere and biosphere, with important consequences for human well-being. Consequently, accurate and timely information is necessary for understanding the impact of LU/LCC variations on the structure and functioning of ecosystems, as well as provision, support and regulation of goods and services [29,43,44]. However, it is recognized that inadequate information on LU/LC and its change over time is a recurrent and common problem that prevents policymakers from making sound, informed decisions [27,45–47]. Currently, the official LU/LC information in Switzerland (Arealstatistik) is updated approximately every 6 to 8 years and derived by visual interpretation of aerial photographs where an LC and an LU category are assigned to each intersection point of a regular 100 m grid [8]. Although this data set is very useful thanks to its thematic richness, neither its low spatial resolution nor the update frequency allow for providing accurate and timely information to depict and understand the dynamic of LU/LCC and the related impact across the country [48,49]. Accurate LU/LC change assessment and effective LU/LCC projections require higher spatial (e.g., 30 m) and temporal (e.g., yearly) data products to build consistent time series [47,50].

#### *1.3. Downscaling as a Possible Approach for High-Resolution LU/LC Data*

Besides traditional remote sensing approaches, such as unsupervised, supervised, or object-based classifications [51–54], more advanced techniques include machine or deep learning [55–58], new sensors with higher spatial and spectral resolution such as Sentinel-2 [59,60], and automated procedures to reduce the time-consuming process of manual verification of data [61,62]; a possible alternative to generate high-resolution LU/LCC data [63] is represented by downscaling techniques [64]. In many disciplines, downscaling is used to derive local scale maps from information available at coarser resolution. Climatologists refer to statistical downscaling [65] to describe this general approach that has been widely used not only for temperature and precipitation information [66,67] but also for wind speed [68] and air humidity [69]. In turn, downscaled climatic information is used in many different applications, such as hydrological modeling [70–72], species distribution modeling [73], and geological risk assessments [74]. However, downscaling of LU/LCC data is not very common and has not been widely applied [64,75,76].

Several statistical approaches have been used for downscaling. For instance, Barodssy et al. [77] used fuzzy rule-based models to predict frequency distributions of daily precipitation; Bürger and Chen [78] compared regression methods to derive river runoff from large-scale climatic scenarios; Biau et al. [79] used geostatistical methods (kriging) to estimate rainfall; and Coulibaly et al. (2005) [67] investigated the use of temporal neural networks to downscale temperature and precipitation.

Downscaling is not restricted to climatic data and has been used with remote sensing data to derive, for instance, soil moisture maps [80]. Species distribution modeling can also be defined as a general downscaling approach that predicts species distributions from point observations combined with spatially explicit environmental predictors [81]. The term "downscaling" can also be used when creating a land use map from combined input layers at various scales. For example, Remm [82] used case-based predictions to map the distribution of habitat classes from Landsat 7 ETM imagery, grayscale and color orthophotos, an elevation model, a digital base map, and a soil map.

Case-based algorithms are problem-solving methods that learn from experiences at a low level of generalization [83]. They can be considered as an Artificial Intelligence (AI) method that derives results from the data as directly as possible, without the formulation of an intermediate model. Machine learning (ML) specialists distinguish between lazy learning, which typically combines information during the problem-solving phase, and eager learning, which tends to derive a generalization and forget about raw observations after the learning phase [56,83,84]. Remm [82] argued that case-based methods represent a promising alternative for a large range of downscaling problems such as habitat mapping and the prediction of species' potential distributions, especially with large and complex datasets where generalization is difficult.

Based on these considerations, the aim of this paper is to present a lazy learning method, which could be assimilated to a case-based approach, for downscaling LU/LC information for Switzerland from a 100 m lattice of points to a 25 m resolution grid, taking advantage of an existing 1:25,000 digital base map and building an expert system defining possible correspondences between the base map and the land use categories. This method is then applied for three different periods of time to assess land use and land cover change.

#### **2. Materials and Methods**

While land use and land cover maps are commonly derived from remote sensing or photointerpretation [85], traditional base maps of Switzerland have been available for more than a century and have been provided in digital format since the year 2000. LU/LC maps derived from the classification of remotely sensed images often have a "salt and pepper" appearance that does not meet end-user demand [26]. Land use derived from aerial photo interpretation can define many different classes of land use/cover categories, but the production is rather time-consuming. National base maps usually lack the thematic details that can be obtained from aerial photointerpretation but generally have an excellent geographic precision to define landscape patches and linear features. Hereafter will be presented the data inputs (Table 1), the downscaling methodology, and the expert system, together with its implementation and validation strategy.

**Table 1.** Data input sources.


#### *2.1. Data Inputs*

#### 2.1.1. Land Use Statistics

LU/LC data are generated by visual interpretation of aerial photographs taken from a Federal Office of Topography (swisstopo)'s aircraft flying at an altitude of 5000 m and taking photos to regularly cover, over a 6-year period, the entire surface of Switzerland [8]. LU/LC maps are obtained by visually interpreting and assigning a LU/LC category to each point of a regular 100 m lattice laid over the Swiss territory, for a total of more than 4 million points over the country.There are three official classifications: (1) the standard nomenclature NOAS04 (72 basic categories which are a combination of LC and LU, 17 and 27 aggregation classes and 4 main domains); (2) the Land Cover nomenclature NOLC04 (27 basic categories and 6 main domains); and (3) the Land Use nomenclature NOLU04 (46 basic categories, 10 aggregation classes and 4 main domains). Three time periods are currently available (1979/85, 1992/97, 2004/09), and the latest version has just been finalized (2013/18) [7].

Strictly speaking, this dataset is not a LU/LC map, because its categories are assigned to points at the intersection of a 100 m grid rather than indicating the predominant LU/LC within each hectare square (Figure 1). It was developed as land use statistics over relatively large zones rather than as an LU/LC map per se. It tends, however, to be often used as an LU/LCC map in many applications [86] and remains the most exhaustive source of LU/LCC information for Switzerland. Even if this dataset is thematically more precise than the classification commonly used in Europe—the Coordination of Information on the Environment Land Cover (CORINE Land Cover, CLC), which has 44 classes [87]—it suffers from a low spatial and temporal resolution. Indeed, LU/LC units are coarse with a spatial resolution of 1 hectare, and a lot of information is therefore aggregated with a large degree of generalization. Consequently, various landscape features, qualities, particularities, and configurations cannot be correctly represented.

**Figure 1.** Original LU/LC statistics with 72 classes (Arealstatistik) at a 100 m resolution from the 2013–2018 period.

#### 2.1.2. National Base Map (Land Cover)

The Swiss Federal Office of Topography (swisstopo) provides digital versions of national topographic base maps at a 1:25,000 scale as a landscape model in a vector format. The data model used until 2011 was called Vector 25 and was later replaced by the Swiss Topographic Landscape Model (TLM3D). Both include millions of natural and artificial landscape features, together with their position, shape, type and many other attributes [88]. This land cover information is defined in 29 categories (Figure 2). Data on linear features such as rivers, roads and rails can also be obtained separately. The national base maps (TLM3D) are the geographically most precise source of land cover information for the entire country with a geometric accuracy for different landscape features between 0.2 m and 3 m that is partially updated every year.

#### 2.1.3. Resolution

For the analysis, all available datasets were rasterized either at a 25 m and/or at a 100 m resolution to allow raster overlays. With a surface of 42,000 km2, Switzerland can be described with approximately 4 million pixels at a 100 m resolution and with 64 million pixels at a 25 m resolution.

#### 2.1.4. Data Quality of Inputs

The two main data inputs of this study are of the highest possible quality from the two main national producers of geospatial data. The first one, land use statistics, has been developed by the Swiss Federal Office of Statistics with state-of-the-art photointerpretation methods, resulting in a very high quality of thematic resolution in the selection of land use classes on each hectare point of the country. The second input, the national base maps, are produced by the Swiss Federal Office of Topography and represent the official geographic representation of the country at a 1:25,000 scale with very high spatial accuracy but a less developed thematic resolution than the first input.

**Figure 2.** LC information with 29 classes from the National Base Map swissTLM3D for 2021.

*2.2. Downscaling Algorithm and Expert System*

#### 2.2.1. Downscaling Algorithm

The approach used for downscaling the existing land use information from 100 to 25 m relies on both the geographic precision of the 1:25,000 national topographic base map rasterized at 25 m and the detailed LU/LC categories obtained from the land use statistics available at 100 m. The developed algorithm uses inverse distance weighting combined with an expert system to assign reasonable land use categories at a finer scale (Figure 3).

**Figure 3.** Downscaling land use statistics (2013–2018) in the area of Zermatt from 100 m to 25 m resolution using inverse distance weighting, expert knowledge, national base map, transport and river information at 25 m. (**a**) Expert system (details in Figure 4), (**b**) 1:25,000 base map, (**c**,**c***-* ) hectare land use information, (**d**,**d***-* ) downscaled land use, (**e**) 1:25,000 linear features (rivers, roads and rails), (**f**,**f***-* ) overlay of linear features to downscaled land use at 25 m resolution. Downscaling process (p1) and linear features addition (p2).


decision. Ten land use categories remain unmatched (shaded). Three categories correspond to punctual or linear features (dark grey), and seven categories (light grey) correspond to different types of buildings that were merged with their surroundings (green). A full version of the expert table is available in the Supplementary Material Table S1.

The main steps of the downscaling methodology are the following (data preparation: 1–2–3; process [1] of downscaling: 4–10; process [2] for linear features: 11):


The Inverse Distance Weighting (IDW) calculates a scaled distance to each of the 36 Landuse100 neighbors around the pixel under investigation, and does so in two spatial dimensions (Equation (1)):

$$dj(i) = \text{sqrt((x25i - x100j)^2 + (y25i - y100j)^2)/maxrange} \tag{1}$$

where *j* spans from 1 to 36 nearest neighbors and *i* represents each visited pixel at a 25 m resolution and *maxrange* the maximum distance between these pixels.

Then the inverse distances are summed up by land use category (Equation (2)):

$$D\_k(i) = \sum\_{j=1}^{\infty} \left( \frac{\lambda\_j(k)}{d\_j(i) + s} \right) \tag{2}$$

where *λj*(*k*) = 1 *if landusetypej* = *k* <sup>0</sup> *otherwise* , and *<sup>s</sup>* is smoothing factor.

Finally, the category scoring the highest sum of inverse distances is assigned to the pixel under investigation (Equation (3)):

$$K(i) = \text{greatest}\_k(D\_k(i))\tag{3}$$

IDW is generally used for interpolating cardinal discrete values (e.g., temperatures, precipitations) but can be also used on nominal values to spatially interpolate missing LC data [89], create super-resolution LC maps [90], or rescale LC data [91]. IDW assumes that values that are close to one another are more similar than those that are at a greater distance. In the proposed method, measured values surrounding the prediction location are considered but then they are spatially weighted (i.e., taking into account the distance of each pixel and the frequency of classes) and constrained (i.e., by the expert system). Consequently, the proposed downscaling method is filling a geographically very precise map (the base map) with information on the most plausible LU/LC category found in the neighborhood by reducing the possible choices with an expert table. This approach has the advantage of fully respecting the geographical quality of the base map while enriching it with a better thematic resolution.

#### 2.2.2. Expert System

An expert system was used to constrain the possible choices when using information from BaseMap25 to select the most appropriate Landuse100 category (Figure 4). For instance, if a forest patch is defined for a particular site in BaseMap25, the expert system constrains the choice of Landuse100 categories to include only those related to forest. Furthermore, the expert system can also select a default Landuse100 category when the distance-based algorithm fails to make a clear selection because of the lack of eligible land use categories within the searching neighborhood.

#### *2.3. Implementation*

Although the algorithm used is relatively simple, it requires extensive calculations on several large grids (25 m resolution grid for Switzerland contains approximately 64 million pixels), which is time-consuming. To ensure the portability of the code on different platforms, the algorithm has been implemented using a set of open-source software and libraries. The programming language is Python 3 [92] using the PyCharm Community Edition (CE) Integrated Development Environment (IDE) [93]. The implemented algorithm relies on the following libraries: (1) the Geospatial Data Abstraction Library (GDAL) [94] to handle raster data; (2) Numpy [95] and math for mathematical operations on large multi-dimensional arrays; and (3) xlrd [96] and pandas [97] for reading and formatting information from Excel files.

To ensure a fast and efficient processing, the algorithm has been parallelized and executed on the High-Performance Computing cluster at University of Geneva [98], allowing to process the entire 64 million pixels grid significantly faster. The analyzed area is subdivided in tiles and the processing of each tile is performed on a different node on the cluster. The surface of Switzerland was divided in 900 tiles (30 × 30), which were associated to 900 jobs in the cluster, using SLURM job arrays. The processing time is on average less than 2 h per tile, but the actual processing time of a single tile is strongly dependent on the user priority and on the cluster load at the execution time. When the results of all the jobs are available, a final script is launched to assemble the 900 processed tiles together in a single image.

The code is freely available on GitHub [99] under an Apache 2.0 license, and a static version can be downloaded from the University of Geneva digital repository [100] under a CC-BY 4.0 license.

#### *2.4. Validation and Accuracy Assessment*

The categories of the downscaled LU/LC maps at 25 m (with or without linear features) for the three periods were compared with the categories of the original LU/LC point statistics. A random sample of 500,000 out of 4,129,070 points was used for this assessment. Multi-class classification metrics were used to assess the efficiency of the algorithm to classify each LU/LC category [101–103]. Results are shown as F1-score, which is obtained by harmonic mean from the recall and precision, and Cohen's kappa coefficient [104]. Results were expressed as class F1-scores and weighted means for each dataset, and kappa coefficient for each dataset. Barplots and boxplots are used to display these metrics and the main misclassifications per category. The relative surface area per category was also compared to that of the original Landuse100 dataset [64,105].

#### **3. Results**

Results from the land use downscaling are best appreciated at a regional scale (Figures 5 and S1), where the map produced by the algorithm not only preserves the geographic precision of the original BaseMap25 map but also inherits the higher definition of categories (66 categories in the downscaled map compared to 29 categories in BaseMap25). Seventy-two categories were originally found in Landuse100, but a few groupings had to be performed on some less important categories (e.g., buildings and their surroundings). The new map has a much finer grain and can be used for Geographical Information System (GIS) overlays at much finer scales, where the 16-fold increase in the density of pixels gives substantially improved results. Fine-scale details on roads and river networks were also maintained, whereas they are often difficult to represent adequately in coarse resolution raster formats (Figure 6).

**Figure 5.** Land use downscaling at 25 m of Switzerland for the 2013–2018 period. Swiss regional parks boundaries in red. (See Figure 6 for legend interpretation.)

The average percentage of similarity for the three periods evaluated between the original points found in the Landuse100 statistics and the resulting Landuse25 classifications is 71% when assessed without roads, rivers, and rails, and 68% when assessed with these linear features. Several categories have a correspondence of 70% or greater. Overall F1-scores and kappa coefficient are high (0.69 and 0.67, respectively). Individual classes F1-scores range from 0.007 to 0.98 depending on the category (Figure 7). Weighted mean F1-scores per dataset range from 0.67 to 0.70. Kappa scores range from 0.65 to 0.69. Categories from the Landuse100 classification that show zero correspondence are those that were not predicted, such as flood protection structures (63), field fruit trees (38), or grooves and hedges (58), and were thus discarded from the validation analysis. The categories that are best respected by the downscaling are those that have a large spatial representation and those that were imposed from the topographic map.

**Figure 6.** Land use downscaling from 100 m to 25 m and overlay of linear features for state of Geneva for the three periods.

**Figure 7.** Weighted boxplots showing the F1-scores from the 72 downscaled classes (Landuse25) compared to the original classes (Landuse100), for each period, with ("all") and without linear

features. Numbers along the boxplots represent the land use categories and their size the proportion of each category in the dataset. Blue represents categories attributed based on "swisstopo" original classes; black represents categories attributed based on IDW of the 36 nearest neighbors. Values under each boxplot are overall weighted mean F1-score (F1) and kappa coefficient (k) for the considered year and dataset.

The percentages of misclassification for the downscaling of 2018 without linear features are represented in Figure 8. The major misclassifications concern unproductive grasslands (15%), farm pastures (12%), and meadows (10%). The classes in which the pixels were misattributed are also presented (color in the box). This figure will help in understanding the behavior of the downscaling algorithm to further improve it.

**Figure 8.** Percentage of all misclassifications for the 17 main classes of the 2018 downscaled LU without linear features. Legend lists categories that were misattributed.

The absolute surface area difference from the downscaled map (Landuse25) categories compared to the original Landuse100 map ranges from −87,585 ha for the category "unproductive grass and shrubs" to +118,174 ha for the category "normal dense forest" (Figure 9). On average, a difference of 532 ha per category is observed. Underlaying reasons for large discrepancies in surface areas (misattributed categories) are further displayed in Figure 8. The relative surface area difference (proportion from downscaled map compared to original) ranges from −100% for small categories that were not included in the algorithm to +355% for cat. 14, "surroundings of unspecified buildings".

**Figure 9.** Absolute surface area difference from the downscaled map (Landuse25) categories compared to the original Landuse100 map.

#### **4. Discussion**

The method presented in this study is a generalizable approach that can be used to downscale different types of geographic information such as results from photo interpretations, classifications of remotely sensed images, or vegetation and soil maps. The general idea is to use the nearest precise point observations to define the attribute of a target pixel at a finer resolution within an area defined by a high-resolution land cover map. Other

co-variables such as remote sensing indices (e.g., Normalized Difference Vegetation Index (NDVI)), topographic position, slope, or orientation could also be used to help model the most probable LU/LC category and could certainly further improve the present method.

In the case study presented here, the use of a national topographic base map at a 1:25,000 scale to define main land cover patches guarantees a perfect overlay with official maps. Indeed, topographic base maps are now available in digital format and are widely used as reference maps in most studies and field work. The possibility of matching the geometry defined by base maps with detailed land use categories reinforces the chance of uptake from end-users. However, some recent changes in LU/LC that might be visible from remote sensing images could be lost if involving changes between incompatible land use categories as defined by the expert system.

One of the main strengths of the proposed approach is that it avoids the "salt and pepper" effect generally resulting from the classification of remotely sensed images [26], except in the case of object-oriented classifications [106]. We believe that our approach could be particularly useful in this respect and could be used to smooth out land use maps obtained from remote sensing classifications. In such a case, land use statistics would be replaced by the land use classification obtained from supervised or unsupervised classification of the remotely sensed image, whereas the land cover information would still be retrieved from the national base map. Moreover, with the increased spatial resolution and the removal of the "salt and pepper" effect, LU/LC changes are more evident and can be assessed more easily (Figure 10).

**Figure 10.** Comparison of the 1997 and 2018 LU/LC data with original 100 m resolution (**above**) and the downscaled outputs at 25 m (**below**) for the Bulle and Freiburg cities. Urban expansion (as well as other land cover changes and the smoothing effect) of this economically active region can be more easily visualized with the improved spatial resolution.

The proportion (71%) of exact matches between the land use categories recorded at the 100 m resolution lattice points of the Landuse100 statistics and the corresponding points in the downscaled map (25 m resolution) is particularly encouraging. Discrepancies appear to result mainly from the change in scale and geometric precision brought in by the 1:25,000 base maps. The fundamental difference in approach between the land use statistics (punctual) and the topographic map (surface) is also important. A visual check confirms that most divergences mainly occur along boundaries of patches of different land use and along linear and small features. By using an inverse distance calculation—which gives a higher weight to close-by information—to choose the best land use category for each pixel, we greatly favored the retaining of input land use categories at an observed location in the result. This effect can be modulated by the smoothing factor in Equation (2).

The approach could be used at even finer scales by rasterizing for instance the topographic maps at 10 m instead of 25 m or using a regional map at a finer scale (1:10,000). Getting accurate land use information is crucial for calculating landscape indices such as those obtained from FRAGSTAT [107]. As demonstrated by Uuemaa et al. [108], changing grain size can have a significant effect on many landscape metrics. Therefore, one should pay attention to the scale at which a given landscape metric is calculated, and what grain size of land use maps is best adapted for the purpose, as there is no single land use scale which is "best" for observing and managing changes in land use patterns.

The implementation of a case-based reasoning algorithm in Python code proved to be a very efficient approach given the very large number of pixels (64 million) to handle but required the development of a purpose-written software. The availability of tools for implementing case-based approaches within existing GIS packages would certainly be very useful for many applications. One such tool has been developed by Remm [82] using Microsoft Visual Studio.NET, and this can predict several response types (binomial, multinomial, continuous, and complex) based on continuous and categorical predictors.

Riitters [109] mentions the practical trade-off that exists between the generality, precision, and realism of a method, as previously proposed by Levins [110], and emphasizes that generality and realism should be maximized at national scales, while precision should be optimized at local scales. Indeed, the use of the same raw data to develop models at different scales will enhance thematic and geographic comparisons between disciplines and across scales. Methods to upscale and downscale data are therefore central to local, national, and global environmental assessments.

The presented approach contributes to tackle the need of increased spatial resolution of LU/LCC information. However, the temporal issue remains a major problem. Indeed, land cover changes are caused by either natural or anthropogenic sources such as climate change, demographic growth, and economic growth [111,112]. Therefore, the state of land cover is highly dynamic and involves an inherent challenge for its mapping and monitoring that remains not adequately addressed [85]. Timely and reliable information on land cover change is crucial to efficiently mitigate the negative impact of environmental changes [113,114]. Traditional environmental data collection (e.g., field data collection) suffers from many shortcomings, among them the data inconsistency caused by changes in reporting methodologies through time, the gaps (missing measurements) in data series, and the fact that it is notoriously time-consuming and labor-intensive [115]. One of the main advantages of remotely sensed data collection is that it can provide a synoptic and repetitive view of a given area/region. With the opening of different Earth Observations (EO) data archives such as Landsat [112,116,117], it becomes possible to build consistent time series (i.e., image of the same location at regular intervals) to compare different periods of time and derive trends [118–120].

At the national scale, various attempts have been made to use satellite imagery, but unfortunately, no processes are currently in place to routinely generate accurate, consistent, and regular LU/LCC data. These large volumes of freely and openly available EO data are still underutilized and are not effectively used for national environmental monitoring. Therefore, mapping and monitoring LU/LC changes remain as challenges that are not

adequately addressed at the national scale, and new methodologies are required to produce consistent and reliable yearly, medium-to-high-resolution (spatial, temporal, thematic) time series of LU/LCC data and projections of their (future) change across Switzerland to inform national and regional environmental policies and planning.

With the opening of different medium-to-high-resolution satellite EO data archives such as Landsat or Copernicus and the development of advanced data science techniques (e.g., Big Data, Artificial Intelligence, High-Performance Computing), it now becomes possible to build consistent time series of LC, to investigate the spatio-temporal dynamics of LC, and to perform quantitative assessments of LC dynamics, by comparing different periods of time, deriving trends, and determining environmental trajectories [121]. This can generate a consistent and reliable LU/LC time series that will help to understand the evolution of LU/LCC in Switzerland over the last 30 years and to model future changes according to plausible scenarios by 2050 [122,123].

To reach this objective, an essential pre-condition to support user applications and generate usable information products is to facilitate data access, preparation, and analysis. The systematic and regular provision of Analysis-Ready Data (ARD) can significantly reduce the burden of EO data usage. To be considered as ARD, data should be processed according to a minimum set of requirements (e.g., radiometric and geometric calibration; atmospheric correction; metadata description) and organized in a way that allows immediate analysis without additional effort [124,125]. In Switzerland, more than 37 years of satellite EO Analysis-Ready Data over Switzerland are made available by the Swiss Data Cube [126,127]. The increasing availability of EO data together with the improved computing and storage capacities allow monitoring, mapping, and assessing LU/LC and its change over time on large areas in a consistent and reliable manner. This favors the development of annual, high-quality LU/LC products based on time series data and can inform on class stability and transitions [128].

#### **5. Conclusions**

The proposed downscaling approach allowed for combining the geographic precision of existing topographic base maps with the thematic details of photo-interpreted land use statistics. Improved land use maps open the door to more accurate watershed analyses, species and habitat distribution modeling, and species dynamic models of migration. Accurate land use information is also the base for developing sustainable development indicators defining the structural and functional attributes of landscapes.

The proposed approach could be implemented for three time periods. It allowed for an efficiently downscaled LU/LC at a spatial resolution that is more suitable for environmental change monitoring. The increased spatial resolution removed the "salt and pepper" effect, and consequently, LU/LC changes are more evident. However, the changing definition of the base map input across the years resulted in discrepancies in the resulting downscaled maps that are refraining their use for land use change analyses. We therefore recommend using the original land use statistics data for trend analyses and our downscaled data for GIS analyses at each time period. The presented approach contributes to tackling the need for increased spatial resolution of LU/LC information, but the temporal issue is still a major problem. The use of dense time series of satellite data such as those provided in the Swiss Data Cube can be a promising solution to investigate to obtain high spatial and temporal LU/LC data over the country.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/land11050615/s1. Figure S1. Final result of the downscaled Land Use/Land Cover at 25 m for the period 2013/18. Table S1. Full version of the expert table.

**Author Contributions:** Conceptualization, A.L. and R.M.; methodology, A.L.; software, A.L., G.G. and D.R.; validation, A.L. and N.K.; formal analysis, A.L., G.G. and D.R.; data curation, A.L. and N.K.; writing—original draft preparation, G.G., A.L. and R.M.; writing—review and editing, all; visualization, A.L. and N.K.; supervision, A.L.; project administration, A.L.; funding acquisition, A.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** Support by the Swiss Federal Office for the Environment with the ValPar.ch project is gratefully acknowledged.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The results presented in this study are openly available in Yareta, the University of Geneva Digital Repository, at: https://doi.org/10.26037/yareta:dlx3hu54jfa3ne3 c2xjfcnqpxm (accessed on 21 February 2022) and in the ValPar project Spatial Data Infrastructure at: http://valpar.unige.ch:8080/geonetwork/srv/eng/catalog.search#/metadata/da78fa9f-a756-4 538-a28a-8be7d53d4676 (accessed on 21 February 2022). The source code of the algorithm is available at: https://github.com/ggiuliani/LULCdown (accessed on 21 February 2022). Publicly available datasets were analyzed in this study. These data can be found here: Land Use Statistics (https: //www.bfs.admin.ch/bfs/de/home/statistiken/raum-umwelt/erhebungen/area.html, accessed on 21 February 2022); Base Maps (SwissTLM3D—https://www.swisstopo.admin.ch/en/geodata/ landscape/tlm3d.html, accessed on 21 February 2022).

**Acknowledgments:** Swisstopo and the Swiss Federal Office for Statistics are gratefully acknowledged for making the input data freely available.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

