Exploring the Applicability of Self-Organizing Maps for Ecosystem Service Zoning of the Guangdong-Hong Kong-Macao Greater Bay Area

: Sustainability is one of the major challenges in the 21st century for humanity. Spatial zoning of ecosystem services is proposed in this study as a solution to meet the demands for the sustainable use of ecosystem services. This study presented a workﬂow and performed an exploratory analysis using self-organizing maps (SOM) for visualizing the spatial patterns of the ecosystem service value (ESV) of the Guangdong-Hong Kong-Macao Greater Bay Area (GBA). The zoning was performed based on 11 types of ecosystem services, resulting in 11 ecosystem service zones. Each of the zones derived has its unique characteristics in terms of the dominating ecosystem service types, ESV, land use/land cover patterns, and associated human activity levels. It is recommended that reasonable and effective utilization of the ecosystem services in the GBA should be based on its zonal characteristics rather than haphazard exploitations, which can contribute to the sustainable economy and environment of the region. The applicability of SOM for the GBA ecosystem service zoning has been demonstrated in this study. However, it should be stressed that the method and workﬂow presented in this study should mainly be used for supporting decision-making rather than used for deriving gold-standard zoning maps.


Introduction
A majority of the mega-city regions in the bay areas of the world (e.g., the San Francisco Bay Area and the Tokyo Bay area) are highly developed in terms of economy and technology [1,2]. Mega-city regions account for more than 60% of the global economy and there has been an increasing trend of people migrating from inland regions to these highly developed coastal areas; therefore, the sustainable development of such regions is highly important [2]. China State Council proposed the development of its Guangdong-Hong Kong-Macao Greater Bay Area (GBA) in 2015. Developing this emerging mega-city region is a national strategy of China to create a model for its full-scale opening-up.
The GBA covers an entire area of 65,000 km 2 and the population is around 70 million as of the end of 2018. The GBA was once not as integrated and economically important as it is nowadays. The GBA has evolved to an economy led by manufacturing in the last decade of the 20th century from an economy dominated by agriculture, and now it is moving toward an economy based on strong innovation, advanced technologies, and an internationally oriented business environment [3]. The GBA has become the biggest metropolitan region in the globe measure specific attributes of a system, reducing the complexity of the system and describing it in a simplified manner [21]. Indicators enable the inclusion of expert knowledge in decision-making processes and establish a linkage between actions and targets, the fluctuations of which demonstrate the variation in the phenomena that they represent [21]. Indicators can be adopted individually or as a combined indicator set [22].
With the 11 ESV indicators, the objective of this study is to perform an exploratory analysis using self-organizing maps (SOM, a machine learning method based on artificial neural networks [23,24]) for spatial zoning and visualization of the GBA ecosystem services. Zoning-dividing a study area into multiple homogeneous regions based on regional characteristics-is an important direction of geographic studies and can be useful in regional planning and development [25]. Zoning of the ecosystem services tailors targeted monitoring, policy-making, decision-making, and management that enable effective regional planning for the GBA's sustainable development. As a result of the increasing complexity of geographic data and geospatial analysis, SOM has been proven to be an effective approach to clustering analysis and pattern detection [26,27]. Various studies have proven the ability of SOM to perform zonings of geospatial phenomena, such as water quality zoning [28], land use zoning [29], streetscape zoning [30], tectonic zoning [31], and landslide susceptibility zoning [32]. SOM is also an approach that can be readily used by researchers and practitioners. Note that as the ecosystem services in the GBA are spatially and temporally dynamic, we do not intend to provide fixed and long-lasting zoning results, but only to explore the applicability of SOM for such a zoning purpose. The zoning results and maps must be updated and fine-tuned depending on changes in the input ESV indicators and zoning requirements.

