Decadal changes in the rice-cropping system in the Ayeyarwady Delta using a large archive of satellite imagery from 1981 to 2020

In the last 40 years, the rice-cropping system has considerably changed in the Ayeyarwady Delta. The large archive of satellite imagery provides a history of how land and water resource managements have changed in the face of growing populations, resource demand, and climate change. This study aimed to assess the decadal changes in the rice-cropping system in the Ayeyarwady Delta by using the large archive of satellite imagery for the last 40 years (1981 − 2020). The long-term NDVI dataset provided various information on rice cultivation. Signal processing techniques were used to detect on the historical changes in the rice-cropping system, and the impact of climate change was assessed by using trend analysis. Until the 1980s, single-cropping of summer rice was dominant in the Delta. To enhance the grain yield of rice, the irrigation facilities were introduced in 1992 under an initiative of the Myanmar government. As a result, the annual cropping intensities increased from 1.087 ± 0.390 in the 1980s to 1.422 ± 0.499 in the 2010s. The information on historical change in the rice-cropping system would be useful to consider the practical and cost-effective utilization of remaining land and water resources. Moreover, the trend analysis of NDVI time-series showed negative trends in coastal areas. This indicates that the rice production in coastal areas has been constrained by the saline intrusion. The salt-affected areas are expected to expand under future climate change scenarios. Government support is highly required for sustainable rice production in the Delta.


