Paddy Rice Phenological Mapping throughout 30-Years Satellite Images in the Honghe Hani Rice Terraces

Abstract

The Honghe Hani Rice Terraces represent the coexistence between natural and cultural systems. Despite being listed as a World Heritage Site in 2013, certain natural and anthropogenic factors have changed land use/land cover, which has led to a reduction in the size of the paddy rice area. It is difficult to accurately assess these changes due to the lack of historical maps of paddy rice croplands with fine spatial resolution. Therefore, we integrated a random forest classifier and phenological information to improve mapping accuracy and stability. We then mapped the historical distribution of land use/land cover in the Honghe Hani Rice Terraces from 1989–1991 to 2019–2021 using the Google Earth Engine. Finally, we analyzed the driving forces of land use types in the Honghe Hani Rice Terraces. We found that: (1) forests, shrubs or grasslands, and other croplands could be discriminated from paddy rice during the flooding and transplanting period, and water bodies and buildings could also be discriminated from paddy rice during the growing and harvesting period. (2) Inputting phenological feature data improved mapping accuracy and stability compared with single phenological periods. (3) In the past thirty years, 10.651%, 8.810%, and 5.711% of paddy rice were respectively converted to forests, shrubs or grasslands, and other croplands in the Honghe Hani Rice Terraces. (4) Lower agricultural profits and drought led to problems in identifying the driving mechanisms behind paddy rice distribution changes. This study demonstrates that phenological information can improve the mapping accuracy of rice terraces. It also provides evidence for the change in the size of the rice terrace area and associated driving forces in Southwest China.

1. Introduction

Paddy rice is an important caloric source that feeds more than half of the world’s population [1]. Rice fields are also vital temporary and anthropogenic wetlands that support habitats for wildlife [2]. Mapping temporal-spatial changes in paddy rice plays a crucial role in food security and ecological conservation [3]. However, applying satellite images to map paddy rice distribution presents several challenges. The similarity of spectral features between paddy rice and other crops or vegetation types throughout the year is one of the challenges for mapping paddy rice [4]. On the one hand, the spectral changes of deciduous plants, other crops, and paddy rice show a consistent positive correlation trend with precipitation and temperature changes [5]. On the other hand, different cultivation conditions, e.g., the input of nitrogen fertilizers or irrigation, also affect the spectral change of paddy rice [6]. These problems bring challenges when discriminating paddy rice from other crops or vegetation types. Therefore, determining the differences in spectral features between other crops and vegetation types at different growing phases or cultivation conditions is the key to paddy rice mapping.

The availability of satellite images with suitable time windows is another challenge for paddy rice mapping. In general, paddy rice phenology can be split into three periods: (1) the flooding and transplanting rice period (FTP); (2) the rice growing period (GP); and (3) the fallow period after harvest (FP) [7,8]. Previous studies used phenological features to map and discriminate paddy rice from other land use types by detecting its surface water and green vegetation change [8,9,10,11]. However, due to cloud cover, cloud shadow, and fog contamination, it is difficult to acquire high-quality images or suitable time window images for detecting the phenological phase of paddy rice [3,11,12]. This also limits the application of classifiers to map paddy rice based on spectral bands and vegetation indices of phenological features [3]. Therefore, a few previous studies have integrated classifiers and phenological information to map paddy rice distribution and assess mapping performance.

Historical mapping is an effective tool to monitor or assess the temporal-spatial change in land use/land cover (LULC) and ecosystem services [13,14]. Many historical LULC maps of paddy rice croplands have been produced based on satellite images. However, coarse-resolution data might lead to uncertainty in the LULC change process [14,15], making it challenging to map paddy rice and LULC at the landscape scale. Landsat data have provided finer spatial resolution and temporal coverage from 1984 to the present day [16], so Landsat has been widely used to map paddy rice at regional and national scales. We found that most current historical mappings of paddy rice based on Landsat images lacked information about other land use types, leading to problems in identifying the driving mechanisms behind paddy rice distribution changes.

Ethnic minorities, primarily the Hani people, have shaped the steep, mid-mountain slopes of the Ailao Mountains into rice terraces in a subtropical monsoon climate over the past 1300 years [17,18]. The Hani culture believes that natural environments have their divine owners, and that the Hani people serve as custodians with the right to use these natural environments, but also a duty to protect them [18,19]. Based on the traditional cultural system of “vernacular knowledge”, the Hani people maintain the balance between the local natural environment and traditional agricultural systems [20]. Thus, the Hani people have developed a sustainable landscape composed of forests, water systems, villages, and terraces [21,22]. This harmonious coexisting system of the Honghe Hani Rice Terraces (HHRT) maintains local ecosystem services and rice productivity while providing aesthetic diversity and cultural value for humans [18,20,22]. Consequently, the HHRT was designated as a Globally Important Agricultural Heritage System by UNFAO in 2010 and as a World Heritage Site by UNESCO in 2013 [20,21,23]. Mapping temporal-spatial changes of paddy rice in the HHRT is vital to local landscape planning and supporting heritage conservation.

Traditional rice terraces worldwide face various challenges from natural and anthropogenic factors [24,25]. Economic benefits, labor shortages, and food habit changes have led to rice terraces decreasing by 40% in Japan [26,27]. Soil erosion is also threatening the quality of most terraces [28]. Meteorological disasters, organic fertilizer misuse, pests, and soil erosion are the primary reasons that have led to rice terrace degeneration and loss at Ifugao in the Philippines [29]. The HHRT is also facing a series of problems: conversion of paddy rice areas to dry crops [30], farmer outmigration [31], landslide risk [17], climate change [18], etc. These problems are causing traditional culture and paddy rice area losses in the HHRT [31,32]. Most of these studies were conducted by questionnaire surveys or field monitoring, resulting in a lack of temporal-spatial quantified evidence to reveal driving forces in the rice terraces.

Here, we will investigate the traditional ecological knowledge used to acquire phenological information on paddy rice in HHRT and phenological-based paddy rice to produce historical mapping in HHRT using Landsat time-series images and the Google Earth Engine (GEE). In this study, we will address four objectives: (1) assess the separability between paddy rice and other land use types at different phenological periods; (2) assess the accuracy and stability of paddy rice mapping between three inputs of phenological information; (3) update the historical maps of LULC for HHRT from 1989–1991 to 2019–2021 using GEE and phenological information; and (4) analyze the driving forces of LULC in HHRT from 1989–1991 to 2019–2021.

2. Materials and Methods

2.1. Study Area

The HHRT (22°30′–23°30′N, 101°30′–103°30′E), distributed across four counties (Yuanyang, Honghe, Jinping, Lvchun), covers 978.071 km² in the southern Yunnan Province, China (Figure 1). The terraces are located on the south edge of the Ailao Mountains, which range in elevation from 339 m to 2865 m. Due to the lower latitudes, the study area has a typical ‘stereoscopic’ climate. The average temperature in the study area is 25 °C. The Honghe hot-dry valley in the northern part of the Ailao Mountains has the highest recorded temperature (42 °C), and the lowest recorded temperature on the highest mountain is 11.6 °C [21]. The major climate in the study area is the subtropical monsoon. Precipitation is mainly influenced by the southwest monsoon from the Indian Ocean, and the southern area of the Honghe River has the highest annual rainfall (1397.6 mm) in Yunnan Province [21]. River valleys with high evaporation and abundant rain often cause the study area to have dense fog and clouds, seriously affecting paddy rice mapping.

2.2. Phenological Features of Paddy Rice in the HHRT

2.2.1. Traditional Ecological Knowledge Investigation

