Spatial patterns of sandy beaches in China and risk analysis of human infrastructure squeeze based on multi-source data and ensemble learning

Sandy beaches provide essential ecological and economic services, but their functions are increasingly threatened by human activities. Analyzing the spatial distribution of China's sandy beaches and the impacts of human activities offers valuable insights for coastal resource management and ecological protection. However, remote sensing technologies face challenges such as limited data sources and tidal influences, which affect recognition accuracy. Therefore, integrating multi-source remote sensing data and reducing the impact of tidal fluctuations to improve recognition accuracy remains a key challenge. This study proposes an innovative approach utilizing multi-source data and an ensemble learning model to identify sandy beaches in China (2016–2024). By integrating Sentinel-1/2 satellite data, terrain data, and nighttime light data, along with spectral, index, terrain, texture, and polarization features, sandy beaches were identified across multiple years, and the results were consolidated into a unique dataset to analyze spatial patterns and risks from human infrastructure squeeze. (1) High-precision classification identified 3347 sandy beaches in China, covering a total area of 320.50 km². Guangdong had the largest number, area, and perimeter, while Hebei had the widest sandy beaches. (2) In Fujian, Guangdong, and Taiwan, the identified sandy beaches covered 54.57, 78.88, and 46.60 km², with perimeters of 1435.89, 2849.39, and 1324.98 km, and widths of 54.91, 38.92, and 57.17 m, respectively. These results were significantly better than those from published datasets. (3) From 1990 to 2024, the area at risk from human infrastructure squeeze increased from 134.39 to 181.42 km², a rise of 47.03 km², with the most significant increase occurring between 1995 and 2000. Guangdong and Fujian showed growth rates of 0.38 and 0.32 km² yr⁻¹, respectively. This study provides an up-to-date dataset on China's sandy beaches. It assesses their spatial patterns and human impact risks, contributing to research and policy for the sustainable development of coastal zones (https://doi.org/10.5281/zenodo.15307240, Meng et al., 2025).

1 Introduction

As transitional zones between land and sea, sandy beaches play a critical ecological role, providing essential habitats for various flora and fauna and offering crucial ecosystem services such as climate regulation, wave attenuation, and storm surge protection (Temmerman et al., 2013; Lansu et al., 2024; Mentaschi et al., 2018). Moreover, sandy beaches support a wide range of economic activities including tourism, fisheries, marine resource development, and land use, thereby holding considerable economic value (Mejjad et al., 2022; Brščić and Šugar, 2020; Ranasinghe, 2016). However, with the development of coastal areas, sandy beach ecosystems are increasingly threatened by human activities. Urban expansion, pollution, and tourism-related disturbances have transformed many natural sandy beaches into residential, commercial, or industrial land, leading to habitat degradation and the loss of ecological functions (Lu et al., 2018b; Barbier et al., 2011). Therefore, accurate identification and mapping of sandy beach areas are essential for ecological conservation, resource management, and sustainable coastal development.

Efficient and precise identification of sandy beach areas is a core component of sandy beach research (Ponte et al., 2016). Traditional field-based surveys, while offering high accuracy, are often constrained by high labor costs, low efficiency, and environmental limitations (Rodríguez-Padilla et al., 2024; Valderrama-Landeros et al., 2022; Bai et al., 2021; Blenkinsopp et al., 2021). In contrast, remote sensing technologies based on satellites, and aerial platforms – including unmanned aerial vehicles (UAVs), with their high temporal resolution and broad spatial coverage, enable effective data acquisition across different regions and climatic conditions and thus have become the primary method for large-area sandy beach identification (Xu et al., 2024; Sun et al., 2024; Salameh et al., 2019; Nativí-Merchán et al., 2021; Splinter et al., 2018). Current remote sensing-based sandy beach identification methods are generally categorized into visual interpretation and automated classification (Sekovski et al., 2014; Toure et al., 2019). While visual interpretation relies heavily on expert experience and high-quality imagery, it is susceptible to tidal variations and subjectivity, making it unsuitable for large-area or long-term monitoring. Automated classification methods, on the other hand, extract features such as spectral and texture information to achieve more consistent and efficient identification, gradually becoming the mainstream approach in remote sensing-based sandy beach studies. These methods are particularly suited for large-area monitoring, providing efficient and accurate results over extensive areas (Masria, 2024; Yasir et al., 2020; Wei et al., 2021).

Currently, remote sensing-based automated sandy beach recognition methods encompass a variety of techniques, including threshold segmentation, machine learning methods, and deep learning methods (Sekar et al., 2022; Vicens-Miquel et al., 2022; Rabehi et al., 2023). Threshold segmentation divides the image into sandy beach and non-sandy beach areas by analyzing pixel value thresholds of remote sensing images. This method is simple and efficient, but the segmentation performance is often poor in complex backgrounds or areas with subtle grayscale differences (Bao et al., 2018). Machine learning and deep learning methods are similar in principle, as both build classification models through training datasets, enabling them to handle complex images and exhibit strong pattern recognition capabilities. Machine learning methods typically require integrating multiple features from remote sensing images, such as spectral, spatial, and texture information, while deep learning methods automatically extract features and perform classification using multi-layer neural networks, capable of handling large-scale and complex datasets. For example, an exploratory survey and monitoring of sandy beaches was carried out by comparing Sentinel-2 and Landsat images, using random forests and various spectral indices (Latella et al., 2021), while a binary image segmentation method based on the U-Net model in convolutional neural networks was applied to accurately delineate the sandy beach outline of the southeastern coast of Australia (Yong et al., 2024). Despite significant progress in improving recognition accuracy and spatial scalability, machine learning and deep learning methods still face several challenges under complex conditions and tidal interference. Specifically, due to the single spectral feature of optical imagery, the recognition accuracy of machine learning methods still needs improvement in certain situations, particularly in areas where feature boundaries are unclear or the background is complex. For deep learning methods, the influence of tidal variations remains a significant challenge. To overcome the impact of tidal variations on sandy beach extraction from remote sensing images, multi-year data needs to be introduced, which requires a large amount of labeled data for training. However, high-quality labeled data are often difficult to obtain, and the training process is time-consuming, increasing computational resource consumption. Moreover, during the training of multi-layer neural networks, overfitting may also become an issue (Zhu et al., 2023; Chen et al., 2022; Zhou et al., 2023b; Wang et al., 2024b).

The coastal regions of China are among the most economically developed and dynamic areas in the country, encompassing several key economic belts and contributing to over 57 % of the national gross income (Tian et al., 2016; Wang et al., 2024; Feng et al., 2025). These areas not only have a solid industrial foundation, vibrant trade, and advanced service sectors but also serve as the core engines of national economic growth. However, despite the rapid economic development in these regions, there is a lack of nationwide dynamic monitoring tools for sandy beaches, and the risks posed by human infrastructure squeeze – particularly due to urban development and coastal expansion – are not well understood. Past research has mainly relied on single optical imagery data and has not fully considered the use of multi-year data to address tidal changes, which limits the comprehensive understanding and monitoring of coastal beach changes. To address these challenges, this study is divided into two main components: (1) This study proposes a novel method for sandy beach identification using multi-source remote sensing and ensemble learning. Sentinel-1/2 imagery, nighttime light, and terrain data are integrated to extract spectral, index, terrain, polarization, and texture features. A stacking ensemble framework enhances classification accuracy and robustness. Multi-year sandy beach data extracted using ensemble learning from 2016 to 2024 are merged to construct an annual representative sandy beach dataset, reducing tidal fluctuation impacts and supporting precise mapping and monitoring; (2) Based on human infrastructure data, buffer zones are created and combined with this dataset to identify areas at risk of encroachment. Spatial analysis evaluates anthropogenic pressure on beaches across regions. The results offer a scientific basis for delineating ecological red lines and inform strategies for beach conservation and coastal zone management.

2 Study area and materials

2.1 Study area

To identify sandy beaches along China's coast from 2016 to 2024, this study adopts the 2024 coastline of China as the baseline. A buffer zone extending 10 km inland and 20 km offshore is applied to ensure the integrity of the sandy beach area (Li and Niu, 2022). The study area spans from 107°54^′ to 124°21^′ E longitude and from 17°58^′ to 41°8^′ N latitude, encompassing 14 regions (Liaoning, Hebei, Shandong, Jiangsu, Zhejiang, Fujian, Guangdong, Hainan, Taiwan, Guangxi, Tianjin, Shanghai, Hong Kong, and Macao) (Fig. 1). For statistical convenience, the coastal regions of Hong Kong and Macao are included within the Guangdong coastal region.

Figure 1 Location of China's coastal zone and distribution of partial typical land feature points for results verification.

2.2 Remote sensing data

This study used Sentinel-1, Sentinel-2, and Google Earth imagery, with data collected from 1 January 2016 to 31 December 2024. The Sentinel-1 images are Level-1 Ground Range Detected (GRD) products, acquired in Interferometric Wide (IW) mode, including dual-polarization VV (vertical transmit, vertical receive) and VH (vertical transmit, horizontal receive) data (Tian et al., 2022). A total of 23 859 Sentinel-1A images were selected, and annual mean composites were generated for each polarization channel (Figs. 2a and S1 in the Supplement). For Sentinel-2, five spectral bands were used: four visible and near-infrared bands (Bands 2, 3, 4, and 8) with a resolution of 10 m, and one shortwave infrared band (Band 11) with a resolution of 20 m. A total of 80 759 L2A and L1C images were selected, with cloud cover below 10 % (Figs. 2b and S1) (Li et al., 2022). After atmospheric correction, the QA60 quality band was applied to remove clouds and cloud shadows. Since atmospheric correction is not essential for applications not involving complex spectral analysis, L1C images were used for years lacking L2A data (Wang et al., 2023). Finally, cloud-free median composites were generated and resampled to a 10 m resolution. Google Earth imagery was sourced from the “Google Earth Engine” plugin and online map services on the QGIS platform, with a spatial resolution of 1 to 5 m, enabling clear visualization of beach areas with high geographical accuracy and temporal consistency.

Figure 2 Spatial distribution of Sentinel-1 and Sentinel-2 images used in this study from 2016 to 2024: (a) Image Count of Sentinel-1, (b) Image Count of Sentinel-2.

2.3 Other data

The NASADEM dataset, provided by NASA and based on C-band radar with a 30 m resolution, accurately reflects global elevation. This study used the 2000 NASADEM data, resampled to 10 m, as terrain data for the period from 2016 to 2024. VIIRS nighttime light data (VNTL), provided by NOAA and NASA with a 500 m resolution, were used to represent global nighttime illumination. The data from 2016 to 2024 were selected, resampled to 10 m, clipped to [0,63] to remove outliers, normalized, and then used as the nighttime light dataset for this study. The land cover data selected for this study is the Landsat-derived annual China Land Cover Dataset (CLCD) from 1990 to 2024, with a resolution of 30 m, used to study changes in the area of human infrastructure squeeze risk (Yang and Huang, 2021). The per capita GDP, resident population, built-up area, and road area of the study area from 1990 to 2024 were obtained from the National Bureau of Statistics and local statistical bureaus (https://www.stats.gov.cn/, last access: 15 November 2025). In addition, data from 1990, 2000, 2010, and 2020 were extracted for Jinjiang, Licheng, and Fengze Districts of Quanzhou City to assess human infrastructure squeeze risk. GDP and population density data were obtained from the China GDP and Population Spatial Distribution Kilometer Grid Datasets (Xu, 2017a; Xu, 2017b). Built-up area data were sourced from the GHS-BUILT-H R2023A dataset published by the European Commission Joint Research Centre in 2022–2024 (Pesaresi and Politis, 2023). Road network data were derived from the national basic geographic information data (https://www.webmap.cn/, last access: 8 April 2025) and Harvard University (Baum-Snow et al., 2017). If data for a specific year were unavailable, the data from the nearest available year were used instead.

2.4 Published datasets

In this study, we used three datasets to evaluate our identified sandy beach dataset (Table 1): (1) The China sandy beach dataset, data directly obtained from the OpenStreetMap (OSM) database; (2) The 2022 China 10 m sandy beach dataset was identified using a support vector machine method based on Sentinel-2 imagery (Ni et al., 2024); (3) The 2020 China coastal land use dataset at 10 m resolution was identified using an object-oriented classification method based on Sentinel-2 imagery (Miao et al., 2022).

Table 1 Datasets used to evaluate the accuracy of sandy beach classification.

3 Methodology

To accurately monitor the current status of sandy beaches in China, this study integrates multi-source data from 2016 to 2024, including Sentinel-1/2 imagery, terrain data, and nighttime light data. By processing these multi-source datasets, spectral, index, terrain, texture, and polarization features were extracted to construct a multi-source, multi-feature, and multi-model framework based on a stacking ensemble learning approach, aiming to delineate the spatial distribution of sandy beaches in China. On this basis, 100 m buffer zones were established based on impervious surface data from 1990 to 2024 and combined with the sandy beach dataset for overlay analysis to identify areas at risk of human infrastructure squeeze and to analyze their spatiotemporal dynamics. The technical framework of the study is illustrated (Fig. 3).

Figure 3 The technical framework of the study.

3.1 Sample generation via visual interpretation

The main land cover types in the non-beach areas of the study area include farmland, forest, water, and impervious, which are defined based on typical coastal land cover characteristics in China (Yang et al., 2021). In the sample library, these non-sandy beach types, in addition to the sandy beach samples, play an important role. Therefore, we define the land use types in the study area as sandy beaches and non-sandy beaches, with non-sandy beach types including farmland, forest, water, and impervious. Based on field surveys, we conducted annual visual interpretation using Google Earth and Sentinel-2 imagery to label samples from 2016 to 2024. Through this approach, a candidate sample library was built, containing 341 057 samples, with 249 098 for model training and 91 959 for testing (Figs. 4 and S2 in the Supplement). The samples were randomly divided into training and testing sets at an approximate 7:3 ratio, while ensuring consistent sample proportions across the 12 regions of the study area and maintaining balanced ratios of sandy beach and non-sandy beach samples within each region. This partitioning strategy effectively avoided regional bias and ensured the reliability and representativeness of the model evaluation.

Figure 4 Annual counts of sandy beach and non-sandy beach samples in the training and testing sets.

3.2 Constructing sandy beach classification input features using multi-source data

This study selected 24 indicators from five dimensions – spectral, index, terrain, polarization, and texture – as input features for the ensemble learning model (Table 2). First, the VV and VH polarization modes were selected as polarization features. Next, five bands – Band 2, Band 3, Band 4, Band 8, and Band 11 – were chosen as spectral features, while eight indices – NDWI, NDVI, NDSI, NDBI, EVI, SAVI, BSI, and NDUI – were treated as index features. Then, grayscale images were generated using Band 3, Band 4, and Band 8, and six texture metrics – contrast, correlation, variance, angular second moment, inverse difference moment, and entropy – were derived from the grayscale images as texture features. Finally, elevation, slope, and aspect were selected as terrain features.

Table 2 Input features for sandy beach classification.

3.3 Sandy beach classification using an ensemble learning model

Stacking is a powerful ensemble learning method that uses predictions from multiple base learners as inputs to a meta-learner for final prediction. It combines the strengths of different models to overcome individual limitations, enhancing accuracy, stability, and generalization (Chen et al., 2024). In this study, Random Forest (RF), Gradient Boosting Decision Tree (GBDT), eXtreme Gradient Boosting (XGB), and Light Gradient Boosting Machine (LGBM) were selected as base learners. Their output classification probabilities were used as input features for the meta-learner, which adopted Logistic Regression (LR) to integrate probabilities, calculate distances to target classes, and produce final results (Table 3). For performance comparison, Support Vector Machine (SVM), K-Nearest Neighbor (KNN), RF, Classification and Regression Tree (CART), and Convolutional Neural Network (CNN) models were also employed as benchmark models.

Table 3 Parameters of models.

3.4 Model results processing and analysis

3.4.1 Post-processing

Pixel-based classification algorithms inevitably produce salt-and-pepper noise because sandy beaches share similar spectral characteristics with certain areas, such as buildings and bare land, making them difficult to distinguish (Mattson et al., 2024). To reduce noise and remove isolated areas, we applied spatial filtering and connectivity analysis. A 100 m coastal buffer based on administrative boundaries was used to clip and filter non-sandy beach areas. Based on sandy beach patches, holes smaller than 2000 m² (20 pixels) were filled, and sandy beach patches smaller than 1000 m² (10 pixels) were removed. Finally, the sandy beach boundaries were refined through a thorough review and manual interpretation based on visual assessment.

3.4.2 Accuracy assessment

This study assessed the accuracy of sandy beach classification results using 91 959 test samples from 2016 to 2024, calculating Accuracy, Precision, Recall, sandy beach F1-score, and Area Under the Curve (AUC) through the confusion matrix (Table 4). These metrics were chosen to comprehensively evaluate different types of classification errors and overall model performance. Accuracy measures the proportion of correctly classified samples, reflecting overall correctness. Precision indicates the reliability of positive predictions, highlighting false positives, while Recall reflects the model's ability to detect all positive samples, highlighting false negatives. The F1-score balances Precision and Recall, providing a single metric when both false positives and false negatives are important. AUC evaluates the model's discrimination ability across different thresholds, indicating stability and robustness.

Table 4 Accuracy assessment metrics.

TP (True Positive) refers to the number of samples that are actually positive and correctly classified by the model, while TN (True Negative) refers to the number of samples that are actually negative and correctly classified. FP (False Positive) represents the number of samples that are actually negative but incorrectly classified as positive, and FN (False Negative) represents the number of samples that are actually positive but incorrectly classified as negative. TPR (True Positive Rate, or Recall) measures the model's ability to identify positive samples and is calculated as $\text{TPR}=\frac{\text{TP}}{\text{TP}+\text{FN}}$ whereas FPR (False Positive Rate) measures the proportion of negative samples incorrectly classified as positive, calculated as $\text{TPR}=\frac{\text{FP}}{\text{FP}+\text{TN}}$.

3.4.3 Analysis of model variable importance

SHAP (Shapley Additive exPlanations) values are a model interpretation method based on cooperative game theory, used to quantify the contribution of each feature to the model's prediction. They calculate the marginal contribution of each feature across all possible feature combinations and take the average, thereby fairly attributing the “responsibility” for the prediction. This approach can explain the prediction behavior for individual samples or the entire model. In this study, SHAP values were used to interpret the contributions of a stacked ensemble model. Specifically, SHAP values were first calculated for each base learner (RF, GBDT, XGB, LGBM), then these values were normalized and weighted according to the coefficients of the meta-learner (LR) to obtain the contribution of each base learner to the final ensemble prediction. Finally, the weighted SHAP values of all base learners were summed to derive the SHAP values for the entire stacked ensemble model, thereby quantifying the influence of each base model on the final prediction and providing a measure of interpretability for the ensemble learning.

3.4.4 Analysis of human infrastructure squeeze risk trends

Approximately 33 % of global sandy beaches have less than 100 m of infrastructure-free buffer space, and in this study, human infrastructure squeeze risk is defined as the risk that occurs when human infrastructure (such as buildings, roads, ports, etc.) is located within 100 m of a sandy beach area. This threshold was chosen because infrastructure within this range can significantly impact the sandy beach area, leading to changes in its size or degradation of the ecosystem (Lansu et al., 2024). Therefore, this study established a 100 m buffer zone around the Chinese land use dataset and conducted an overlay analysis using the obtained sandy beach data. The analysis generated annual human infrastructure squeeze risk areas to evaluate the squeeze effects of infrastructure on sandy beaches from 1990 to 2024. To further quantify the spatiotemporal changes in human infrastructure squeeze risk areas across regions, the study applied the Sen's slope method to calculate trend values (Zhang et al., 2023b; Qiao and Xia, 2025). By applying this method, the study quantified the annual change rate and long-term trend of the risk zones across different regions from 1990 to 2024.

4 Results and analysis

4.1 Accuracy Assessment and Mapping of Chinese Sandy Beaches

The ensemble learning algorithm's performance was evaluated using key metrics such as Accuracy, Precision, Recall, sandy beach F1-score, and AUC, with the results presented in Table 5. From 2016 to 2024, all classification results had values above 0.85 for these metrics, with the best results in 2024, confirming the algorithm's stability, reliability, and accuracy in sandy beach classification.

Table 5 Accuracy assessment of ensemble learning results.

Based on the maximum-extent composite of sandy beach recognition results over 8 years, the spatial distribution of the number, length, width, and area of sandy beaches in China is shown in Fig. 5. Nationwide, a total of 3347 sandy beaches were identified, with a total length of 10 343.57 km, an average width of 46.42 m, and a total area of approximately 320.50 km² (Fig. 5).

Figure 5 Distribution of sandy beaches in China and the morphology of typical coastal areas. The maps/aerial images shown are Sentinel-2 imagery obtained from the public GEE platform (© Google Earth Engine).

From a provincial perspective, Guangdong hosts the largest number of sandy beaches, with 1096 sandy beaches, and also ranks first in both total length (2849.39 km) and total area (78.88 km²). In terms of sandy beach width, Hebei has the widest sandy beaches, with an average width of 68.77 m. Other regions, such as Fujian (sandy beach area 54.57 km², total length 1435.89 km, width 54.91 m) and Hainan (sandy beach area 51.65 km², total length 1977.96 km, width 41.02 m), also show significant sandy beach resources. This highlights the notable spatial differences in coastal sandy beach resources across China. The results indicate that coastal sandy beach resources in China exhibit significant uneven distribution, with considerable differences in the number, shape, and scale of sandy beaches across different regions (Fig. 6).

Figure 6 Spatial distribution of sandy beaches in China: (a) sandy beach numbers, (b) sandy beach length, (c) sandy beach width, and (d) sandy beach area. The histograms in panels (b), (c), and (d) are based on the results of each individual beach rather than regional statistics.

4.2 Feature Importance Analysis, Model Evaluation, and Dataset Comparison

We evaluated the importance of model variables for the years 2016–2024 (Fig. 7). Based on the average importance ranking from 2016 to 2024, the top five features were $\text{Elevation}\left(\mathrm{0.1270}\right)>\text{B3}\left(\mathrm{0.1242}\right)>\text{B2}\left(\mathrm{0.1049}\right)>\text{NDVI}\left(\mathrm{0.0783}\right)>\text{VV}\left(\mathrm{0.0600}\right)$. When classified by feature category, the highest average importance was found in spectral features (0.3403), followed by index features (0.2458), topographic features (0.1847), texture features (0.1317), and polarization features (0.0975). The most important feature in each category was as follows: Elevation for topographic features, B3 for spectral features, NDVI for index features, VV for polarization features, and VAR for texture features.

Figure 7 Analysis of model variable importance.

We evaluated the performance of our model and several other machine learning algorithms (SVM, KNN, RF, CART, and CNN) using Accuracy, Precision, Recall, F1-score, and AUC on the validation set (Table 6). Our model outperformed all others, achieving the highest accuracy (0.9335), precision (0.9160), and AUC (0.9802). It showed significantly better classification performance, especially in reducing misclassification. In comparison, the RF model had good precision but lower recall and accuracy, while KNN performed the worst across all metrics. To further assess the accuracy and reliability of the dataset, this study compared three published datasets in selected areas of Fujian, Shandong, and Guangdong (Fig. 8). Dataset 1, directly obtained from the OpenStreetMap (OSM) database, is highly subjective and tends to misclassify non-beach areas as sandy beaches, while also missing some actual sandy beach areas. Dataset 2, constructed using a support vector machine on Sentinel-2 imagery, has high consistency with our dataset, but still misses some sandy beach areas. Dataset 3, created using an object-oriented approach, demonstrates high accuracy for other land cover types but faces significant misclassification issues with bare land and urban areas. The results show that our dataset provides higher accuracy in sandy beach classification, significantly reducing misclassification.

Table 6 The accuracy of each model.

Figure 8 Comparative analysis of sandy beaches identified in this study and three published datasets. The maps/aerial images shown are Sentinel-2 imagery obtained from the public GEE platform (© Google Earth Engine).

According to the comparison results, our dataset shows significant advantages over published datasets 1, 2, and 3 in several key metrics, particularly in terms of sandy beach area, perimeter, width, and number (Fig. 9). Overall, our dataset provides a certain supplementary advantage in sandy beach area coverage, with larger areas in all regions: Fujian (54.57 km²), Guangdong (78.88 km²), and Taiwan (46.60 km²), significantly surpassing published datasets 1 (30.35, 20.59, and 19.07 km²) and 2 (29.17, 58.35, and 25.51 km²). Regarding perimeter, our dataset closely matches actual sandy beach boundaries as shown in Figs. 8 and 9: Fujian (1435.89 km), Guangdong (2849.39 km), and Taiwan (1324.98 km), compared to published datasets 1 (581.95, 826.40 and 509.22 km) and 2 (756.92, 1856.62, and 906.63 km). In terms of width, our dataset also outperforms published datasets 1 (52.21, 20.94 and 48.42 m) and 2 (45.18, 32.58 and 27.56 m), with values of 54.91 m in Fujian, 38.92 m in Guangdong, and 57.17 m in Taiwan. The number of identified sandy beaches in our dataset is higher than in published datasets 1 and 2, further highlighting the reduced misclassification and noise in our results (Fig. 9). Moreover, as shown in Figs. 8 and 9, published dataset 3 has far higher area, perimeter, width, and number values than the other three datasets and the actual values, leading to many non-sandy beaches being incorrectly identified as sandy beaches.

Figure 9 Provincial statistical comparison of sandy beaches in this study and three published datasets (regions sorted from south to north): (a) Length comparison, (b) Width comparison, (c) Area comparison, (d) Number comparison.

4.3 Spatial-temporal changes in human infrastructure squeeze risk

We calculated and analyzed the temporal changes in the risk areas affected by human infrastructure squeeze in the study area (Fig. 10). From 1990 to 2024, the risk area increased annually due to human infrastructure development (e.g., ports, dams, and coastal development). The area rose from 134.39 km² in 1990 to 181.42 km² in 2024, a total increase of about 47.03 km². The most significant increases occurred between 1995–2000, with a nationwide increase of approximately 19.34 km² in human infrastructure squeeze risk areas during this period. However, the growth rate has slowed in recent years, with only a 0.80 km² increase from 2020 to 2024, indicating stabilization due to enhanced regulation or resource saturation.

Figure 10 Interannual changes in human infrastructure squeeze on sandy beaches: (a) Risk area of human infrastructure squeeze from 1990 to 2024, (b) Changes in risk area of human infrastructure squeeze from 1990 to 2024.

We conducted a systematic assessment and comparative analysis of the current human infrastructure squeeze risk area in coastal regions of China in 2024, as well as its trend from 1990 to 2024 (Figs. 11 and S3 in the Supplement). In 2024, the distribution of human infrastructure squeeze risk areas showed significant regional differences. The top three regions with the largest risk areas were Guangdong (37.62 km²), Shandong (33.02 km²), Taiwan (32.52 km²), and Fujian (32.13 km²), all exceeding 30 km², accounting for a significant proportion of the national coastal squeeze risk area. These regions are densely populated, economically developed, and have high coastal development density, making them typical areas of coastal urbanization and industrial concentration. Hainan (16.42 km²) and Liaoning (14.18 km²) also showed high-risk areas, reflecting ongoing interference with natural sandy beaches due to port development and coastal infrastructure. In contrast, Tianjin (0.40 km²), Shanghai (0.40 km²), and Jiangsu (1.23 km²) had smaller risk areas, likely due to limited sandy beach resources. The trend analysis from 1990 to 2024 revealed significant growth rates in Guangdong (0.38 km² yr⁻¹) and Fujian (0.32 km² yr⁻¹), indicating rapid and substantial sandy beach squeeze risks.

Figure 11 Spatial changes of human infrastructure squeeze risk area: (a) current risk area, (b) changes in the risk area from 1990 to 2024.

5 Discussion

5.1 Selection and combination of input features

The influence of different input features and their combinations on recognition accuracy was analyzed (Table 7), taking 2024 as an example. Results show that spectral (S) features outperforms other single features, achieving the highest accuracy (Accuracy: 0.8976, F1-score: 0.8907). Index (I) features perform slightly worse, while terrain (T), texture (Tr), and polarization (P) features alone yield relatively low accuracies. Multi-feature fusion substantially improves model performance, with combinations involving spectral features showing the greatest gains. Among them, the “$\text{Spectral}+\text{Index}+\text{Texture}+\text{Terrain}$” combination achieves high accuracy (Accuracy: 0.9455, F1-score: 0.9416), while the five-feature combination (“$S+I+T+\text{Tr}+P$”) delivers the best overall results (Accuracy: 0.9476, F1-score: 0.9438, AUC: 0.9873). We compared the values of sandy beach and non-sandy beach areas for polarization features (VV), texture features (Correlation), and terrain features (Elevation) to explore the factors contributing to improved model accuracy (Fig. 12). In terms of polarization features (VV), the polarization reflectance of sandy beach areas is lower, while non-sandy beach areas, such as water bodies and buildings, show higher reflectance (Novelli et al., 2016). This may be due to the different scattering characteristics of sand compared to other surfaces, such as water, which is smooth and reflects stronger signals (Lu et al., 2018a). This contrast might be influenced by factors like data source and acquisition angle and requires further investigation. For texture features (Correlation), sandy beach areas exhibit consistent textures due to uniform sand grain arrangements, while non-sandy beach areas show more complex textures influenced by diverse land cover (Hasituya et al., 2016). In regions with complex land cover, texture differences may diminish, affecting recognition accuracy. Regarding terrain features (Elevation), sandy beaches are typically flat and near sea level, while non-sandy beaches exhibit more significant elevation differences (Li et al., 2023). In regions with minimal terrain variation, terrain features may have a weaker influence (Wang et al., 2024a).

Table 7 Variations in input features and their combinations.

Figure 12 Comparison of polarization, texture, and terrain features between sandy beaches and non-sandy beaches.

5.2 Effectiveness and robustness of ensemble learning

In this study, the Stacking model integrates four base classifiers – RF, GBDT, XGB, and LGBM – effectively combining their strengths in feature extraction and classification. This ensemble approach enhances the model's adaptability and generalization, improving its robustness in complex sandy beach area recognition tasks (Chen et al., 2021; Chrysafis et al., 2017). Each base model has its advantages: RF effectively handles high-dimensional features well and reduces overfitting, GBDT has strong nonlinear fitting capability, XGB achieves a good balance between optimization efficiency and accuracy, while LGBM performs with higher training speed and precision under large-scale and high-dimensional data conditions (Tian et al., 2025a; Mugunthan et al., 2023). By integrating these models, Stacking adapts better to challenges such as varying lighting, tidal conditions, and complex terrain, thus reducing misclassification. Although XGB demonstrates a slight advantage in Accuracy (0.9508), Precision (0.9433), and AUC (0.9894), and LGBM also shows certain strengths in Accuracy (0.9491), Precision (0.9414), and AUC (0.9890), the Stacking model outperforms all base learners in Recall (0.9549) and F1-score (0.9472), while exhibiting comparable performance to the base learners across the remaining metrics (Fig. 13), highlighting its robustness and reliability in overall performance. This demonstrates its stable performance under different data distributions. Overall, Stacking not only ensures competitive accuracy but also offers robust generalization for remote sensing tasks like sandy beaches change monitoring.

Figure 13 Comparison of F1-score and overall accuracy across models.

5.3 Factors shaping the spatial pattern of China's sandy beaches

This study produced a 10 m resolution sandy beach dataset for China's coastal regions. The spatial distribution of sandy beaches shows a clear pattern: more in the north and the south, and fewer in the central coast. This pattern is shaped by a combination of geological structures, marine dynamics, climate conditions, and river sediment input (Tian et al., 2025b; Ma et al., 2024; Wang and Wu, 2010). The northern coast features flat terrain and large rivers like the Yellow River and Liao River, which deliver abundant sediment, supporting sandy beach formation (Zhang et al., 2022). In the south, strong marine forces and tropical climates allow rivers to carry large amounts of sediment to the coast, where wide sandy beaches form under suitable conditions (Zhou et al., 2023a; Mentaschi et al., 2018). In contrast, the central coast is dominated by bedrock and muddy coastlines. Major rivers such as the Yangtze mainly contribute to tidal flats and muddy shores rather than typical sandy beaches (Kuai et al., 2021; Liu et al., 2020). Overall, the distribution of sandy beaches in China is shaped by sediment supply, coastal type, and hydrodynamic conditions, with more beaches in the north and south, and fewer along the central coast.

5.4 Exploration of the factors affecting the change in human infrastructure squeeze risk area

Based on the 1990–2024 Sen's slope trend analysis, several patterns emerge when analyzing the relationship between human infrastructure squeeze risk areas and related factors across regions (Fig. 14). Economically developed regions like Shanghai, Jiangsu, and Tianjin show slower growth in squeeze risk areas. This can be attributed to their smaller sandy beach areas, which limit the available space for human activities, and their higher per capita GDP, which is associated with more efficient urban planning and infrastructure development, leading to slower increases in squeeze risk. In contrast, regions with larger sandy beach areas, such as Guangdong, Shandong, Fujian, Taiwan, and Hainan, show faster growth in squeeze risk, as more space for development leads to higher risk over time. Additionally, regions with larger resident population, like Guangdong, Shandong, and Hebei, experience a positive correlation between population growth and squeeze risk. As the resident population increases, urbanization accelerates, putting more pressure on infrastructure and increasing the squeeze risk. Both built-up area and road area are positively correlated with squeeze risk, as their expansion reflects urban growth, intensifying competition for space.

Figure 14 Spatial changes of various factors from 1990 to 2024: (a) per capita GDP change, (b) resident population change, (c) built-up area change, (d) road area change.

We conducted a comprehensive analysis of the 1990, 2000, 2010, and 2020 data for the typical areas of human infrastructure squeeze risk in Jinjiang District, Quanzhou District, and Fengze District of Quanzhou, as detailed in the Supplement (Fig. S4). The analysis focused on the changes in indicators such as GDP, population density, built-up area, and road length. The results show that all these indicators have significantly increased over the past few decades, which is closely related to the risk of human infrastructure squeeze on sandy beaches. Specifically, with the growth of GDP, the expansion of economic activities has led to a continuous concentration of population in these areas, increasing the potential pressure on sandy beaches and the surrounding ecosystem. The expansion of built-up areas reflects the progress of urbanization, which, especially in vulnerable ecological regions such as sandy beaches, is often accompanied by a squeeze on natural resources. Although the increase in road length has promoted economic development, it has also intensified the risk of the squeeze on coastlines and surrounding sandy beaches.

5.5 Limitations of the proposed method and dataset

This study combines stacking ensemble learning models with remote sensing imagery to create a 10 m resolution sandy beach dataset, integrating multi-dimensional features such as spectral, terrain, polarization, and texture features. However, the study has some areas for improvement: (1) The sandy beach area is significantly affected by tidal changes, and the remote sensing images acquired at different tidal levels may exhibit considerable spatial differences. The current study has not fully addressed the issue of tidal phase matching. Future research should combine tidal prediction models or measured tidal data to select multi-temporal remote sensing images that are consistent with the target tidal phase, reducing the impact of tidal variations on beach extraction and improving the spatiotemporal accuracy of sandy beach extraction. (2) This study has preliminarily applied the Stacking ensemble learning method; however, there is still room for optimization in model selection, fusion strategy, and feature contribution evaluation. Future research can further optimize the base model structure, explore the synergistic effects between different models, and incorporate feature importance analysis to enhance the stability and classification accuracy of the ensemble model. (3) Currently, the assessment of human infrastructure squeeze risk in this study primarily relies on impervious surface expansion and buffer zone analysis. However, with relatively low spatial and temporal resolution, it is difficult to accurately depict the dynamic process of human infrastructure squeeze on sandy beaches. Future research should introduce higher spatial and temporal resolution data (e.g., high-frequency remote sensing monitoring, annual building change data, etc.) to enable continuous dynamic monitoring and analysis of the human infrastructure squeeze process.

6 Data availability

The dataset on the distribution of sandy beaches in China and the distribution of human infrastructure squeeze risk is available at: https://doi.org/10.5281/zenodo.15307240 (Meng et al., 2025).

7 Conclusions

This study used ensemble learning algorithms to identify sandy beach areas between 2016 and 2024 and validated the superiority of the method through various classification metrics. The results demonstrate that the method exhibits high stability and accuracy in the beach classification task, effectively distinguishing sandy beaches from other land types. Additionally, there are significant regional differences, with some regions having abundant sandy beach resources while others have relatively fewer sandy beach areas, indicating an uneven distribution of resources. By comparing with several published datasets, the proposed dataset in this study showed higher accuracy and reliability, especially in some typical coastal areas, ensuring precise classification of sandy beaches. This provides more accurate data support for sandy beach resource management and ecological protection, particularly in the context of increasingly scarce beach resources.

In addition to identifying sandy beach resources, this study also focused on the risk of infrastructure squeeze on sandy beach resources along China's coast. The study found that as infrastructure projects such as port construction and seawall building continue to advance, coastal areas are facing the gradual squeeze of their sandy beaches. Although the rate of expansion in some regions has slowed, the pressure from infrastructure development on beach resources continues to persist. The impact of infrastructure squeeze varies between regions, with some regions facing more severe challenges, while others experience relatively minor changes in sandy beach resources. These findings provide important theoretical support for future sandy beach resource management and protection, highlighting the balance between infrastructure development and ecological preservation. This is particularly important in the context of rapid urbanization and economic development. Achieving sustainable sandy beach resource protection will be a crucial focus for future coastal region policy planning.