Introduction
Agriculture is the largest economic sector in Myanmar, and thus the sustainable growth strategy for rice production is the main economic policy (Maung 1998). Successive Myanmar governments have attempted to expand the arable land area to maximize the rice production. Since 1988, the Myanmar government has encouraged farmers to make full use of unused and fallow lands (Kurosaki 2008). However, the rice production per hectare is low, mainly because of the traditional practice of single-cropping. Infrastructure improvements are indispensable to enhance the land-use efficiency for rice production. Irrigation facilities make it possible to plant rice in seasons when rainwater is limited (Maung and Satoh 2004). Improving productivity through more efficient utilization of scarce land resources is the best and highly effective way. The Myanmar government has enthusiastically promoted the provision of irrigation facilities into agricultural land as a way of boosting rice production since 1992 (Matsuda 2009). Overall, with the combination of supportive policies and technologies, rice production has increased over several decades (Garcia et al. 2000). However, the implementation of irrigation facilities did not necessarily guarantee higher returns to farmers because they incurred more cost for water resources (Oo et al. 2017). Competition for water resources is common during the dry season, especially where there is a limited regulation and an absence of cooperative water management. Irrigation scheduling, which is defined as the decision of when and how much water need to apply to a field, is crucial 1 3 in intensive paddy fields (Pardossi et al. 2009). Therefore, it is necessary to understand the distribution pattern of the rice-cropping system for apposite water resource management. Moreover, several studies suggest that increasing temperature, rising sea level, and changes in patterns of rainfall and its distribution under climate change could lead to a modification in land and water resources for rice production (Wassmann et al. 2009;Zabel et al. 2014;Chun et al. 2016;Shrestha and Htut 2016). There are also serious implications for food security if no effort is made to adapt rice production to climate change.
Satellite remote sensing has enabled long-term observation of land-use changes and vegetation dynamics because of the consistent and repeatable data (Herrmann et al. 2005;Hill et al. 2008;Bégué et al. 2011;Schucknecht et al. 2013;Jiang et al. 2017). The Advanced Very High-Resolution Radiometer (AVHRR) on the National Oceanic and Atmospheric Administration (NOAA) polar-orbiting satellites provides the longest record with global coverage since 1981 . AVHRR and alternative sensors have continued to provide records for 40 years. In general, the normalized difference vegetation index (NDVI) is used as the most elementary indicator of vegetation activity (Ruimy et al. 1994). NDVI is strongly associated with the biophysical properties of vegetation, such as biomass, leaf area index, or crop yield (Box et al. 1989;Quarmby et al. 1993;Fensholt et al. 2004). Therefore, the time-series of NDVI approximates the stage of the growth cycle, including germination, leaf emergence, and start of senescence. Many studies have focused on extraction of phenology types, phenology change detection, trends in phenology, and the relationships between phenology and climate change (Fuller 1998;Chen et al. 2001;Hermance et al. 2007;Julien and Sobrino 2009).
In the time-series analysis, signal processing is often used to distinguish between signal and noise. Fourier analysis is widely known as a typical method because of its prowess and simplicity (Sakamoto et al. 2005;Canisius et al. 2007;Mingwei et al. 2008). Fourier analysis presupposes the stationarity and periodicity of data. However, the frequency components are not always consistent. In Myanmar, the ricecropping system changed from single to double-cropping, and the timing varied with the locations considerably. It is impossible to understand accurately the frequency characteristics from such data with strong nonstationarity. Wavelet analysis is a method of obtaining the time change of the intensity of each frequency component by enlarging, reducing, and translating a basis function called a mother wavelet (Sakamoto et al. 2005(Sakamoto et al. , 2006Martínez and Gilabert 2009). Wavelet analysis is effective for detecting the frequency characteristics of data with strong nonstationarity. However, wavelet analysis has the inconvenience that it is difficult to select the optimum mother wavelet (Tse et al. 2004). Due to the deficiencies of wavelet analysis, a new type of time-frequency analysis called Hilbert-Huang Transform (HHT) was proposed for analyzing nonstationary signals (Huang et al. 1996). In HHT, data is first decomposed into a simple Intrinsic Mode Function (IMF) by a method called Empirical Mode Decomposition (EMD), and then Hilbert Spectrum analysis is performed on each IMF to obtain the instantaneous frequency. HHT has higher time and frequency resolution than conventional analysis methods, and its usefulness has been shown through analysis results for many natural phenomena at various spatial scales from local to global. Some of the applications of HHT in the geophysical studies have been reviewed by Huang and Wu (2008). However, there have been few studies still using remote sensing data (Chen et al. 2012;Son et al. 2017).
The objectives of this study are (1) to understand the decadal changes in the rice-cropping system in the Ayeyarwady Delta, Myanmar, by using signal processing techniques, (2) to assess the impact of climate change on vegetation dynamics by using trend analysis with the long-term archive of NDVI time-series. To ensure the efficient utilization of the limited land and water resources, monitoring the rice-cropping system is a major research priority in Myanmar. From a historical point of view, reviewing long-term changes in the rice-cropping system is important to predict and prepare for future rice production and water resource requirements. The large archive of satellite-derived dataset provides a history of how land and water management has changed over the past 40 years in the face of the growing population, resource demand, and climate change.

Study area
In Myanmar, there are three principal rice zones: the delta zone, the dry zone, and the hill zone (Matsuda 2009). In this study, we focused on the delta zone as the paddy field occupies most of the agricultural land. The Ayeyarwady Delta (hereafter the Delta) is located at the center of Myanmar, between 94° and 97° E longitude, and 15.6° and 18.6° N latitude. Figure 1 shows the land cover map from the United Nations Environment Programme (UNEP) 2000 data. The Delta is a flat alluvial plain with a network of tributaries of the Ayeyarwady River. The Ayeyarwady River deposits enormous quantities of sediments. Therefore, the Delta is characterized by high soil quality, in comparison with other rice zones of Myanmar. The Delta belongs to a tropical monsoon climate with two seasons: rainy season and dry season. The annual precipitation ranges from 2500 to 4000 mm. More than 1 3 90% of precipitation falls during the rainy season. During the dry season, rice and other crops cannot grow without irrigation facility.
According to the Ministry of National Planning and Economic Development of Myanmar, rice production is 88.8% of the total crop production in the Delta (Zaw et al. 2011, Supplementary Figure S1). Although the contribution of rice production is high in the Delta, pulses and beans are also cultivated in agricultural lands. The main purpose of this study is to investigate the utilization situation of irrigation facilities in the Delta. Therefore, we divided into two crop types from the perspective of the cropping system. Crops grown in the rainy season (from July to December) and the dry season (from January to June) are defined as the monsoon crop and summer crop, respectively.
The low elevations of the Delta are exceptionally vulnerable to saline intrusion, storm surge, sea-level rising, and other impacts of climate change. To investigate the effects of salt damage on crop management, a detailed analysis was conducted in Kangyidaunt and Labutta Townships (Fig. 1b). Labutta Township faces the coast and has been already affected by the saline intrusion, while Kangyidaunt Township is located approximately 80 km inland from the coast and has not yet been affected by saline intrusion.

Large archive of satellite imagery
We used the NDVI dataset in the Global Vegetation Health Products provided by NOAA (https:// www. star. nesdis. noaa. gov) for monitoring the changes of the cropping system in the Delta. The dataset is blended by large archives of satellite imagery from AVHRR onboard afternoon polar-orbiting satellites (the 9,11,14,16,18 and 19) from 1981 to 2012 and the Visible/Infrared Imager and Radiometer Suite (VIIRS) from the Suomi National Polar-orbiting Partnership (Suomi-NPP) satellite from 2013 to the present. NDVI was calculated from the red and near-infrared (NIR) spectral bands, according to the formula in Eq. (1).
The bandwidths of the red (640-680 nm) and NIR (850-880 nm) bands for Suomi-NPP are much narrower than the red (580-680 nm) and NIR (725-1000 nm) bands for AVHRR. Since NDVI has high-frequency noises related to variable transparency of the atmosphere, bidirectional reflectance, orbital drift, and sensor degradation, which makes it difficult to use. The noises are removed by applying statistical techniques (Kogan 1997). The NDVI dataset is no noise (smoothed) 7-day composite with a resolution of 4 × 4 km. Although other datasets with higher spatial resolution, such as Landsat or MODIS, are also extremely useful for mapping the agricultural land, we used the long-term high-frequency dataset to detect historical changes in the cropping system over last 40 years.

Peak analysis
Peak analysis was used as a classical method to evaluate the cropping system. Although we used the dataset of the smoothed 7-day NDVI, there were still some noises on the curves. The estimation of special points, such as maxima and inflection points, requires smooth phenology curves. Therefore, a Savitzky-Golay filter was applied to fit the time-series of NDVI (Savitzky and Golay 1964). In Myanmar, the cropping intensity is single or double, and thus two parameters for the Savitzky-Golay filter (i.e., flame length and order) were adjusted so that there were only one or two peaks in each year. As a result, the flame length was set to 9, and order was set to 1. After smoothing, the maxima was easily detected, as the first derivative equals 0 and changes from positive to negative at the maxima point. The number of maxima in phenology curves was then counted to distinguish between single-cropping and double-cropping.

Time-frequency analysis (HHT)
A novel time-frequency analysis method based on the HHT was applied to characterize of vibration signals (Huang et al. 1996). The HHT consists of two parts: EMD and the Hilbert spectral analysis. EMD is the core method of the HHT. According to EMD, the original time-series data can be decomposed into a series of IMF components. Once we have obtained the IMFs according to EMD, then the Hilbert spectral analysis, which is based on the IMFs and the Hilbert transform, can be obtained. Firstly, y(t) is calculated by the Hilbert transform, as shown in Eq. (2).
where P indicates the Cauchy principal value. With the Hilbert transform y(t) of the function x(t), we obtain the analytic function, where a is the instantaneous amplitude, and θ is the instantaneous phase function. The instantaneous frequency is simply With both amplitude and frequency being a function of time, we can express the amplitude (or energy, the square of amplitude) in terms of function of time and frequency, H(ω, t). The marginal spectrum can then be defined as where [0, T] is the temporal domain within which the data are defined.

Trend analysis (Mann-Kendall test and Sen's slope)
A trend analysis was performed in this study by using the Mann-Kendall test (Mann 1945;Kendall 1975). The Mann-Kendall is a nonparametric method that is robust against nonnormality of a dataset. The Mann-Kendall test statistic (S) is calculated as: where n is the length of the sample, x i and x j are the data values in years i and j (j > i), respectively, and sgn(x jx i ) is the sign function as: The variance is calculated as: where q is the number of tied groups and t i is the number of observations in the ith group. Then, the standardized test statistic value (Z) is calculated by utilizing the following equations.
The presence of a statistically significant trend is evaluated using the Z value. Positive values of Z indicate increasing trends while negative Z values show decreasing trends. In this study, this statistic is used to identify pixels, for which the p-value was less than 0.05. Besides, the magnitude of trends in the time-series data was evaluated by a simple nonparametric procedure developed by Sen (1968). The trend is calculated by using the following equation: where β is Sen's slope. Sen's slope is insensitive to outliers, and much more accurate than the slope of simple linear regression for skewed and heteroskedastic data. β > 0 indicates increasing trends in a time series. Otherwise, the data series presents decreasing trends during the period. Figure 2 shows examples of typal phenological cycle of four land cover types, i.e., agriculture, scrubland, deciduous forest and evergreen forest, in Fig. 1. The agriculture found two types of phenological cycle, expressed as single and double-cropping. There was no triple-cropping in Myanmar. The time-series of NDVI in agriculture was collected from the same one-pixel location within Kangyidaunt Township in the 1980s and 2010s. In the 1980s, one peak of NDVI, ranging between 0.2 and 0.6, appeared during the monsoon season, and NDVI values were close to zero during the summer season. In the 2010s, the phenology found two seasonal cycles in one calendar year. The peaks appeared during the monsoon season and the summer season, and valleys appeared between the peaks. The valley from the summer

Phenological cycle using NDVI
season to the monsoon season was deeper with NDVI values of less than 0.2, while the valley from the monsoon season to the summer season was relatively shallow with NDVI values of around 0.3. The characteristic difference in the phenological cycle between single-cropping and double-cropping was noticeable in maximum NDVI during the summer season. By comparing the phenological cycles of the 1980s and the 2010s in Kangyidaunt Township, the change in the cropping system could be determined. Phenological cycles of scrubland and forests also showed a distinctive seasonal patterns on NDVI variation. The shapes of the phenological cycle were basically similar to agriculture with one seasonal peak, but the peak and valley values and its timing are different among land cover types. The peak NDVI of forests was the highest, ranging from 0.5 to 0.7. NDVI values tended to saturate at 0.7 in the forests with high LAI. Wide ranges of NDVI were shown at the valley due to the difference in evergreen and deciduous forests. Evergreen forest had higher NDVI values at the valley. NDVI values of deciduous forest were close to zero at the valley. NDVI values of scrubland were between agriculture and forests. In this study, we focused only on change in the cropping system in agriculture, and other land cover types were ignored. Figure 3 shows the time-series of NDVI from 1981 to 2020 in Kangyidaunt and Labutta Townships. NDVI was spatially aggregated by using polygons of each Township. In Kangyidaunt Township, the phenological cycle was once a year until 1993, but shifted to twice a year since 1994 (Fig. 3a). The turning point from once to twice almost coincided with the year when irrigation facilities were comprehensively introduced in 1992 by the Myanmar government.  On the other hand, the phenological cycle was once a year in Labutta Township from 1981 to 2020, although sometimes two cycles appeared (Fig. 3b). Closed and open circles in Fig. 3 show peak NDVI detected by peak analysis during the monsoon and summer seasons, respectively. Figure 4 shows five IMFs obtained from initial data of Kangyidaunt and Labutta Townships in Fig. 3. The initial data were decomposed into five IMFs from high-frequency IMF (C 1 ) to low frequency IMF (C 5 ) by EMD. The sum of all IMFs matches the initial data. The first IMF (C 1 ) was more symmetric with simple vibration mode than the initial data. Therefore, it became easier to visually distinguish the difference in phenological cycle between the single-cropping and double-cropping. After decomposing the initial data into multiple IMFs with EMD, the HHT analysis was performed on each IMF. Figure 5 shows the instantaneous frequency spectrum of the first IMF (C 1 ) in Kangyidaunt and Labutta Townships. Most frequencies were less than 0.2 in Labutta Township with phenological cycle of singlecropping; however, the instantaneous frequencies in Kangyidaunt Township showed high values of approximately 0.4 or more since 1994 because of shift from single-cropping to double-cropping. It might be possible to confirm changes in cropping system with the frequency in the Hilbert Spectrum. Figure 6 shows the histogram of the accumulated annual frequencies. There were two distinct peaks of the accumulated annual frequency at 6.3 and 12.6 due to the differences in phenological cycle. Histogram with a valley between the two peaks is useful for a threshold value. In this study, we assumed to be the double-cropping when the accumulated annual frequencies were more than 10. Figure 7 shows the cropping intensity in the Delta by using peak analysis (counting the number of maxima, Fig. 7a-d) and time-frequency analysis (distinguishing by difference in frequency, Fig. 7e-h). In Myanmar, the cropping intensity is 1 or 2. However, as we calculated the decadal average for each pixel, the values were in the range of 1-2. For example, a cropping intensity of 1.5 means that single-cropping and double-cropping were conducted for 5 years each. The distribution patterns of cropping intensity obtained from both analyzes were almost same. In the 1980s, single-cropping was conducted in most agricultural lands of the Delta. Since the 1990s, double-cropping rapidly increased along the main rivers. However, double-cropping was not conducted on the southwest coast of the Delta. The overall cropping intensities of the Delta in the 1980s, 1990s, 2000s, and 2010s were 1.044 ± 0.327, 1.244 ± 0.365, 1.387 ± 0.448 and 1.386 ± 0.483 (mean ± SD) by peak  analysis, and 1.087 ± 0.390, 1.225 ± 0.455, 1.367 ± 0.498 and 1.422 ± 0.499 by time-frequency analysis, respectively. Figure 8 shows the time-series of cropping intensity derived from satellite images and the sown area of rice derived from agricultural statistics by Department of Agriculture during 1988 − 2019. These showed similar fluctuation patterns. The cropping intensities and the sown area were low until around 1992, and then increased. The sown area of rice increased 1.65 times from 1.28 to 2.09 × 10 4 km 2 in the Delta from 1988 to 2019. Figure 9 shows the distribution maps of maximum NDVI (NDVI max ) during the monsoon and the summer seasons in the 1980s and 2010s. The distribution patterns of NDVI max in the 1980s were obviously different between the seasons (Fig. 9a, b). NDVI max in the monsoon season was higher for most areas in the Delta (0.397 ± 0.049), while NDVI max in the summer season was substantially lower (0.243 ± 0.073). In the Yangon region in the southeast of the Delta, NDVI max was almost zero. In the 2010s, the difference in NDVI max between the seasons became smaller (Fig. 9c, d). NDVI max during the monsoon and summer seasons were 0.376 ± 0.071 and 0.341 ± 0.091, respectively. The trend of NDVI max for 40 years was analyzed by using the Mann-Kendall test together with Sen's slope estimator (Fig. 9e, f). There were statistically significant positive trends of NDVI max during the summer season at large part of area (p < 0.05), except for coastal areas (Fig. 9f), although not significant positive trends during the monsoon season (Fig. 9e). The increase in NDVI max during the summer season seems to be the result of the changes in cropping system. The Sen's slope during the summer season tended to be positively higher in some areas with double-cropping (Fig. 7). In contrast, the Sen's slope during the monsoon season tended to be negative in the double-cropping areas, because the growing period for monsoon crop became shorter by the initiation of double-cropping.

Trends of NDVI max for 40 years
During both the monsoon and summer seasons, the Sen's slope showed negative trends along the coastline in the southwest Delta. It seems that the cause is the deterioration of saline intrusion due to the sea-level rise. The salt-affected areas analyzed by Fee et al. (2017) corresponded to the areas where the Sen's slope was negative (Fig. 9e, f). The effect of the saline intrusion was shown in Labutta Township but not in Kangyidaunt Township. Therefore, we investigated in detail whether there were differences in NDVI between the two Townships (Fig. 3). In Kangyidaunt Township, NDVI max during the monsoon season was constantly high, ranging from 0.3 and 0.4, and did not show a statistically significant trend. On the other hand, the time-series of NDVI max during the monsoon season showed a negative trend in Labutta Township, and NDVI max showed relatively high values of more than 0.3 in the 1980s but values declined to less than 0.3 from the late 2000s onwards. In addition, in Labutta Township, double-cropping was conducted only in some parts of the north where there was no influence of saline intrusion (Fig. 9). When NDVI data were calculated for the whole Township, the time-series of NDVI in Labutta Township showed one peak during the monsoon season over the years.

Change in cropping system
In this study, two analytical methods of signal processing, i.e., peak analysis and time-frequency analysis, were examined to better understand the decadal changes in cropping system in the Delta by using archives of the 7-days NDVI dataset for 40 years. In the 1980s, the phenological cycle in agriculture had one peak during the monsoon season ( Figs. 2 and 3), because monsoon crop can be cultivated using only rainwater but summer crop cannot be cultivated without irrigation facilities (Kurosaki 2008). The Myanmar government has actively introduced irrigation facilities since 1992 (Matsuda 2009), and thus the phenological cycle with two peaks appeared in the area where irrigation facilities were introduced. In peak analysis, the number of maxima in the phenological cycle was counted to distinguish between single-cropping and double-cropping. This method is simple but widely used as the phenological index technique in remote sensing community (Tao et al. 2017). A time-frequency analysis known as HHT was also applied to detect changes in the cropping system. The HHT is designed for analyzing nonlinear and nonstationary data, while Fourier analysis only give a meaningful interpretation to linear and stationary processes. Despite its wide range of applications, there are few studies by using remote sensing dataset (Chen et al. 2012;Son et al. 2017). The timely fluctuated data of NDVI was decomposed into five IMFs by EMD (Fig. 4). The first IMF (C 1 ) showed clearer phenological cycles than originals in Fig. 3 by noise removal. Then, the instantaneous frequency defined using the Hilbert-Huang transform showed the local phase change (Fig. 5). The instantaneous frequency of double-cropping had higher spikes than singlecropping. In this study, the accumulated annual frequency of 10 was used as the threshold value to distinguish the ricecropping system (Fig. 6). The time change of the frequency could be clearly detect by applying HHT, although the conventional Fourier analysis cannot provide information on the time change of frequency.
The distribution patterns of cropping intensity were similar between two analytical methods (Fig. 7). This indicates the robustness of each analytical method. The rapid 1 3 increase in cropping intensities occurred in the 1990s. It coincided with the time when the irrigation facility was introduced. As a result, the sown area of rice also increased at the same time (Fig. 8). There was a relationship between the cropping intensity and sown area. According to agricultural statistics, the sown area of rice increased 1.65 times 1 3 during 1988 − 2019. This value was slightly larger than the cropping intensities in the 2010s by using peak analysis (1.386 ± 0.483) and time-frequency analysis (1.422 ± 0.499). One of the reasons for the difference is due to a mismatch in agricultural lands. This study was based on the UNEP land cover map in 2000 with a spatial resolution of 1 km. However, land cover in Southeast Asia changed significantly in the past decades, especially for agricultural lands. The influence of land cover change is ignored in this study. Moreover, other land covers, such as scrubland and forests, with one phenological cycle within a year (Fig. 2) are also included as the mixed pixel, the so-called mixel. It is difficult to extract information only on agricultural lands from the satellite data with a spatial resolution of 4 km. Using higher resolution satellite data, such as MODIS, would be desired to get more accurate estimates. Son et al. (2017) assessed the rice-cropping systems with very high accuracy in the Mekong Delta, Vietnam from the MODIS time-series data using EMD. Besides, the cropping intensity showed greater fluctuations than the sown area of rice (Fig. 8). It is desirable to use smooth time-series data in signal processing analysis. We applied the noise removal by the Savitzky-Golay filter, but it might not have been enough. There were over-or under-corrections by pixel even with the same parameters. The improvement of noise removal method for signal processing would be an issue for the future work.

Land and water resource managements
The archive of satellite dataset provided information on the historical changes in cropping system. Double-cropping system was increased along the main rivers in the Delta (Fig. 7). The information can help us understand where land and water resources are still available for constant development, improving cost-effectiveness. Despite favorable agricultural conditions in the Delta, the profitability of rice production is still lower than that in neighboring countries (Coelli et al. 2002;Khan et al. 2016) because the single-cropping is still predominant in Myanmar. Kangyidaunt is one of the townships where rice production has greatly improved because of the introduction of irrigation facilities. Irrigation increases the production capacity and provides benefit to famers in terms of the investment costs per hectare. However, insufficient rice production was reported among farmers immediately after the launching of irrigation facilities, and the lowest rice production occurred at the tail-end farms with difficult access to water resources (Oo et al. 2017). Generically, summer crop can be more productive than monsoon crop when grown under irrigated conditions because of the increased hours of sunshine. However, the peak NDVI values of summer crop were lower than those of monsoon crop in the 1990s and 2000s (Fig. 3a). There might be weaknesses on irrigation management skills of both water users and managers. So, training for efficient utilization of water resources should be promoted. After mastering the management skills, peak NDVI values of summer crop were higher and more stable in the 2010s (Fig. 3a). The government should pay attention to the stability of summer rice production to manage the operation status of irrigation facilities. Educational assistance on irrigation management to farmers is vital for the sustainable success of the agricultural transition process.

Impact of climate change in the Delta
Since 1992, the area of double-cropping consistently increased through the introduction of irrigation facilities (Matsuda 2009), but single-cropping continued to be conducted in the coastal area (Figs. 3, 7). During the monsoon season, crop cultivation was possible even in the coastal area, as the salinity of soil and river water remained low due to the heavy rainfall (~ 3000 mm). However, there was almost no rain during the summer season, and thus seawater traveled many tens of km inland (Driel and Nauta 2014;Sakai et al. 2021). River water contaminated with seawater was not suitable for irrigation. Therefore, crop cultivation was difficult in the coastal area because of the pressures due to both drought and salt stress during the summer season. Besides, the effect of the saline intrusion on crop production was exacerbated by sea-level rise associated with climate change (Rao et al. 2013). Seawater with higher salinity penetrated further inland, and salt accumulation in soils increased in the coastal area (Clarke et al. 2015). Although the monsoon rainfall could remove the accumulated salt in the soil, it could not sufficiently leach the salt deposits when the salinity level exceeded a certain threshold, resulting in a decrease in crop production (Clarke et al. 2015). The negative trends of NDVI max in Labutta Township and other coastal areas were the result of a decrease in crop production due to the exacerbation of salt accumulation (Figs. 3b and 9e, f). According to Fee et al. (2017), the influence of saline intrusion extends through the entire Labutta Township by 2050 (Fig. 9). The saline intrusion has an impact not only for the sustainability of coastal agriculture, but also for the availability of water to meet the growing nonagricultural demand. To recognize the regional changes in the food supply in the face of climate change, the long-term monitoring is required thorough remote sensing.

Conclusion
The large archive of satellite imagery during the period of 1981-2020 provided valuable information for understanding spatiotemporal changes in cropping system and the trends in the Delta of Myanmar. The changes in the cropping system were identified by using the peak analysis and the time-frequency analysis. There were two types of phenological cycles: single-cropping (once a year) and double-cropping systems (twice a year). Although singlecropping system was distributed in mostly in the Delta in the 1980s, summer crop was cultivated since 1992 onwards as a result of the introduction of irrigation facilities. The cropping intensities by the peak analysis and the time-frequency analysis were 1.386 ± 0.483 and 1.422 ± 0.499 in the 2010s, respectively. The double-cropping system enhances crop production efficiency, but efficient expansion strategy of irrigation facilities is important for a sustainable growth of crop production. The potential production differs among locations, depending on the presence of fertile soils and ease of access to the river. The historical maps will provide the information on the practical and cost-effective utilization of remaining land and water resources.
The time-series of NDVI showed negative trends in the coastal area where problems of salt damage were deteriorated. Crop production in the coastal area was constrained by saline intrusion. During the summer season, crop could not be cultivated due to the drought and salt stress. Even during the monsoon season with heavy rain, the accumulated salt in soil reduced the crop production. The salt-affected areas are expected to expand under future climate change scenarios. Government support would be more essential than ever for the sustainable crop production in the Delta.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.