We used semi-structured interviews based on the recall method to collect traditional ecological knowledge (TEK) of paddy rice phenology. We investigated TEK from eight villages in Honghe and Yuanyang Counties (Figure 1), which mainly included five Hani ethnic minority villages, two Yi ethnic minority villages, and one Dai ethnic minority village in the HHRT. Each local village head invited one or two elders to answer our TEK questions. These elders were all above 60 years old and could speak Mandarin or the local dialects. Some previously served as the heads of paddy rice production or were priests (咪谷) in their villages. Each interview was composed of four people in a relatively closed room, including an interviewer who was responsible for asking questions, an interviewer who recorded and supplemented questions, a participant who answered our questions, and a translator who was responsible for conveying the correct meaning between Mandarin or the local dialect and their ethnic minority language. The questions included the planting history of rice species in the HHRT, the traditional paddy rice farming methods in the HHRT, the phenological features and dates of each paddy rice production sector, and the conventional technology and equipment used in paddy rice production. The interviewer did not follow the list of TEK questions during the semi-structured interviews but allowed the interviewees to freely answer our questions according to their experiences and perceptions [33]. This helped to increase reciprocal interactions between interviewers and participants [34].

2.2.2. Traditional Paddy Rice Ecological Knowledge

The Hani people have plenty of TEK for paddy rice production. Local governments encourage HHRT farmers to plant traditional rice species to protect traditional cultures. Although the Hani Terrace zones have high enough temperatures and rainfall, they only plant single rice [23]. Therefore, Table 1 and Figure 2 show that the rice phenology in HHRT could be separated into three periods: (1) FTP; (2) GP; and (3) harvesting period (HP). The phenological dates used by the interviewees to answer our questions are based on the lunar calendar. In the past 20 years, the Gregorian calendar has been, on average, 36 days earlier than the lunar calendar.

As shown in Table 1, the lunar calendar month of January is the awakened month of creatures in the Hani nationality [35]. In January of the lunar calendar, the day of the year (DOY) is from 36 to 66, and the Hani people nurse rice seedlings and flood the paddies. February and March of the lunar calendar (DOY from 66 to 126) are the time for transplanting seedlings. The DOY of the FTP in the lunar calendar is from 36 to 126. In this period, the paddy rice NDVI values are lower than 0.4, and the EVI, NDSVI, and LSWI values are lower than 0.2. At the same time, the paddy rice values of NDVI, EVI, and NDSVI are lower than those of forests and shrubs or grasslands, but the paddy rice value of NDSVI is lower than that of other croplands. The paddy rice values of LSWI were lower than those of forests but higher than those of shrubs or grasslands, and other croplands (Figure 2).

The work of the Hani people from April to June of the lunar calendar is hoeing, repairing farm tools, and preparing to harvest (Table 1). This period can be defined as the GP of paddy rice and the DOY is from 126 to 216. In this period, the paddy rice value of NDVI increases from 0.4 to 0.8, EVI increases from 0.2 to 0.7, NDSVI increases from 0.2 to 0.55, and LSWI increases from 0.2 to 0.4 (Figure 2). At the same time, the vegetation indices of forests, shrubs or grasslands, and other croplands also show an increasing trend during GP.

The major rice harvesting period is from late July to August of the lunar calendar. A part of the autumn harvest will last until early September of the lunar calendar (Table 1). This period can be defined as HP and the DOY is 216 to 300. Therefore, the paddy rice values of NDVI, EVI, NDSVI, and LSWI show a decreasing trend but are higher than those of the FTP.

After the harvest period, the paddy rice in the HHRT enters the FP. The Hani people celebrate the rice harvest from the first Loong Day to the first Monkey Day in October of the lunar calendar. In addition, the paddy lands are fallowed from October to December of the lunar calendar, and the Hani people plough the paddy lands and repair ridges and farm tools (Table 1). However, the HHRT paddy flooding differs greatly from other places. The local interviewees said that most rice paddies are flooded after harvest to the transplant period. The DOY of the fallow and flooding period is from 300 to 36 of the following year. In this period, the paddy rice values of NDVI, EVI, and NDSVI are lower than those of forests, shrubs or grasslands, and other croplands. The paddy rice value of LSWI is lower than that of forests but higher than that of shrubs or grasslands, and other croplands (Figure 2).

2.3. Landsat Data Pre-Processing

Image pre-processing, phenological detection, classification, and validation were analysed by Google Earth Engine (GEE) Python API and geemap [36]. Given that Landsat satellites have several advantages, such as a long observation history, fine spatial resolution, and free access, Landsat satellite datasets are widely used to monitor LULC change and map historical datasets [37]. We selected Surface Reflectance Tier 1 (SRT1) products of the Thematic Mapper Sensor (TM) and the Operational Land Imager (OLI) to map the land use types from 1989 to 2011 and 2013 to 2021, respectively. Because the band numbers, wavelength ranges, and spatial resolution of the Enhanced Thematic Mapper Plus (ETM+) are similar to those of TM, a previous study used TM images and ETM+ images to map the historical distribution of paddy rice in Japan [4]. In this study, we combined TM images from 2010–2011 and ETM+ images from 2012 into an image collection from 2010 to 2012 to make up for missed Landsat images in 2012. A total of 1456 Landsat images were used in this study from 1989 to 2021 (Table S1), and the path/row of Landsat satellites in the HHRT are 128/44, 128/45, 129/44, 129/45, and 130/44.

Both Landsat datasets underwent atmospheric correction and orthorectification. Then, the CFmask algorithm was used to reduce the effect of clouds and cloud shadows on the images. Previous studies have proven that the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and land surface water index (LSWI) can help improve the mapping accuracy of paddy rice [8,9,38]. The normalized difference senescent vegetation index (NDSVI) can help detect vegetation senescence caused by leaf water loss [39]. Finally, we calculated four spectral indices and added these indices into spectral bands based on the GEE platform, including NDVI [40], EVI [41], NDSVI [39], and LSWI [8]:

(1)$\mathrm{NDVI}=\frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+{\rho }_{red}}$

(2)$\mathrm{EVI}=2.5\times \frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+6\times {\rho }_{red}-7.5\times {\rho }_{blue}+1}$

(3)$\mathrm{NDSVI}=\frac{{\rho }_{swir1}-{\rho }_{red}}{{\rho }_{swir1}+{\rho }_{red}}$

(4)$\mathrm{LSWI}=\frac{{\rho }_{nir}-{\rho }_{swir1}}{{\rho }_{nir}+{\rho }_{swir1}}$

where ${\rho }_{blue}$, ${\rho }_{red}$, ${\rho }_{nir}$, and ${\rho }_{swir1}$ are the reflectance of the blue, red, near-infrared (NIR), and shortwave-infrared 1 (SWIR1) bands in Landsat images, respectively.

2.4. Ground Reference Data

We collected reference data using a stratified sampling method based on Yang’s products [14]. A total of 1000 points of interest (POIs) were collected between 1999 and 2019. The POIs in 1999 and 2019 were used to map the land use types of Landsat TM images and Landsat OLI images, respectively. In this study, we defined six land use types, including paddy rice, forests, shrubs or grasslands, other croplands, water bodies, and buildings. To ensure that each land use type had enough reference data for the training and testing datasets, each one was equally allocated 50 POIs. In addition, we assigned the remaining 700 POIs based on the land use area of Yang’s dataset in 1999 and 2019. Yang’s dataset was based on GEE, and Landsat images mapped annual LULC data in China from 1990 to 2019 [14]. This dataset has six land use types in the study area: forests, shrubs, grasslands, croplands, buildings, and water bodies. Therefore, we merged the shrubs and grasslands of Yang’s dataset as a new land use type in our classification scheme, and the POIs of paddy rice were incorporated into the croplands. We randomly allocated POIs of each land use type in 1999 and 2019 by Yang’s dataset and the number of POIs. Because the average overall accuracy of Yang’s dataset is 79.30% [14], we double-checked all of the POIs on Google Earth with high spatial resolution. In the end, there were 118 and 97 POIs for paddy rice, 555 and 528 POIs for forests, 127 and 126 POIs for shrubs or grasslands, 97 and 141 for other croplands, 52 and 53 POIs for water bodies, and 51 and 55 POIs for buildings in 1999 and 2019, respectively.