Spatial Zoning for Facilitating Sustainable Ecosystem Services
Sustainability is one of the major challenges in the 21st century for humanity. Many researchers have criticized that the existing studies on ecosystem system services have not led to more significant contributions to decisions about sustainable development [33]. Bennett and Chaplin-Kramer [33] explored the "Science for the sustainable use of ecosystem services" and suggested three directions for improving the use of ecosystem services for long-term sustainability purposes, including (1) focusing on ecosystem service provision over longer time frames [34,35], (2) understanding how current decisions will affect future ecosystem service provisions [36,37], and (3) ensuring dynamic planning for sustainable use of ecosystem services [38]. They also suggest that decision-makers should improve their understanding of indicators and their effective use to enhance the decision-making of the sustainable use of ecosystem services. Mapping multiple ecosystem services to understand associations among them has been one of the focuses of ecosystem service studies, which can be an effective tool to meet the demands of the above-suggested directions [33,39], given the geospatial nature of ecosystem services. It, however, has not been further studied in Bennett and Chaplin-Kramer [33].
Spatial zoning of ecosystem services is proposed in this study as a solution to meet the demands of the suggested directions above, as (1) an effective ESV zoning map facilitates a more sustainable long-term use of the ecosystem services in an area, (2) mapping and understanding the current zonal ESV characteristics provides a basis for understanding future ESV distributions, (3) an ESV zoning map can be dynamically updated, (4) a spatial zoning task takes various indicators into account and explains the associations among the indicators.
There exist various approaches for conducting spatial zoning. For example, Xu et al. [40] presented a zoning approach to account for the correlation between the changes in ecosystem service supply and demand, based on hotspot analysis and geographic detector techniques. Abdollahi et al. [41] adopted a spatial multi-criteria evaluation model to perform zoning related to river pollution. Machine learning methods have also been used for spatial zoning. For instance, Yang et al. [42] proposed a supervised classification algorithm to perform climate zoning. Additionally, Liu et al. [43] compared six common methods that can be used for landscape zoning based on ecosystem service bundles, including K-means++, SOM, Iterative Self-organizing Data Analysis, Expectation Maximization, Density-Based Spatial Clustering of Applications with Noise, and Natural Breaks (Jenks); the results show that the performances of the algorithms vary from area to area, while SOM mostly stayed at the top tier in terms of the four landscape pattern indices used (i.e., landscape division index, aggregation index, Shannon's diversity index, and perimeter-area fractal dimension index).
SOM is a machine learning approach for unsupervised clustering [23,24]. Unlike a multi-criteria evaluation model and supervised machine learning, prior knowledge is not required by SOM, and thus an ESV zoning map derived using SOM can be readily updated, which specifically meets the demand of the third direction mentioned above (i.e., ensuring dynamic planning for sustainable use of ecosystem services). Pérez-Hoyos et al. [44] used SOM to identify ecosystem functional types in Spain and indeed highlighted the potential of this approach for monitoring the spatial pattern of the ecosystem service status. Therefore, this study leverages SOM for the ecosystem service zoning of the GBA.

Study Area
The GBA ( Figure 1) is of monsoon climate mainly, crossing the tropics and the subtropics. Its annual temperature is around 22.3 degrees Celsius on average and annual precipitation ranges from 1300 to 2500 mm on average. It is abundant in illumination and water resources. The vegetation types in the GBA are diverse; the main type is the evergreen broad-leaved forest. The elevation is higher in the northern part (the mountainous side) while lower in the southern part (the coastal side).

Ecosystem Services
In this research, the four types of ecosystem services of China together with the corresponding 11 subtypes introduced in Xie et al. [18] were adopted (Table 1) for the exploratory analysis. There are 11 cities in the GBA, including two Special Administrative Regions (Hong Kong and Macao), two mega-cities (Guangzhou and Shenzhen), five large cities (Dongguan, Foshan, Jiangmen, Huizhou, and Zhaoqing), and two relatively small cities (Zhongshan and Zhuhai). Each city in the GBA has its economic characteristics. Guangzhou is located in the center of the GBA, which is a transportation hub and has both the heavy industry (e.g., auto manufacturing and petrochemical manufacturing) and the light industry (e.g., electronic product manufacturing), and the service industry. Shenzhen is one of the leading high-technology industrial centers in China (e.g., Huawei, Tencent, and ZTE). Hong Kong is an international hub for finance and famous for its wholesale and retail sales. Macao's economy is mainly characterized by its tourism and gaming industries. Dongguan and Foshan are major manufacturing centers in the GBA. Zhuhai mainly develops the light industry and tourism. The other cities in the GBA are relatively less developed. Jiangmen, Zhongshan, and Huizhou are dominated by the light industry. Zhaoqing's main sectors are the light industry and agriculture.

Ecosystem Services
In this research, the four types of ecosystem services of China together with the corresponding 11 subtypes introduced in Xie et al. [18] were adopted (Table 1) for the exploratory analysis. The overall ESV ( Figure 2) representing the total value of the 11 types of level-2 ecosystem services (Table 1) in the GBA was derived from Xu [46]. The data for the 11 types of level-2 ecosystem services provided in Xu [46] are 11 raster files in 1 km cell size (Supplementary Materials Figures S1-S11). Spatial variations in the overall ESV across the GBA can be seen in Figure 2, with generally lower values in the central parts of the GBA (i.e., the more urbanized areas). Spatial variations in the overall ESV across the administrative regions of the individual cities can also be seen ( Figure 2).
The sum of the ESV of the GBA is approximately 109.93 billion United States dollar (USD) and the overall ESV ranked by cities is shown in Table 2. Table 2 also lists the 2018 Gross Domestic Product (GDP) values of the GBA cities. The GDP data were collected from the Guangdong Statistical Yearbook [47], the World Bank Open Data [48], the Census and Statistics Department of Hong Kong [49], and the Worldometer [50]. Those less economically developed cities (e.g., Zhaoqing, Huizhou, and Jiangmen) are associated with greater overall ESV and overall ESV per capita, while those more economically developed cities (e.g., Shenzhen, Zhuhai, and Hong Kong) have lower overall ESV and overall ESV per capita.
The overall ESV (Figure 2) representing the total value of the 11 types of level-2 ecosystem services (Table 1) in the GBA was derived from Xu [46]. The data for the 11 types of level-2 ecosystem services provided in Xu [46] are 11 raster files in 1 km cell size (Supplementary Materials Figure S1 to Figure S11). Spatial variations in the overall ESV across the GBA can be seen in Figure 2, with generally lower values in the central parts of the GBA (i.e., the more urbanized areas). Spatial variations in the overall ESV across the administrative regions of the individual cities can also be seen ( Figure 2).

Figure 2.
The overall ESV of the GBA (spatial resolution: 1 km; temporal stamp: 2018). Note: quantile classification is adopted to classify the cell values, which leads to an effective visualization of the spatial variation in the data.
The sum of the ESV of the GBA is approximately 109.93 billion United States dollar (USD) and the overall ESV ranked by cities is shown in Table 2. Table 2 also lists the 2018 Gross Domestic Product (GDP) values of the GBA cities. The GDP data were collected from the Guangdong Statistical Yearbook [47], the World Bank Open Data [48], the Census and Statistics Department of Hong Kong [49], and the Worldometer [50]. Those less economically developed cities (e.g., Zhaoqing, Huizhou, and Jiangmen) are associated with greater overall ESV and overall ESV per capita, while those more economically developed cities (e.g., Shenzhen, Zhuhai, and Hong Kong) have lower overall ESV and overall ESV per capita.   In addition, the land use/land cover (LULC) of the GBA is illustrated in Figure 3; the LULC data were obtained from the Geographical Information Monitoring Cloud Platform of China [51]. The distribution of the overall and mean ESV in the GBA across different LULC types is shown in Table 3. The three LULC types with the greatest overall ESV are forest land, agricultural land, and water body; the three LULC types with the greatest mean ESV are water body, agricultural land, and forest land.
In addition, the land use/land cover (LULC) of the GBA is illustrated in Figure 3; the LULC data were obtained from the Geographical Information Monitoring Cloud Platform of China [51]. The distribution of the overall and mean ESV in the GBA across different LULC types is shown in Table 3. The three LULC types with the greatest overall ESV are forest land, agricultural land, and water body; the three LULC types with the greatest mean ESV are water body, agricultural land, and forest land.

ESV Zoning
The workflow of the ESV zoning includes three steps in general.
Step one is to extract the cell values from the 11 level-2 ecosystem service layers (Table 1 and Supplementary Materials Figures S1-S11) and generate an 11 indicators × 52,680 cells matrix as input of the SOM.
Step two is the SOM training and clustering process, in which the SOM size and parameters need to be determined.
Step three is to analyze the zonal characteristic in relation to two external variables (i.e., LULC and night-time lights). Details are provided in Sections 3.3.1 and 3.3.2 as follows.

SOM Zoning
This study adopted SOM to perform spatial zoning of the GBA based on the 11 level-2 types of ecosystem services (Table 1). SOM is a machine learning approach for unsupervised clustering [23,24]. It extracts and learns information from a multi-dimensional input space and generates a two-dimensional (2D) output space in which data properties are represented and connected topologically (Figure 4).
A SOM visualization consists of multiple neurons; each neuron has a unique position on the 2D output space (i.e., the SOM neuron grid or lattice) and associated data samples (a value array at a specific cell location of the indicator layers) from the input space ( Figure 4). Each data sample is mapped (linked) to a neuron (one neuron can involve multiple input data samples). Similar input data samples (determined according to the values of input indicators) are placed close together on the 2D output space of the neuron grid.
(t), the vector weights of the BMU and the neighboring neurons are updated using the following equation: where is an n-dimensional weight vector, ( ) is the learning rate, ℎ ( ) is the neighborhood kernel of the BMU, and v is a randomly chosen input vector from the training dataset. Specifically, assume that X = (x 1 , x 2 , . . . x n )' denotes the input indicators and that w l = (w l1 , w l2 , . . . w ln ) represents the weight vector associated with neuron l (each neuron is assigned a random vector of weights with the same dimensionality as the input data), n denotes the number of input indicators, and w lj is the weight assigned to input indicator x j linked to the neuron l, with each data sample of the training data randomly placed into the SOM. Data samples are iteratively introduced to the low-dimensional neuron grid to detect the best matching unit (BMU). The BMU is the "winning" neuron containing the weight vector that minimizes the Euclidean distance between itself and the data sample at hand. The iterative process is described as SOM training. That is, the principle of SOM is to identify the nearest node to each data sample, and then move the winning neuron (i.e., the BMU) closer to the training data sample at a certain learning rate. In each iteration (t), the vector weights of the BMU and the neighboring neurons are updated using the following equation: where w l is an n-dimensional weight vector, α(t) is the learning rate, h cl (t) is the neighborhood kernel of the BMU, and v is a randomly chosen input vector from the training dataset. The 2D output space of SOM can be further projected onto the real world as spatial zones based on the input data sample locations if the data samples are location-based ( Figure 4). As such, SOM serves as a valuable dimension reduction technique that can be used to generate a 2D zonal map that is more readable and comprehensible by decision and policymakers compared with the multi-dimensional input space.
The SOM of this study was performed using the Kohonen package of R [52]. The value vectors of the 11 ecosystem services associated with the 52,680 raster cells across the study area (1 km spatial resolution) were normalized and adopted as the input layers (indicators) of the SOM. For the SOM training, a hexagonal SOM structure was used, using a standard linearly declining learning rate from 0.05 to 0.01 over 300 iterations, with the whole dataset involved in every iteration [32,53]. Determining the size of the (number of neurons) SOM grid can be a qualitative process, depending on the necessary number of zones (i.e., the necessary degree of dimension reduction of the input space) required for policy-making. There is no perfect number of zones, considering the high spatial diversity and heterogeneity of natural resources. For example, a 3 × 3 neuron combination allows for at most nine zones; generating a relatively small number of zones allows policymakers to make relatively coarse ESV management plans by neglecting certain zonal details. A 7 × 7 neuron combination allows for at most 49 zones; generating a relatively large number of zones allows policymakers to make more detailed plans. However, it does not mean that the larger the number of neurons is, the more detailed classification it can result in, as it may end up with a large number of neurons with similar patterns (which are redundant and difficult to interpret without further clustering). Additionally, using a large number of neurons may result in many empty and useless neurons without any data point (data sample).
To minimize redundant neurons after a SOM has been trained, a hierarchical SOM clustering can be performed to further group similar neighboring output neurons based on Pearson correlations and a visual interpretation of the patterns of the neuron weight vectors; each resulting neuron cluster represents an ecosystem service zone of the GBA. Specifically, the Pearson correlation between a pair of neighboring output neurons being greater than 0.75 indicates that the two neurons are highly correlated and can be further clustered into the same zone [54]. Patterns in the distribution of the input indicators in the individual neurons can be observed by illustrating the weight vectors across the SOM, using a "fan diagram", where individual fan representations of the magnitude (importance) of each input indicator in the weight vector are exhibited for each neuron.

Zonal Characteristic Analysis
The ESV zones derived from the SOM were first characterized by their internal ESV compositions and structures; certain types of ecosystem services are the key indicators (dominant types) that differentiate a zone from the others. To further assess the differences and characteristics of the ESV zones derived using the SOM in relation to external variables, two analyses related to LULC and night-time lights were performed. These two factors pertain to the existing human activity situations on the ground and the accessibility to the GBA's ecosystem services. First, based on the 52,680 raster cells of the study area, for each ESV zone derived using the SOM, the LULC distributions across the ecosystem service zones are visualized using a stacked column chart to reveal the zonal LULC patterns. Second, based on the 52,680 cells of the study area, for each ESV zone derived using the SOM, the annual average radiance values of the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band Night-time Lights of the GBA ( Figure 5) across the ecosystem services zones are visualized using a box plot. The night-time dataset can reflect human activity levels [55]. The night-time data in this study were obtained from the Earth Observation Group of the National Oceanic and Atmospheric Administration of the United States of America [56].

Results
The workflow and parameters leading to the result are summarized in Figure 6. Detailed results are described in Sections 4.1 and 4.2 as follows.

Results
The workflow and parameters leading to the result are summarized in Figure 6.  Figure 7 shows the neuron weight vectors generated by the 5 × 5 SOM. It can be observed that neurons with similar weight vectors are close to each other on the neuron grid. Figure 8a demonstrates the number of data samples falling in each neuron of the SOM. Figure 8b presents the 11 ecosystem service zones derived based on the Pearson correlations and a visual interpretation of the patterns of the neuron weight vectors. As mentioned in Section 3.3.1, there is no perfect number of zones; the determination of the number of neurons depends on the degree of dimension reduction of the input space required by policymakers. Figure 9 illustrates the exploratory process toward the selection of a 5 × 5 combination. The 2 × 2 and 3 × 3 neuron combinations (Figures 9a,c) allow for at most four and nine zones, respectively, with which policymakers can make relatively coarse ESV management plans by neglecting certain zonal details. The 4 × 4 neuron combination (Figure 9e) allows for at most 16 zones (Figure 9c), with which policymakers can make more detailed plans. To possibly increase the level of zonal details, it was decided in this study to further explore a 5 × 5 SOM (Figure 7) and a 6 × 6 SOM (Figure 9g). However, it was noticed that the 6 × 6 combination resulted in empty and useless neurons without any data sample ( Figure 9h); therefore, the 5 × 5 combination was adopted eventually.

Ecosystem Service Zones and Their Characteristics
The 11 ecosystem service zones were further projected onto the study area, as shown in Figure 10 (the final zoning map). Zone 2, Zone 5, and Zone 8 appear to be the largest ecosystem service zones across the GBA. The characteristics of the 11 zones (described according to the fan diagram in Figure 7) and their ESV are provided in Table 4. Zone 2, Zone 6, and Zone 11 are the three zones with the greatest overall ESV; Zone 6, Zone 10, and Zone 11 are the three zones with the greatest mean ESV (Table 4). Additionally, Zone 1 and Zone 2 are similar but the latter is more capable of providing the SC supporting  Figure 7 shows the neuron weight vectors generated by the 5 × 5 SOM. It can be observed that neurons with similar weight vectors are close to each other on the neuron grid. Figure 8a demonstrates the number of data samples falling in each neuron of the SOM. Figure 8b presents the 11 ecosystem service zones derived based on the Pearson correlations and a visual interpretation of the patterns of the neuron weight vectors. As mentioned in Section 3.3.1, there is no perfect number of zones; the determination of the number of neurons depends on the degree of dimension reduction of the input space required by policymakers. Figure 9 illustrates the exploratory process toward the selection of a 5 × 5 combination. The 2 × 2 and 3 × 3 neuron combinations (Figure 9a,c) allow for at most four and nine zones, respectively, with which policymakers can make relatively coarse ESV management plans by neglecting certain zonal details. The 4 × 4 neuron combination (Figure 9e) allows for at most 16 zones (Figure 9c), with which policymakers can make more detailed plans. To possibly increase the level of zonal details, it was decided in this study to further explore a 5 × 5 SOM (Figure 7) and a 6 × 6 SOM (Figure 9g). However, it was noticed that the 6 × 6 combination resulted in empty and useless neurons without any data sample ( Figure 9h); therefore, the 5 × 5 combination was adopted eventually. vice than Zone 10 is. Zones 7 to 9 all provide provisioning services, regulating services, and supporting services, while they deviate from each other at the specific level-2 ecosystem services that they mainly provide (Tables 1 and 4). For instance, Zone 7 is more capable of providing the RMP and WS provisioning services than Zone 8 is; Zone 9 is more capable of providing the CR regulating service than Zone 7 and Zone 8 are. Lastly, Zone 11 is mainly characterized by the WS provisioning service and the HR regulating service (Tables 1 and 4).   vice than Zone 10 is. Zones 7 to 9 all provide provisioning services, regulating services, and supporting services, while they deviate from each other at the specific level-2 ecosystem services that they mainly provide (Tables 1 and 4). For instance, Zone 7 is more capable of providing the RMP and WS provisioning services than Zone 8 is; Zone 9 is more capable of providing the CR regulating service than Zone 7 and Zone 8 are. Lastly, Zone 11 is mainly characterized by the WS provisioning service and the HR regulating service (Tables 1 and 4).    The 11 ecosystem service zones were further projected onto the study area, as shown in Figure 10 (the final zoning map). Zone 2, Zone 5, and Zone 8 appear to be the largest ecosystem service zones across the GBA. The characteristics of the 11 zones (described according to the fan diagram in Figure 7) and their ESV are provided in Table 4. Zone 2, Zone 6, and Zone 11 are the three zones with the greatest overall ESV; Zone 6, Zone 10, and Zone 11 are the three zones with the greatest mean ESV (Table 4). Additionally, Zone 1 and Zone 2 are similar but the latter is more capable of providing the SC supporting service (Tables 1 and 4). Both Zone 3 and Zone 4 provide mainly supply and supporting services, while they deviate from each other at the main level-2 ecosystem services that they mainly provide, respectively (Tables 1 and 4). That is, comparing Zone 4 with Zone 3, Zone 4 is more capable of offering the FP and RMP provisioning services; Zone 3 is dominated by the SC supporting service while Zone 4 is dominated by the NCM supporting service. Zone 5 is singly characterized by the WS provisioning service (Tables 1 and 4). Zone 6 and Zone 10 are the most diverse in terms of the dominant ecosystem service types, which cover all four types of level-1 ecosystem services (Tables 1 and 4); however, Zone 6 is more capable of providing the WS provisioning service and HR regulating service than Zone 10 is. Zones 7 to 9 all provide provisioning services, regulating services, and supporting services, while they deviate from each other at the specific level-2 ecosystem services that they mainly provide (Tables 1 and 4). For instance, Zone 7 is more capable of providing the RMP and WS provisioning services than Zone 8 is; Zone 9 is more capable of providing the CR regulating service than Zone 7 and Zone 8 are. Lastly, Zone 11 is mainly characterized by the WS provisioning service and the HR regulating service (Tables 1 and 4).   Provisioning service (RMP), regulating service (GR and CR), and 21.22 1.55 Figure 10. The ecosystem service zoning result of the SOM projected onto the study area (spatial resolution: 1 km; temporal stamp: 2018). Table 4. Qualitative and quantitative characteristics of the ecosystem zones.

Zonal Characteristics Based on External Variables
The LULC patterns across the ecosystem service zones are shown in Figure 11. It is observed that Zones 1 to 3 and Zone 9 are all dominated by forest lands, while they differ from each other in terms of the percentages of the individual LULC types. For instance, Zone 1 has a greater percentage of agricultural lands; Zone 2 has the greatest proportion of forest lands and a tiny percentage of water bodies; Zone 3 has a greater proportion of urban or built-up lands; Zone 9 is high in both the percentage of urban or built-up land and the percentage of agricultural land. In addition, Zone 4, Zone 7, and Zone 8 are mainly characterized by agricultural lands and forest lands. Zone 5 is mainly covered by urban or built-up lands, agricultural lands, and forest lands. Zone 6 is mainly covered by water bodies, agricultural lands, and forest lands. Zone 10 is mainly covered by grass lands, agricultural lands, and forest lands. Zone 11 is mainly described by water bodies, agricultural lands, and urban or built-up lands. These individual zonal LULC characteristics reflect the ecosystem services that the individual zones can provide (Table 4). For example, Zone 11 provides more WS and HR services, given its higher water body and agricultural land coverages; those zones good at providing FP services are generally associated with a large proportion of agricultural lands (e.g., Zones 4, 6, 7, and 8).
The distribution of the 2016 annual average radiance value night-time lights across the ecosystem service zones is shown in Figure 12, which is also distinguishable based on the patterns in the box plot. It is observed that Zone 5 and Zone 11 have relatively high levels of human activities (implying high accessibility), which are mainly covered by urban or built-up lands, grass lands, and agricultural lands ( Figure 11). Zones 1 and 2 have the lowest levels of human activities (suggesting relatively low accessibility), which are dominated by forest lands (Figure 11). Other zones are with moderate levels of human activities. the ecosystem service zones is shown in Figure 12, which is also distinguishable based on the patterns in the box plot. It is observed that Zone 5 and Zone 11 have relatively high levels of human activities (implying high accessibility), which are mainly covered by urban or built-up lands, grass lands, and agricultural lands ( Figure 11). Zones 1 and 2 have the lowest levels of human activities (suggesting relatively low accessibility), which are dominated by forest lands (Figure 11). Other zones are with moderate levels of human activities. Figure 11. LULC patterns across the ecosystem service zones. Figure 11. LULC patterns across the ecosystem service zones.

Use of the ESV Zoning Map
The GBA's economic well-being and prosperity are sustained by its natural capital, including its ecosystems that offer essential goods and services for human beings. Mankind highly relies on well-functioning ecosystems that lay the foundation for a constant flow of ecosystem services from nature to society. Mapping ecosystem services and their zonings is important to understanding how the services contribute to human wellbeing. Mapping ecosystem services and their zonings are also important to policy-makings that have an impact on bringing ecosystem services to practical applications including but not limited to sustainable management of natural resources, land use optimization, nature conservation and restoration, landscape planning, environmental protection, climate protection, and nature-based solutions [57]. Indeed, maps can be adopted to communicate complex geospatial information and are thus essential tools for practical applications. Maps are also useful for raising awareness about spatial ecosystem services supply and demand.
Deriving benefits from natural ecosystems must prioritize sustainability, meeting the United Nations' sustainable development goals [58]. Economy and environment are two important dimensions of sustainability [58], the long-term sustainable interaction of which represent a sustainable economy. Ensuring green development and ecosystem stability has also been strengthened in China's Outline Development Plan for the Guangdong-Hong Hong-Macao Greater Bay Area [59]. Therefore, more specifically, the 11 types

Use of the ESV Zoning Map
The GBA's economic well-being and prosperity are sustained by its natural capital, including its ecosystems that offer essential goods and services for human beings. Mankind highly relies on well-functioning ecosystems that lay the foundation for a constant flow of ecosystem services from nature to society. Mapping ecosystem services and their zonings is important to understanding how the services contribute to human wellbeing. Mapping ecosystem services and their zonings are also important to policy-makings that have an impact on bringing ecosystem services to practical applications including but not limited to sustainable management of natural resources, land use optimization, nature conservation and restoration, landscape planning, environmental protection, climate protection, and nature-based solutions [57]. Indeed, maps can be adopted to communicate complex geospatial information and are thus essential tools for practical applications. Maps are also useful for raising awareness about spatial ecosystem services supply and demand.
Deriving benefits from natural ecosystems must prioritize sustainability, meeting the United Nations' sustainable development goals [58]. Economy and environment are two important dimensions of sustainability [58], the long-term sustainable interaction of which represent a sustainable economy. Ensuring green development and ecosystem stability has also been strengthened in China's Outline Development Plan for the Guangdong-Hong Hong-Macao Greater Bay Area [59]. Therefore, more specifically, the 11 types of level-2 ecosystem services (Table 1) have been adopted in this study as input indicators for the spatial ecosystem service zoning of the GBA. Based on the SOM machine learning approach, this study has characterized and mapped the spatial distributions of the ESV, resulting in 11 zones (Figure 10). By deriving the ESV zones (Figure 10), this study has shed light on using SOM as a tool for facilitating the sustainable exploitation and utilization of the GBA ecosystem services.
Reasonable and effective utilization of the ecosystem services in the GBA should be based on its zonal characteristics rather than haphazard exploitations. For example, this study has revealed that Zhaoqing is mainly covered by Zones 1, 2, and 8; Jiangmen and Huizhou are mainly covered by Zones 1 and 2 as well as Zones 7 and 8; Zhongshan is shared by multiple zones, including mainly Zone 2, Zone 5, Zones 8 and 9, and Zone 11. Each zone has its unique characteristics in terms of its internal ESV composition and structure and its relationships with external variables (i.e., LULC and levels of human activities). With the zonal map, decision-makers will be able to make targeted plans of ESV derivation for the individual city areas based on the spatial distribution and characteristics of the ESV zones ( Figure 10). For instance, Zone 1 and Zone 2 are similar in terms of the supply and regulating services that they can provide; however, Zone 1 is associated with a slightly higher level of human activity, implying slightly higher accessibility and meaning that policymakers can prioritize Zone 1 to ensure more environment-friendly exploitation and utilization of the supply and regulating services, which can minimize the construction of extra manmade infrastructures and reduce costs. Take Zone 6 and Zone 10 as another example; both zones meet the need if the goal is to derive multiple and diverse types of ecosystem services from a single zone; however, policymakers may prioritize Zone 6 if grass land conservation is a key concern because Zone 6 is covered by a much smaller percentage of grass land compared with Zone 10. Rational exploitation and utilization of the ecosystem services in such manners will contribute to the smart economy and smart environment of the GBA, and thus the sustainable economy. The smart economy here is referred to as a diverse, healthy, and free economy that brings benefits to the region; the smart environment here is referred to as intelligent resource management concerning both the natural and built environment [3].

Recommendations
This study concurs with Vollmer et al. [60] that mapping and representing ecosystem services spatially is crucial to the effectiveness of regional spatial planning. This research has presented a workflow to utilize SOM machine learning to perform spatial zoning of ecosystem services. Such zoning results ( Figure 10) have practical implications for the GBA's sustainable development. However, as mentioned in Section 1, as the ecosystem services in the GBA are dynamic, we do not intend to provide fixed and long-lasting zoning results in this work, but only to explore and demonstrate the applicability of SOM for such a zoning purpose. The zoning results and maps must be updated and fine-tuned depending on future changes of the input ESV indicators and zoning requirements. It is worth noting that the method can be applied for ESV zoning in other regions (not limited to the GBA) and can be applied in any zoning analysis for facilitating conservation-related decisionmaking and policymaking. The advantage of SOM is that it allows for some flexibility in determining the useful number of zones (see Section 3.3.1) according to policymakers' dimension reduction needs. However, it is worth stressing that this can be a disadvantage in certain circumstances; for example, the method and workflow are not suitable for circumstances in which one single gold-standard zoning result is required, as the determination of the SOM size is qualitative and can be somewhat arbitrary. Therefore, the method and workflow presented in this study should mainly be used for supporting decision-making rather than used for deriving gold-standard zoning maps.
For future works, first, more types of ecosystem services should be considered and derived as there should be more than the 11 types of ecosystem services ( Table 1) that the GBA can provide. Second, although this study investigated the zoning results in connection with LULC and night-time lights variables, actual ecosystem service accessibility analyses based on a topologically accurate road network dataset can be more useful for facilitating the effective utilization of the ecosystem services (e.g., calculating distances from the raster cells to different types of roads (hierarchical classes), road density within a certain radius, or perhaps even potential accessibility). Additionally, this study has mainly investigated the economic and environmental dimensions of sustainability. It will also be necessary to further explore the social dimension of sustainability. How people and society interact with the natural ecosystems can also play a role in the zoning. Indeed, Bennett and Chaplin-Kramer [33] have recommended some directions for achieving more sustainable socioecological systems, including (1) understanding the implications of ecosystem services provision for human wellbeing [61,62], (2) understanding the role of various forms of human capital in the provision of ecosystem services for human wellbeing [33], and (3) understanding the function of socioecological systems in the provision of ecosystem services, their resilience, and distribution [63,64]. Sustainability needs to be co-developed with practitioners in a way that matters to people and to be dynamically improved over time. Furthermore, the spatial resolution of the input layers of this exploratory study is 1 km, which should be increased in future works to refine the zoning outcome. Corresponding to the high-resolution input layers, a more detailed classification of LULC types would be useful to make the results even more meaningful. Lastly, systematically comparing SOM or exploring SOM in combination with other cutting-edge spatial zoning approaches (either supervised or unsupervised, e.g., deep SOM learning) that are potentially suitable for the stated purposes will be necessary next steps.

Conclusions
Sustainability is one of the major challenges in the 21st century for humanity. Many researchers have criticized that the existing studies on ecosystem system services have not led to more significant contributions to decisions about sustainable development. Spatial zoning of ecosystem services is proposed in this study as a solution to meet the demands for sustainable use of ecosystem services. This study performed a workflow and an exploratory analysis using SOM for visualizing the spatial patterns of the ESV of the GBA. The zoning was performed based on 11 types of ecosystem services, resulting in 11 ecosystem service zones. It is recommended that reasonable and effective utilization of the ecosystem services in the GBA should be based on its zonal characteristics rather than haphazard exploitations, which can contribute to the sustainable economy and environment of the region. The applicability of SOM for the GBA ecosystem service zoning has been demonstrated in this study. However, it should be re-stressed that the method and workflow presented in this study should mainly be used for supporting decision-making rather than used for deriving gold-standard zoning maps.