Finally, we randomly selected 50 POIs of each land use type to extract the time series values of NDVI, EVI, NDSVI, and LSWI in 2019 to map the phenology features (Figure 2) and to calculate the Jeffries–Matusita distance (JMD) separability values between indices and spectral bands of the six land use types [42]. The basic idea of JMD is that if spectral or index values become more different between various land use types, the land use types are discriminated more easily [43]. The value of JMD ranges from 0 to 2, and when the value reaches 2, the land use types can be better discriminated [43]. The JMD is calculated using the equation [43,44]:

(5)$J_{ij}=2\times \left(1-e^{-B_{ij}}\right)$

in which

(6)$B_{ij}=\frac{1}{8}\times {\left(u_i-u_j\right)}^2\times \frac{2}{v_i^2+v_j^2}+\frac{1}{2}\times \ln \left(\frac{v_i^2+v_j^2}{2v_i{v}_j}\right)$

where $J_{ij}$ is the value of JMD, $B_{ij}$ is the value of Bhattacharyya distance, and $u_i$ and $v_i$ are the means and the variance of adjacent segments of class i, respectively. We calculated the JMD in R v.4.0.2.

2.5. Classification & Data Input

Random forest has higher stability, lower noise and overfitting, and higher accuracy than other machine learning algorithms [45,46]. Therefore, the random forest classifier has been widely used in remote sensing classification. The principle of the random forest classifier is based on multiple decision trees and bootstrap aggregation to train reference data [47,48]. The number of trees in the random forest is the number of repeat times in the ensemble or the number of “trees” (decision trees) in the “forest”. To reduce overfitting and maintain mapping accuracy in this study, the number of trees was set to 500 in the random forest classifier [46,49]. A random sub-selection of 70% of the POIs was used as the training dataset to map the land use, and the remaining 30% of the POIs were used as the testing dataset to assess the mapping accuracy.

Previous studies have proven that phenology-based approaches can help improve paddy rice mapping accuracy [3,4,8,9]. However, monsoonal climatic conditions in lower latitude areas that lead to persistent cloud limits the detection of paddy rice phenological features during GP [11]. To ensure we had enough Landsat images with lower cloud covers (<5%) to analyze the driving forces of paddy rice area changes, we calculated median values of Landsat spectral bands and vegetation indices during the FTP (DOY from 1 to 126) and the growing and harvesting period (GHP, DOY from 126 to 300) data input within three-year periods [4], which helped to detect the phenological features of paddy rice. We obtained 11 Landsat image collections of three-year periods: 1989–1991, 1992–1994, 1995–1997, 1998–2000, 2001–2003, 2004–2006, 2007–2009, 2010–2012, 2013–2015, 2016–2018, and 2019–2021. Among them, except for 2001–2003 (cloud cover = 27.455%) and 2013–2015 (cloud cover = 5.771%), the cloud covers of the other Landsat image collections are lower than 4%. Finally, we compared three different input datasets to assess the performance of phenological features for LULC mapping improvement. The phenological input dataset (PID) consisted of the vegetation indices (NDVI, EVI, NDSVI, and LSWI) and Landsat spectral bands (B4–7) in the time series of the FTP and GHP. Another two input datasets consisted of the vegetation indices (NDVI, EVI, NDSVI, and LSWI) and Landsat spectral bands (B4–7) in the FTP (FID) and GHP (GID), respectively. The workflow is shown in Figure 3.

2.6. Validation

To assess the accuracy between different input data, we used K-fold cross-validation of random sub-selections of 70% of POIs to map LULC in 2019–21, and the remaining 30% of POIs to assess the accuracy of thematic maps [45]. Then, we randomly selected 300 samples for accuracy validation from six periods (1989–1991, 1995–1997, 1998–2000, 2004–2006, 2010–2012, and 2013–2015). We checked these samples by visual interpretation with Landsat images and Google Earth [14]. Finally, we utilized the F1 score (FS), producer’s accuracy (PA), user’s accuracy (UA), overall accuracy (OA), and kappa coefficient (KC) to determine the accuracy of the land use maps [14].

2.7. Area Changes and Driving Force Analysis

We divided the Landsat images from 1989 to 2022 into collections of three-year periods to cover the gaps caused by clouds and cloud shadows [4,11]. In this study, we found that PID achieved the highest mapping accuracy and stability. Therefore, PID was used to map the land use distribution from 1989–1991 to 2019–2021. Because cloud cover percentage in 2001–2003 was 27.455%, we did not use the land use map in 2001–2003 to analyze driving forces. Instead, we used ArcGIS 10.5 software (version: 10.5.0.6491) and linear regression to compare the differences in land use change trends and land use transfer between HHRT and outside the HHRT from 1989–1991 to 2019–2021. Previous studies revealed that economic developments and climate change are the major driving forces impacting paddy rice area changes in China [50,51]. Therefore, economic data (GDP and proportion of primary industry) and climate data (annual temperature and annual precipitation) were used to assess the driving forces of land use types in four countries from 1989–1991 to 2019–2021 by Pearson correlation analysis in R. Among them, we collected Yunnan Statistical Yearbook economic data from 1989 to 2022 and climate data from local meteorological departments from 1989 to 2021. Finally, we randomly interviewed 37 local people in 8 villages to support driving force analysis. The average age of interviewees was 47 years old. The questionnaires included household income sources, agricultural production, employment planning, forest coverage changes, and water yield changes.

3. Results

3.1. The Separability Analysis in Two Phenological Periods

The vegetation indices and Landsat spectral bands in the FTP performed better when discriminating paddy rice from forests, shrubs or grasslands, and other croplands (Figure 4 and Table 2). The JMD values of NDVI, EVI, and NDSVI between paddy rice with forests and paddy rice with shrubs or grasslands during the FTP are more than 1. This meant that paddy rice with forests and paddy rice with shrubs or grasslands could be discriminated in the FTP by adding the vegetation indices NDVI, EVI, and NDSVI (Figure 4 and Table 2). In addition, B6 (1.059), NDSVI (0.994), B7 (0.979), and LSWI (0.842) show better separability performances than other indices or bands to discriminate paddy rice and other croplands (Figure 4 and Table 2). Furthermore, during the FTP, only B6 (1.006) can better discriminate paddy rice and water bodies, and LSWI (1.092) and B7 (1.247) can better discriminate paddy rice and buildings (Table 2).

The vegetation indices and Landsat bands in GHP performed better when discriminating paddy rice from water bodies and buildings (Figure 4 and Table 2). The vegetation indices NDVI (1.700), EVI (1.509), NDSVI (1.755) and Landsat bands B5 (1.136) and B6 (1.070) could be used to discriminate paddy rice and water bodies during GHP. At the same time, the vegetation indices of NDVI (1.607), EVI (1.148), LSWI (1.465) and Landsat bands of B4 (1.394) and B7 (1.305) performed better in discriminating paddy rice and buildings during GHP (Figure 4 and Table 2). However, only NDVI, NDSVI, and B4 during GHP could be used to discriminate paddy rice and forests. No vegetation indices or Landsat bands during GHP could be used to discriminate paddy rice with shrubs or grasslands, and paddy rice with other croplands (Table 2).

3.2. Phenological Information Improved Mapping Accuracy

By adding phenological information, the mapping accuracy and stability performed better than the other two single phenological periods. The OA and KC of PID on average are 0.871 (SD = 0.011) and 0.813 (SD = 0.015), which are higher than 0.859 (SD = 0.021) and 0.794 (SD = 0.027) of FID on average and 0.796 (SD = 0.016) and 0.695 (SD = 0.023) of GID on average, respectively. In addition, Figure 5 shows that the PA, UA, and FS of PID are higher than those of FID in paddy rice, forests, shrubs or grasslands, and water bodies. However, the UA and FS of PID are lower than those of FID in other croplands and buildings. The PA, UA, and FS of PID are higher than those of GID in each land use type. Figure 5 also shows that PID resulted in the highest stability of PA, UA, and FS compared with the other two single phenological periods. Moreover, Figure 5 also shows that the GID had better stability of PA, UA, and FS than the FID, and the PA, UA, and FS of water bodies in the GID are higher than those in the FID. Therefore, adding phenological information helps improve mapping accuracy and stability.

3.3. Paddy Rice Decreased from the 1990s to 2020s

We selected the PID using the random forest classifier to map the land use changes from 1989–1991 to 2019–2021. The average OA and KC in six periods were 88.721 ± 2.798% and 80.387 ± 4.716%, respectively (Table S2).

Paddy rice, shrubs or grasslands, and other croplands show a decreasing trend from 1989–1991 to 2019–2021 (Figure 6 and Figure 7). The loss proportions of paddy rice, shrubs or grasslands, and other croplands in HHRT are lower than those out of HHRT, but the increased proportion of forests in HHRT is higher than that outside the HHRT. Figure 7 and Table 3 show that 24.724% and 42.529% of paddy rice, 0.468% and 8.484% of shrubs or grasslands, and 51.943% and 56.879% of other croplands decreased in HHRT and outside the HHRT from 1989–1991 to 2019–2021, respectively. The forests increased to 31.572% and 17.027% in HHRT and outside the HHRT from 1989–1991 to 2019–2021, respectively. In addition, the area of paddy rice in the HHRT decreased from 13.659% in 1989–1991 to 5.566% in 2010–2012, then increased to 13.880% in 2013–2015, but decreased to 10.282% in 2019–2021. Although the HHRT was listed as a World Heritage Site in 2013, the paddy rice area decreased from 2013–2015 to 2019–2021 (Figure 7). The area of paddy rice outside the HHRT decreased from 6.347% in 1989–1991 to 2.511% in 2010–2012, then increased to 5.675% in 2013–2015, but decreased to 3.648% in 2019–2021.

Paddy rice is the major land contributor converted to other land use types from 1989–1991 to 2019–2021. Figure 6 and Table 3 show that 10.651% and 22.963%, 8.810% and 11.886%, and 5.711% and 6.640% of paddy rice areas were converted to forests, shrubs or grasslands, and other croplands in the HHRT and outside the HHRT from 1989–1991 to 2019–2021, respectively. At the same time, 34.602% of other croplands and 36.517% of shrubs or grasslands were converted to forests outside the HHRT, and 30.597% of other croplands and 38.976% of shrubs or grasslands were converted to forests in the HHRT. In addition, 23.822% and 24.615% of other croplands were converted to shrubs or grasslands in the HHRT and outside the HHRT from 1989–1991 to 2019–2021, respectively. Therefore, paddy rice and other croplands were the primary sources of the forests and shrubs or grasslands.

3.4. The Driving Factors of LULC

The Pearson correlation analysis and questionnaires revealed that climate and economic factors affected paddy rice area changes (Figure 8). Our questionnaires show that 7 household interviewees rely on agriculture for their income with an average income of 11,125 yuan per year. Working as a migrant worker is the main income source for 27 household interviewees with an average income of 50,763 yuan per year. Regarding willingness to cultivate paddy rice, 16 interviewees hoped the next generation of their family would continue to grow rice and inherit rice terrace culture. Due to the lower agricultural profits, 11 interviewees were unwilling for the next generation to engage in agricultural production. In addition, Pearson correlation analysis shows a negative correlation between GDP and paddy rice (r = −0.493) but a positive correlation between the proportion of primary industry and paddy rice (r = 0.704). The results of the questionnaires indicated that searching for higher incomes might be an important reason for local farmers abandoning agricultural lands. At the same time, we investigated the reasons for rice production changes. Fifteen interviewees thought new rice varieties improved rice production while four interviewees thought that precipitation reduction and water shortages decreased rice production. Figure 8 shows a negative correlation between annual temperature and paddy rice (r = −0.728) but a positive correlation between annual precipitation and paddy rice (r = 0.587).

4. Discussion

4.1. Phenological Features Improve Paddy Rice Mapping

Our study illustrated that adding phenological information can help improve paddy rice mapping accuracy [3,4,8,9,10,11]. Forests, water bodies, and buildings did not have obvious phenological features (Figure 2 and Figure 4), and vegetation indices or spectral bands can help to discriminate these land use types from paddy rice at different phenological periods. The vegetation indices of forests, e.g., NDVI and NDSVI, are much higher than those of paddy rice in any period (Figure 4). In particular, the JMD of forests and paddy rice in NDVI and NDSVI are higher than 1.7 during the FTP (Table 2), which can be used to discriminate forests and paddy rice at FTP. With paddy rice growing at GHP, the vegetation indices of paddy rice are higher than those of water bodies and buildings. Therefore, vegetation indices can be used to identify water bodies and paddy rice, and buildings and paddy rice at GHP [8,9].

However, deciduous vegetation, such as shrubs or grasslands, other croplands, and paddy rice, show a similar annual variation [52]. Therefore, it is difficult to identify paddy rice and other deciduous vegetation types at GHP using vegetation indices and Landsat spectral bands (Table 2 and Figure 4). Because shrubs or grasslands have plant cover during the FTP, the vegetation indices of shrubs or grasslands are higher than those of other croplands and paddy rice (Table 2 and Figure 4) [5]. Previous studies used vegetation indices or bands that could capture water information at transplanting periods to identify paddy rice and other croplands [8,9,10,53,54,55]. In our study, Figure 4 shows that SWIR (B6–7) has higher separability to discriminate other croplands and paddy rice because SWIR can detect water and water absorption [39,56]. Therefore, the calculated NDSVI and LSWI based on SWIR, can be used to discriminate other croplands and paddy rice at the FTP by detecting water differences in leaf tissues and surface moisture, respectively [8,39].

Our study found that single phenological images show a lower accuracy but had overfitting problems when mapping paddy rice or land use types. This limited feature in input datasets might be related to the lower accuracy [57]. Schulz et al. (2021) reported that incorporating time series data of phenological information can help address the overfitting problems of crop monitoring or mapping [58]. By adding phenological information of vegetation indices and Landsat spectral bands, we achieved higher mapping accuracy and stability of paddy rice and land use types than either FID or GID. Therefore, we recommend using phenological information to map paddy rice or deciduous vegetation.

4.2. The Area Changes of Paddy Rice in the Hani Terraces

The loss of paddy rice is caused by several factors [50,51], including economic development, climate, policy, and technology. Economic development has been a major factor impacting paddy rice area losses in southern China since the early 1980s [50,59]. In contrast, industrialization and urbanization have replaced paddy rice in built-up and residential areas in South China in the past 30 years [60,61]. Lower agricultural profits might be the major driver of paddy rice area losses in the HHRT [32]. The development of the secondary and tertiary industry brings a large number of employment opportunities and increased incomes. This leads to farm laborers migrating to cities for higher profits. Zhang et al. (2017b, 2020) reported that young adults prefer to find a nonfarm job in cities because a family’s annual paddy rice profit is much lower than the income from nonfarm jobs in cities [31,62]. Therefore, parts of paddy rice areas have been transferred to other croplands or left uncultivated (Table 3). The remaining labor force in rice terraces mainly consists of elderly people or uneducated women [23]. The HHRT has diverse tourism resources that could recall young entrepreneurs back to the HHRT, and support rural tourism development [20,22]. In addition, the rice-fish-duck integrated farming model not only maintains paddy rice areas and ecosystem services but can also help to improve local economic incomes [23]. Therefore, we recommend that the HHRT develop the rural tourism industry and rice-fish-duck integrated farming model to increase incomes and maintain paddy rice areas.

The climate also has an impact on the paddy rice area changes in the HHRT. Previous studies reported that climate change impacted the growing period and distribution of paddy rice [51,63]. However, continuous severe drought events are an important driving factor impacting paddy rice losses in HHRT [18]. Yang et al. (2019) reported that continuous drought in the Yunnan Province mainly occurred in 1988–1990, 2003–2007, and 2010–2015 [64]. Jiao et al. (2012) also reported that severe drought events caused terraces to dry-up in 2005 and 2010 [18]. In addition, Abbas et al. (2014) and Nichol and Abbas (2015) reported that croplands and shrublands suffered more seriously in dry conditions than forests in the Yunnan Province due to the severe drought event in 2010 [65,66]. These studies confirmed our research that a reduction in rainfall led to the paddy rice area in the HHRT in 2004–2006 (94.607 km²) and 2010–2012 (54.517 km²) being smaller than that in other periods (Figure 6 and Figure 8). The improvement of hydraulic facilities can help reduce the impact of drought on paddy rice.

Other factors could also lead to paddy rice losses in HHRT. For example, the implementation of the Grain for Green Program might lead to paddy rice being converted to forests. Landslide risk might lead to temporary or permanent loss of paddy rice areas [17].

4.3. Applications

Our method can help researchers map the crops and vegetation distribution in finer-scale Landsat data. Phenological dates, especially rice flooding and transplanting dates, are widely used to map paddy rice distribution based on Landsat images [4,10,11,12]. In this study, we acquired paddy rice phenological dates from semi-structured interviews and questionnaire interviews. This not only helped to acquire the specific phenological dates but also helped to explain the reasons for the paddy rice area change. Moreover, because of cloud and fog contamination, it is difficult to acquire enough high-quality images in different phenological periods for mapping paddy rice in the study area [4,12]. GEE provides users with high-performance computing capabilities, massive raster data, and advanced algorithms [16]. Therefore, we merged the images of the three-year time intervals based on GEE to acquire the specific phenological dates with lower cloud and fog contamination. Finally, our results only rely on Landsat spectral bands and vegetation indices, so they can be widely used to map the historical distributions of paddy rice with phenological features.

5. Conclusions

Historical maps of paddy rice are an important tool to support heritage conservation and management in the HHRT. Here, we acquired the phenological date of paddy rice by semi-structured interviews, and compared the separability of paddy rice and other land use types. We found that paddy rice could be discriminated from other land use types in different phenological periods of paddy rice. By adding phenological features, the PID achieved a higher mapping accuracy and stability than any single phenological period. Therefore, we recommend adding phenological information to map paddy rice or other crops. Moreover, we mapped the LULC with three-year periods in HHRT from 1989–1991 to 2019–2021 based on GEE and phenological information. The results show that paddy rice is the major source that was converted to other land use types, and the area losses of paddy rice in the HHRT are lower than those outside the HHRT. The questionnaires and Pearson correlation analysis revealed that lower agricultural profits and drought are the major drivers leading to paddy rice losses in the HHRT. The development of the rural tourism industry and rice-fish-duck integrated farming model could increase employment opportunities and incomes which can help to recall young entrepreneurs back to the HHRT. Thus, our thematic maps can help researchers and local governments in landscape planning and heritage conservation, and support the driving force analysis of traditional rice terrace areas.

Author Contributions

Conceptualization, J.Y. and J.X.; Methodology, J.Y., Y.Z. and D.Z.; Software, J.Y. and H.C.; Validation, J.Y., H.C. and Q.L.; Investigation, J.Y., J.X., Y.Z., Q.L. and G.Z.; Writing—Original Draft Preparation, J.Y.; Writing—Review & Editing, J.Y. and D.Z.; Supervision, J.X. and Y.Z.; Funding Acquisition, J.X. and Y.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

This research is part of Honghe Agroforestry Ecosystem Research Station of Yunnan Province, and supported by China Agricultural Green Development Association. We thanks Zhihao Qin for his valuable comments, thanks Beinn Purvis, Eyalira Jacob Okal, and Hui Cao for language editing, and thanks Haiqin Qin, Yunfen Xu, and Rui Xiao for questionnaire investigation.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The location of the HHRT in China and its elevation ranges.

Figure 2. The time series of NDVI, EVI, NDSVI, and LSWI indices for paddy rice (PR), forests (FS), shrubs or grasslands (SG), and other croplands (OC) calculated from Landsat 8 images in 2019: 50 points for each land use type.

Figure 3. The land use/land cover mapping workflow in the HHRT.

Figure 4. The annual mean of different land use types in four indices and spectral bands. FTP is the flooding and transplanting rice period, GHP is the growing and harvesting period, B4 is red band, B5 is NIR band, B6 is SWIR 1 band, B7 is SWIR 2 band, PR is paddy rice, FS is forests, SG is shrubs or grasslands, OC is other croplands, WB is water bodies, and BS is buildings.

Figure 5. Accuracy assessment with producer’s accuracy (PA), user’s accuracy (UA), and F1 score (FS) in 2019–2021 from three kinds of data input: phenological features input dataset (PID); FTP input dataset (FID); and GHP input dataset (GID). PR is paddy rice, FS is forests, SG is shrubs or grasslands, OC is other croplands, WB is water bodies, and BS is buildings.

Figure 6. The land use cover changes in the study area from 1989–1991 to 2019–2021.

Figure 7. Area changes and trends of paddy rice (PR), forests (FS), shrubs or grasslands (SG), and other croplands (OC) from 1989–1991 to 2019–2021.

Figure 8. The correlation coefficient (r) between paddy rice and its driving forces. (a–d) are the correlation coefficients between rice paddy with annual temperature, annual precipitation, GDP, and proportion of primary industry in four counties from 1989–1991 to 2019–2021, respectively.

Table 1

The Hani lunar calendar and related farming activities.

Month of Lunar Calendar	Meaning of HaniNationality	Farming Activities	Vegetation Indices Changes of Paddy Rice
January	The month of creatures awakening	Nursery rice seedlings	NDVI lower than 0.4; EVI, NDSVI, and LSWI lower than 0.2
February	The season for transplanting seedlings	Transplant seedlings and hoeing in terraces
March
April	The season of leisure	hoeing, hunting, repair farm tools	NDVI from 0.4 increased to 0.8; EVI, NDSVI, and LSWI from 0.2 increased to 0.7, 0.55, and 0.4, respectively.
May
June	The harvest preparation month	hoeing, repair farm tools, prepare for the harvest
July	The month of rice growing	Autumn harvest, fallow	NDVI decreased to lower than 0.4; EVI, NDSVI, and LSWI decreased to lower than 0.2.
August	The month of rice maturation
September	The alternate month of the new year and the last year	Harvest late rice, fallow
October	The first month of new year	Fallow, plough the paddy lands, repair ridges and farm tools	NDVI lower than 0.4; EVI, NDSVI, and LSWI lower than 0.2
November	The month of creature hibernation
December	The month of seed germination

Notes: October of the lunar calendar is the first month of the new year for the Hani nationality.

Table 2

Results of the Jeffries–Matusita distance (JMD) between paddy rice (PR) and forests (FS), shrubs or grasslands (SG), other croplands (OC), water bodies (WB), and buildings (BS).

	FTP					GHP
Indices or Bands	PR vs. FS	PR vs. SG	PR vs. OC	PR vs. WB	PR vs. BS	PR vs. FS	PR vs. SG	PR vs. OC	PR vs. WB	PR vs. BS
NDVI	1.934	1.034	0.094	0.710	0.338	1.029	0.068	0.148	1.700	1.607
EVI	1.444	1.379	0.368	0.751	0.117	0.547	0.158	0.026	1.509	1.148
NDSVI	1.707	1.445	0.994	0.991	0.338	1.188	0.516	0.030	1.755	0.850
LSWI	0.458	0.193	0.842	0.833	1.092	0.398	0.111	0.504	0.088	1.465
B4	1.118	0.062	0.339	0.237	0.462	1.038	0.057	0.650	0.793	1.394
B5	0.531	0.952	0.618	0.968	0.186	0.175	0.209	0.155	1.136	0.325
B6	0.035	0.856	1.059	1.006	0.749	0.007	0.428	0.803	1.070	0.594
B7	0.037	0.527	0.979	0.814	1.247	0.304	0.269	0.798	0.764	1.305

Table 3

The transfer matrix of land use and land cover in four counties, the HHRT, and outside the HHRT in 1989–1991 and 2019–2021 (km²).

		1989–1991
		Land Use Types	Paddy Rice	Forests	Shrubs or Grasslands	Other Croplands	Water Bodies	Buildings
2019–2021	Four counties	Paddy rice	260.243	124.580	16.231	41.091	5.538	14.199
		Forests	282.653	5629.271	1033.406	685.266	1.672	16.710
		Shrubs or Grasslands	102.447	438.720	421.858	529.110	0.764	18.378
		Other Croplands	90.305	137.487	135.016	326.121	1.158	13.950
		Water Bodies	8.045	8.346	0.926	3.036	13.283	4.629
		Buildings	16.243	16.889	11.022	23.117	1.485	17.310
	HHRT	Paddy rice	68.315	17.112	2.508	9.658	0.562	2.225
		Forests	31.323	367.917	91.377	77.892	0.166	1.392
		Shrubs or GrasslandsOther Croplands	14.261	31.550	40.941	66.429	0.032	1.540
			17.277	13.469	16.272	52.263	0.054	1.639
		Water Bodies	0.693	0.234	0.062	0.173	0.207	0.184
		Buildings	1.550	2.443	2.339	4.134	0.027	2.843
	Out ofHHRT	Paddy rice	191.935	107.459	13.722	31.425	4.976	11.962
		Forests	251.326	5261.489	942.090	607.408	1.506	15.317
		Shrubs or Grasslands	88.190	407.126	380.884	462.603	0.732	16.830
		Other Croplands	73.027	124.008	118.721	273.767	1.105	12.308
		Water Bodies	7.352	8.111	0.864	2.863	13.078	4.444
		Buildings	14.690	14.443	8.682	18.974	1.456	14.457

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs15092398/s1, Table S1: The Landsat image numbers in different phenological periods and Landsat sensors from 1989–1991 to 2019–2021; Table S2: The mapping accuracy in six periods; Supplementary information A: Paddy rice mapping code in 1998–2000.

References

1. Kuenzer, C.; Knauer, K. Remote Sensing of Rice Crop Areas—A Review. Int. J. Remote Sens.; 2013; 34, pp. 2101-2139. [DOI: https://dx.doi.org/10.1080/01431161.2012.738946]

2. Elphick, C.S. Functional Equivalency between Rice Fields and Seminatural Wetland Habitats. Conserv. Biol.; 2000; 14, pp. 181-191. [DOI: https://dx.doi.org/10.1046/j.1523-1739.2000.98314.x]

3. Dong, J.; Xiao, X. Evolution of regional to global paddy rice mapping methods: A review. ISPRS J. Photogramm. Remote Sens.; 2016; 119, pp. 214-227. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2016.05.010]

4. Carrasco, L.; Fujita, G.; Kito, K.; Miyashita, T. Historical mapping of rice fields in Japan using phenology and temporally aggregated Landsat images in Google Earth Engine. ISPRS J. Photogramm. Remote Sens.; 2022; 191, pp. 277-289. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2022.07.018]

5. Ibrahim, S.A.; Kaduk, J.; Tansey, K.; Balzter, H.; Lawal, U.M. Detecting phenological changes in plant functional types over West African savannah dominated landscape. Int. J. Remote Sens.; 2021; 42, pp. 567-594. [DOI: https://dx.doi.org/10.1080/01431161.2020.1811914]

6. Song, S.; Gong, W.; Zhu, B.; Huang, X. Wavelength selection and spectral discrimination for paddy rice, with laboratory measurements of hyperspectral leaf reflectance. ISPRS J. Photogramm. Remote Sens.; 2011; 66, pp. 672-682. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2011.05.002]

7. Toan, T.L.; Ribbes, F.; Li-Fang, W.; Floury, N.; Kung-Hau, D.; Jin Au, K.; Fujita, M.; Kurosu, T. Rice crop mapping and monitoring using ERS-1 data based on experiment and modeling results. IEEE Trans. Geosci. Remote Sens.; 1997; 35, pp. 41-56. [DOI: https://dx.doi.org/10.1109/36.551933]

8. Xiao, X.; Boles, S.; Liu, J.; Zhuang, D.; Frolking, S.; Li, C.; Salas, W.; Moore, B. Mapping paddy rice agriculture in southern China using multi-temporal MODIS images. Remote Sens. Environ.; 2005; 95, pp. 480-492. [DOI: https://dx.doi.org/10.1016/j.rse.2004.12.009]

9. Xiao, X.; Boles, S.; Frolking, S.; Li, C.; Babu, J.Y.; Salas, W.; Moore, B. Mapping paddy rice agriculture in South and Southeast Asia using multi-temporal MODIS images. Remote Sens. Environ.; 2006; 100, pp. 95-113. [DOI: https://dx.doi.org/10.1016/j.rse.2005.10.004]

10. Dong, J.; Xiao, X.; Kou, W.; Qin, Y.; Zhang, G.; Li, L.; Jin, C.; Zhou, Y.; Wang, J.; Biradar, C. et al. Tracking the dynamics of paddy rice planting area in 1986–2010 through time series Landsat images and phenology-based algorithms. Remote Sens. Environ.; 2015; 160, pp. 99-113. [DOI: https://dx.doi.org/10.1016/j.rse.2015.01.004]

11. Kontgis, C.; Schneider, A.; Ozdogan, M. Mapping rice paddy extent and intensification in the Vietnamese Mekong River Delta with dense time stacks of Landsat data. Remote Sens. Environ.; 2015; 169, pp. 255-269. [DOI: https://dx.doi.org/10.1016/j.rse.2015.08.004]

12. Jiang, M.; Xin, L.; Li, X.; Tan, M.; Wang, R. Decreasing Rice Cropping Intensity in Southern China from 1990 to 2015. Remote Sens.; 2019; 11, 35. [DOI: https://dx.doi.org/10.3390/rs11010035]

13. Johnson, D.M. Using the Landsat archive to map crop cover history across the United States. Remote Sens. Environ.; 2019; 232, 111286. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111286]

14. Yang, J.; Huang, X. The 30m annual land cover dataset and its dynamics in China from 1990 to 2019. Earth Syst. Sci. Data; 2021; 13, pp. 3907-3925. [DOI: https://dx.doi.org/10.5194/essd-13-3907-2021]

15. Sulla-Menashe, D.; Gray, J.M.; Abercrombie, S.P.; Friedl, M.A. Hierarchical mapping of annual global land cover 2001 to present: The MODIS Collection 6 Land Cover product. Remote Sens. Environ.; 2019; 222, pp. 183-194. [DOI: https://dx.doi.org/10.1016/j.rse.2018.12.013]

16. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ.; 2017; 202, pp. 18-27. [DOI: https://dx.doi.org/10.1016/j.rse.2017.06.031]

17. Gao, X.; Roder, G.; Jiao, Y.; Ding, Y.; Liu, Z.; Tarolli, P. Farmers’ landslide risk perceptions and willingness for restoration and conservation of world heritage site of Honghe Hani Rice Terraces, China. Landslides; 2020; 17, pp. 1915-1924. [DOI: https://dx.doi.org/10.1007/s10346-020-01389-4]

18. Jiao, Y.; Li, X.; Liang, L.; Takeuchi, K.; Okuro, T.; Zhang, D.; Sun, L. Indigenous ecological knowledge and natural resource management in the cultural landscape of China’s Hani Terraces. Ecol. Res.; 2012; 27, pp. 247-263. [DOI: https://dx.doi.org/10.1007/s11284-011-0895-3]

19. Yang, J.; Wang, Y.-C.; Wang, D.; Guo, L. Application of Traditional Knowledge of Hani People in Biodiversity Conservation. Sustainability; 2018; 10, 4555. [DOI: https://dx.doi.org/10.3390/su10124555]

20. Chan, J.H.; Iankova, K.; Zhang, Y.; McDonald, T.; Qi, X. The role of self-gentrification in sustainable tourism: Indigenous entrepreneurship at Honghe Hani Rice Terraces World Heritage Site, China. J. Sustain. Tour.; 2016; 24, pp. 1262-1279. [DOI: https://dx.doi.org/10.1080/09669582.2016.1189923]

21. UNESCO. Cultural Landscape of Honghe Hani Rice Terraces; UNESCO: Paris, France, 2013.

22. Gu, H.; Jiao, Y.; Liang, L. Strengthening the socio-ecological resilience of forest-dependent communities: The case of the Hani Rice Terraces in Yunnan, China. Forest Policy Econ.; 2012; 22, pp. 53-59. [DOI: https://dx.doi.org/10.1016/j.forpol.2012.04.004]

23. Yuan, Y.; Xu, G.; Shen, N.; Nie, Z.; Li, H.; Zhang, L.; Gong, Y.; He, Y.; Ma, X.; Zhang, H. et al. Valuation of Ecosystem Services for the Sustainable Development of Hani Terraces: A Rice-Fish-Duck Integrated Farming Model. Int. J. Environ. Res. Public Health; 2022; 19, 8549. [DOI: https://dx.doi.org/10.3390/ijerph19148549] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35886399]

24. Zhang, Y.; He, L.; Li, X.; Zhang, C.; Qian, C.; Li, J.; Zhang, A. Why are the Longji Terraces in Southwest China maintained well? A conservation mechanism for agricultural landscapes based on agricultural multi-functions developed by multi-stakeholders. Land Use Pol.; 2019; 85, pp. 42-51. [DOI: https://dx.doi.org/10.1016/j.landusepol.2019.03.046]

25. Yu, M.; Li, Y.; Luo, G.; Yu, L.; Chen, M. Agroecosystem composition and landscape ecological risk evolution of rice terraces in the southern mountains, China. Ecol. Indic.; 2022; 145, 109625. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.109625]

26. Sarker, A. Economics of underproduction: A polycentric approach for a depopulated commons in Japan. Ecol. Econ.; 2020; 171, 106597. [DOI: https://dx.doi.org/10.1016/j.ecolecon.2020.106597]

27. Yasushi, M.; Masaya, S.; Eisei, M.; Kumiko, T. When do rice terraces become rice terraces?. Paddy Water Environ.; 2019; 17, pp. 323-330. [DOI: https://dx.doi.org/10.1007/s10333-019-00727-0]

28. Arnáez, J.; Lana-Renault, N.; Lasanta, T.; Ruiz-Flaño, P.; Castroviejo, J. Effects of farming terraces on hydrological and geomorphological processes. A review. Catena; 2015; 128, pp. 122-134. [DOI: https://dx.doi.org/10.1016/j.catena.2015.01.021]

29. Joshi, R.C.; Matchoc, O.R.O.; Bahatan, R.G.; Pena, F.A.D. Farmers’ knowledge, attitudes and practices of rice crop and pest management at Ifugao Rice Terraces, Philippines. Int. J. Pest Manag.; 2000; 46, pp. 43-48. [DOI: https://dx.doi.org/10.1080/096708700227561]

30. Yang, L.; Liu, M.; Lun, F.; Yuan, Z.; Zhang, Y.; Min, Q. An Analysis on Crops Choice and Its Driving Factors in Agricultural Heritage Systems—A Case of Honghe Hani Rice Terraces System. Sustainability; 2017; 9, 1162. [DOI: https://dx.doi.org/10.3390/su9071162]

31. Zhang, Y.; Min, Q.; Zhang, C.; He, L.; Zhang, S.; Yang, L.; Tian, M.; Xiong, Y. Traditional culture as an important power for maintaining agricultural landscapes in cultural heritage sites: A case study of the Hani terraces. J. Cult. Herit.; 2017; 25, pp. 170-179. [DOI: https://dx.doi.org/10.1016/j.culher.2016.12.002]

32. Zhang, Y.; Li, X.; Min, Q. How to balance the relationship between conservation of Important Agricultural Heritage Systems (IAHS) and socio-economic development? A theoretical framework of sustainable industrial integration development. J. Clean. Prod.; 2018; 204, pp. 553-563. [DOI: https://dx.doi.org/10.1016/j.jclepro.2018.09.035]

33. Kallio, H.; Pietilä, A.-M.; Johnson, M.; Kangasniemi, M. Systematic methodological review: Developing a framework for a qualitative semi-structured interview guide. J. Adv. Nurs.; 2016; 72, pp. 2954-2965. [DOI: https://dx.doi.org/10.1111/jan.13031] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27221824]

34. Zhang, J.-Q.; Mammides, C.; Corlett, R.T. Reasons for the Survival of Tropical Forest Fragments in Xishuangbanna, Southwest China. Forests; 2020; 11, 159. [DOI: https://dx.doi.org/10.3390/f11020159]

35. Yang, J. Adaptation Mechanism on Climate Change of Traditional Agriculture Management in Hani Rice Terraces. Ph.D. Thesis; Beijing Forestry University: Beijing, China, 2015.

36. Wu, Q. geemap: A Python package for interactive mapping with Google Earth Engine. J. Open Source Softw.; 2020; 5, 2305. [DOI: https://dx.doi.org/10.21105/joss.02305]

37. Wulder, M.A.; White, J.C.; Loveland, T.R.; Woodcock, C.E.; Belward, A.S.; Cohen, W.B.; Fosnight, E.A.; Shaw, J.; Masek, J.G.; Roy, D.P. The global Landsat archive: Status, consolidation, and direction. Remote Sens. Environ.; 2016; 185, pp. 271-283. [DOI: https://dx.doi.org/10.1016/j.rse.2015.11.032]

38. Zhang, G.; Xiao, X.; Dong, J.; Xin, F.; Zhang, Y.; Qin, Y.; Doughty, R.B.; Moore, B. Fingerprint of rice paddies in spatial–temporal dynamics of atmospheric methane concentration in monsoon Asia. Nat. Commun.; 2020; 11, 554. [DOI: https://dx.doi.org/10.1038/s41467-019-14155-5]

39. Qi, J.; Marsett, R.; Heilman, P.; Bieden-bender, S.; Moran, S.; Goodrich, D.; Weltz, M. RANGES improves satellite-based information and land cover assessments in southwest United States. Eos Trans. Am. Geophys. Union; 2002; 83, pp. 601-606. [DOI: https://dx.doi.org/10.1029/2002EO000411]

40. Defries, R.S.; Townshend, J.R.G. NDVI-derived land cover classifications at a global scale. Int. J. Remote Sens.; 1994; 15, pp. 3567-3586. [DOI: https://dx.doi.org/10.1080/01431169408954345]

41. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ.; 2002; 83, pp. 195-213. [DOI: https://dx.doi.org/10.1016/S0034-4257(02)00096-2]

42. Thomas, I.L.; Ching, N.P.; Benning, V.M.; D’Aguanno, J.A. Review Article A review of multi-channel indices of class separability. Int. J. Remote Sens.; 1987; 8, pp. 331-350. [DOI: https://dx.doi.org/10.1080/01431168708948645]

43. Jung, R.; Ehlers, M. Comparison of two feature selection methods for the separability analysis of intertidal sediments with spectrometric datasets in the German Wadden Sea. Int. J. Appl. Earth Obs. Geoinf.; 2016; 52, pp. 175-191. [DOI: https://dx.doi.org/10.1016/j.jag.2016.06.009]

44. Wang, Y.; Qi, Q.; Liu, Y. Unsupervised Segmentation Evaluation Using Area-Weighted Variance and Jeffries-Matusita Distance for Remote Sensing Images. Remote Sens.; 2018; 10, 1193. [DOI: https://dx.doi.org/10.3390/rs10081193]

45. Yang, J.; Xu, J.; Zhai, D.-L. Integrating Phenological and Geographical Information with Artificial Intelligence Algorithm to Map Rubber Plantations in Xishuangbanna. Remote Sens.; 2021; 13, 2793. [DOI: https://dx.doi.org/10.3390/rs13142793]

46. Belgiu, M.; Drăguţ, L. Random forest in remote sensing: A review of applications and future directions. ISPRS J. Photogramm. Remote Sens.; 2016; 114, pp. 24-31. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2016.01.011]

47. Breiman, L. Random Forests. Mach. Learn.; 2001; 45, pp. 5-32. [DOI: https://dx.doi.org/10.1023/A:1010933404324]

48. Blickensdörfer, L.; Schwieder, M.; Pflugmacher, D.; Nendel, C.; Erasmi, S.; Hostert, P. Mapping of crop types and crop sequences with combined time series of Sentinel-1, Sentinel-2 and Landsat 8 data for Germany. Remote Sens. Environ.; 2022; 269, 112831. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112831]

49. Rodriguez-Galiano, V.F.; Ghimire, B.; Rogan, J.; Chica-Olmo, M.; Rigol-Sanchez, J.P. An assessment of the effectiveness of a random forest classifier for land-cover classification. ISPRS J. Photogramm. Remote Sens.; 2012; 67, pp. 93-104. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2011.11.002]

50. Zhang, G.; Xiao, X.; Biradar, C.M.; Dong, J.; Qin, Y.; Menarguez, M.A.; Zhou, Y.; Zhang, Y.; Jin, C.; Wang, J. et al. Spatiotemporal patterns of paddy rice croplands in China and India from 2000 to 2015. Sci. Total Environ.; 2017; 579, pp. 82-92. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.10.223]

51. Dong, J.; Xiao, X.; Zhang, G.; Menarguez, M.A.; Choi, C.Y.; Qin, Y.; Luo, P.; Zhang, Y.; Moore, B. Northward expansion of paddy rice in northeastern Asia during 2000–2014. Geophys. Res. Lett.; 2016; 43, pp. 3754-3761. [DOI: https://dx.doi.org/10.1002/2016GL068191]

52. Bridhikitti, A.; Overcamp, T.J. Estimation of Southeast Asian rice paddy areas with different ecosystems from moderate-resolution satellite imagery. Agric. Ecosyst. Environ.; 2012; 146, pp. 113-120. [DOI: https://dx.doi.org/10.1016/j.agee.2011.10.016]

53. Pan, X.-Z.; Uchida, S.; Liang, Y.; Hirano, A.; Sun, B. Discriminating different landuse types by using multitemporal NDXI in a rice planting area. Int. J. Remote Sens.; 2010; 31, pp. 585-596. [DOI: https://dx.doi.org/10.1080/01431160902894442]

54. Kordelas, G.A.; Manakos, I.; Aragonés, D.; Díaz-Delgado, R.; Bustamante, J. Fast and Automatic Data-Driven Thresholding for Inundation Mapping with Sentinel-2 Data. Remote Sens.; 2018; 10, 910. [DOI: https://dx.doi.org/10.3390/rs10060910]

55. Verma, A.; Kumar, A.; Lal, K. Kharif crop characterization using combination of SAR and MSI Optical Sentinel Satellite datasets. J. Earth Syst. Sci.; 2019; 128, 230. [DOI: https://dx.doi.org/10.1007/s12040-019-1260-0]

56. Du, Y.; Zhang, Y.; Ling, F.; Wang, Q.; Li, W.; Li, X. Water Bodies’ Mapping from Sentinel-2 Imagery with Modified Normalized Difference Water Index at 10-m Spatial Resolution Produced by Sharpening the SWIR Band. Remote Sens.; 2016; 8, 354. [DOI: https://dx.doi.org/10.3390/rs8040354]

57. Barr, M. A Novel Technique for Segmentation of High Resolution Remote Sensing Images Based on Neural Networks. Neural Process. Lett.; 2020; 52, pp. 679-692. [DOI: https://dx.doi.org/10.1007/s11063-020-10280-1]

58. Schulz, C.; Holtgrave, A.-K.; Kleinschmit, B. Large-scale winter catch crop monitoring with Sentinel-2 time series and machine learning–An alternative to on-site controls?. Comput. Electron. Agric.; 2021; 186, 106173. [DOI: https://dx.doi.org/10.1016/j.compag.2021.106173]

59. Xin, F.; Xiao, X.; Dong, J.; Zhang, G.; Zhang, Y.; Wu, X.; Li, X.; Zou, Z.; Ma, J.; Du, G. et al. Large increases of paddy rice area, gross primary production, and grain production in Northeast China during 2000–2017. Sci. Total Environ.; 2020; 711, 135183. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.135183]

60. Liu, J.; Zhang, Z.; Xu, X.; Kuang, W.; Zhou, W.; Zhang, S.; Li, R.; Yan, C.; Yu, D.; Wu, S. et al. Spatial patterns and driving forces of land use change in China during the early 21st century. J. Geogr. Sci.; 2010; 20, pp. 483-494. [DOI: https://dx.doi.org/10.1007/s11442-010-0483-4]

61. Ghose, B. Food security and food self-sufficiency in China: From past to 2050. Food Energy Secur.; 2014; 3, pp. 86-95. [DOI: https://dx.doi.org/10.1002/fes3.48]

62. Zhang, Y.; He, L.; Liu, L.; Min, Q. Labor productivity of small-scale agriculture and its influence on agricultural landscape conservation in mountainous areas in China: A case study of rice farming in Hani terraced region. Environ. Sci. Pollut. Res.; 2020; 27, pp. 39795-39806. [DOI: https://dx.doi.org/10.1007/s11356-020-08815-z]

63. Wang, X.; Ciais, P.; Li, L.; Ruget, F.; Vuichard, N.; Viovy, N.; Zhou, F.; Chang, J.; Wu, X.; Zhao, H. et al. Management outweighs climate change on affecting length of rice growing period for early rice and single rice in China during 1991–2012. Agric. For. Meteorol.; 2017; 233, pp. 1-11. [DOI: https://dx.doi.org/10.1016/j.agrformet.2016.10.016]

64. Yang, C.; Tuo, Y.; Ma, J.; Zhang, D. Spatial and Temporal Evolution Characteristics of Drought in Yunnan Province from 1969 to 2018 Based on SPI/SPEI. Water Air Soil Pollut.; 2019; 230, 269. [DOI: https://dx.doi.org/10.1007/s11270-019-4287-6]

65. Abbas, S.; Nichol, J.E.; Qamer, F.M.; Xu, J. Characterization of Drought Development through Remote Sensing: A Case Study in Central Yunnan, China. Remote Sens.; 2014; 6, pp. 4998-5018. [DOI: https://dx.doi.org/10.3390/rs6064998]

66. Nichol, J.E.; Abbas, S. Integration of remote sensing datasets for local scale assessment and prediction of drought. Sci. Total Environ.; 2015; 505, pp. 503-507. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2014.09.099]