Evaluation of Geographical and Annual Changes in Rice Planting Patterns with Satellite Images in Flood-Prone Area of Pampanga River Basin, the Philippines

1. Introduction

Floods are some of the most devastating disasters, causing significant damage to crops around the world. The FAO has announced that damages to crops and livestock caused by flooding amounted to USD 21 billion, comprising 19% of the total loss caused by all disasters from 2008 to 2018 [1]. In order to reduce the associated damages, cultivation patterns are varied depending on the flood risk in the floodplains of Southeast Asia [2,3]. Information on cultivation patterns is extremely important for the management and evaluation of flood damage [4]. At present, such information is mainly collected through field surveys. However, field surveys mostly fail to provide the geographical distribution of cultivation patterns, as they are conducted on a points basis. Moreover, surveys ordinarily consider a certain management period or growth stage, for example, the planting period or heading stage of cereals. These information constraints limit the ability to evaluate flood damage based on the crop growth stage or inundation depth.

The Pampanga River Basin is located in the central part of Luzon Island in the Philippines, and is one of the largest floodplains and rice-producing areas in Southeast Asia. Floods are frequently caused by monsoons and typhoons during the rainy season, affecting crops in various ways almost every year. For example, during Typhoon “ULYSSES”, 4490 hectares of farmland were damaged throughout the entire province of Pampanga, with agricultural damage amounting to about 200 million pesos [5]. To avoid these damages, farmers have developed various cultivation patterns including fish raising and vegetable cultivation [6]; however, the geographic distribution and annual change of cultivation patterns has not yet been analyzed. Such an analysis would provide important information, helping to evaluate how farmers have adapted to flooding and to develop countermeasures for the future.

Remote sensing using satellite images has been considered as an alternative evaluation method to field surveys, due to the versatility and wide applicability of the obtained information. Vegetation indices or leaf area index (LAI) obtained from satellite images have been employed for the evaluation of crop growth [7,8], growth stage [9,10], and cultivation management [11]. The coarse resolution and cloud disturbances are major constraints for the utilization of satellite images [12]; however, previous studies have suggested that the use of statistical procedures on long-term time-series data allow for the detection of changes in cultivation patterns. Iwahashi et al. (2021) used time-series data from the MODIS LAI product and revealed conversion from late-matured to early-matured cultivars and extension of dry season rice cropping in Cambodia [13]. Zhao et al. (2022) used long-term time-series data from the MODIS EVI product and revealed the effect of global warming on the phenological changes of wheat in China [14]. Based on the above information, we evaluated annual and geographical differences in cultivation patterns in the Pampanga River Basin, the Philippines by analyzing long-term time-series MODIS data.

2. Materials and Methods

2.1. Study area

The study was conducted in Candaba municipality, located in the floodplain of the Pampanga River Basin (Figure 1). The area is characterized by a tropical monsoon climate. The rainy season is generally from June to October, during which typhoons of various sizes hit the area. In the rainy season, the several rivers flowing into the area often lead to floods, especially in typhoon events. The floods transform some areas, through inundation, into a complex comprising variably saturated wetland, shallow swamps, and deeper freshwater ponds [15].

Candaba belongs to Region Ⅲ (Central Luzon), and is one of the largest rice-producing area in the Philippines [16]. Rice is the dominant crop in Candaba, but other modes of cultivation, including fish and vegetable cultivation, have also been reported [6]. The cultivation period and times for rice vary, depending on flood and dry cultivation.

2.2. LAI data

MODIS data (MCD15A2H) were downloaded for the 16 years from 2007 to 2022. MCD15A2H provides Leaf Area Index (LAI) values with 500 m pixel size every 8 days. In this study, the 2013 version of GLCNMO, a raster data set that contains global land-cover information based on MODIS data with 500 m pixel size [17], was used to identify paddy fields. LAI values were extracted from pixels in MCD15A2H which were identified as paddy fields according to GLCNMO.

2.3. Analysis of LAI dynamics

Hierarchical clustering was performed on the extracted LAI data. The cluster analysis was conducted in two stages. In the first stage, hierarchical clustering was performed on the time-series LAI data for 16 years in each pixel (n = 618), in order to classify areas with similar LAI dynamics. In the second stage, hierarchical clustering was performed on the time-series LAI data for one year in each pixel of each area classified in the first stage of hierarchical clustering (n = number of pixels in the area × 15). The one-year data were divided from October of the previous year to September of the target year, in order to include the planting periods of the dry and rainy seasons. Ward’s method was employed to calculate the distance between clusters. The thresholds to determine the clusters were set at 60% and 70% of all Euclidean distances in the first and second stages, respectively.

2.4. Software

MCD15A2H was downloaded using Google Earth Engine [18]. LAI data extraction from MCD15A2H was performed using Quantum Geographic Information System (QGIS version 3.14). Hierarchical clustering was performed using the SciPy library in Python (SciPy version 1.7.1).

3. Results

3.1. Area classification in Candaba

First, hierarchical clustering was performed on the MODIS LAI product (MCD15A2H) for the 16 years of the study period (2007–2022). The first stage clearly classified Candaba into four areas: The southern area (Area 1), the northern area (Area 2), the lower terrain area (Area 3), and the eastern and western area (Area 4); see Figure 2.

3.2. Annual LAI dynamics for each region

The second stage of hierarchical clustering was performed on the MODIS LAI data to compare the annual LAI dynamics for each area and year.

3.2.1. Area 1 (Southern area)

Area 1 was divided into two clusters. Figure 3 shows the averaged LAI dynamics for each cluster. Both clusters had a large peak during the dry season. Cluster 2 reached its peak earlier than Cluster 1, and the LAI value at the peak was higher than that for Cluster 1. A small peak after the large peak was observed in Cluster 2, and its peak was around the end of the dry season. Figure 4 shows the annual distribution for each cluster. Cluster 1 was frequently observed after 2015, while Cluster 2 was frequently observed before 2014.

3.2.2. Area 2 (Northern area)

Area 2 was divided into three clusters. Figure 5 shows the averaged LAI dynamics for each cluster. All clusters had three peaks throughout the year. Cluster 2 reached its first peak earlier than the other two clusters. Meanwhile, Cluster 3 reached its first peak later than the other two clusters, but its LAI value was significantly higher than those of the other two clusters. Figure 6 shows the annual distribution for each cluster. The northern part of the area was occupied by Cluster 1 and Cluster 3 from 2007 to 2022, while the southern part of the area transitioned from Cluster 1 to Cluster 2 around 2014.

3.2.3. Area 3 (Lower terrain area)

Area 3 was divided into three clusters. Figure 7 shows the averaged LAI dynamics for each cluster. All clusters had a large peak during the dry season (first peak) and an additional peak just before the rainy season (second peak). Cluster 3 reached its first peak earlier than the other clusters. Cluster 2 was the latest and had a small peak before its first peak. The LAI value at the second peak of Cluster 3 was higher than that of Cluster 2. Figure 8 shows the annual distribution for each cluster. Cluster 1 occupied the whole area until 2013, while Cluster 3 was more frequently observed after 2014. Cluster 2 occupied the entire area in 2015 and 2020.

3.2.4. Area 4 (Eastern and western area)

Area 4 was divided into three clusters. Figure 9 shows the averaged LAI dynamics for each cluster. The LAI value of Cluster 1 fluctuated between 1 and 2 throughout the year. Cluster 2 had two peaks: August–September and December–January. Cluster 3 had a single peak, in February–March. Figure 10 shows the annual distribution for each cluster. The western area was occupied by Cluster 1 and Cluster 3 from 2007 to 2021. On the other hand, except for the years between 2010 to 2012, most of the eastern area was occupied by Cluster 2.

4. Discussion

In this study, in order to reveal the differences in the annual LAI dynamics in Candaba, a two-stage clustering analysis was applied: the first stage involved classification of area, while the second stage involved classification of annual LAI dynamics. Direct clustering analysis of annual LAI dynamics without area classification did not yield obvious differences in the LAI dynamics pattern and its geographical distribution (data not shown). As the LAI data for a pixel may include outliers and/or errors, the observation of annual LAI dynamics patterns is often noisy. Increasing the dimensionality by combining 16 years of time-series data allowed for analysis of the geographical differences in LAI dynamics, thus increasing similarity by dividing areas showing similar annual patterns.

In the first stage of classification, Candaba was divided into four areas. As the time-series LAI data were extracted for paddy fields, these areas may reflect local rice cultivation patterns. The classified areas were partly consistent with the geography (Figure 1) and flood areas under the influence of a typhoon (Figure 11). As the timing of cultivation depends on the flood risk in Candaba [6], a map similar to that shown in Figure 2 might be produced by combining the flood risk determined by topography with the experience of farmers. Area 4 presented lower LAI values and obscure peaks, suggesting that rice cultivation was not dominant in this area. This result was expected, as Area 4 occupies relatively higher terrain, where upland crops are sometimes planted (personal communications). Further analysis utilizing finer resolution satellite data may be required for more precise results.

In the second stage of classification, the annual LAI dynamics in each area were divided into two or three clusters. One of the major differences among the clusters was the time of the first (and largest) peak, which seems to reflect the beginning of planting. An earlier first planting generally produced higher second and third peaks. As farmers wait to plant until the flood risk is low, the distribution of clusters might be governed by the flood condition in the rainy season. In particular, Cluster 2 in Area 3, which widely occupied the area in 2015 and 2020, probably indicates that the rice was damaged and delayed by one month, when compared to the other clusters. Notably, 2015 and 2020 were years that suffered from typhoon damage between November and December. In 2015, typhoon “Nona” passed through the area from December 12^th to 17^th, and 17 barangays were damaged by flooding [19]. In 2020, typhoon “Ulysses” passed through the area from November 8^th to 13^th, and 29 barangays were damaged by flooding [5]. According to the agricultural department of the municipal office of Candaba, these typhoons caused a delay in planting during the dry season (personal communications), consistent with the observed LAI dynamics.

Although the largest peak was observed in the dry season, this does not necessarily mean that the area and production of the dry season crop is higher than that in the rainy season. In fact, rainy season crops are also cultivated in areas with low flood risk [6]. The reason why the peak in the dry season appeared to be stronger than that in the rainy season may be that satellite observations are more likely to be disturbed by clouds in the rainy season [20]. Deeply submerged water also disturbs satellite observations [21]. Furthermore, the cultivation in the dry season is more likely to be aligned due to inundation at the end of rainy season, which may be another reason for the enhanced peak observed in the dry season. Melendres (2014) has reported that rice is usually cultivated twice a year in Candaba, but only once in the heart of Candaba Swamp due to flooding in the rainy season [22]. However, the annual LAI dynamics showed 1, 3, and 2 peaks per year in Areas 1, 2, and 3, respectively. This inconsistency might be partly derived from the inadequacy of field surveys, which would not allow for coverage of the whole area.

The annual changes in cluster distribution shown in Figure 4, Figure 6, Figure 8, and Figure 10 suggest that the rice cultivation patterns changed over the 16 years of the study period. For example, Cluster 2 transitioned to Cluster 1 in Area 1, Cluster 1 transitioned to Cluster 2 in the southern part of Area 2, and Cluster 1 partly transitioned to Cluster 3 in Area 3. The transitions in Areas 2 and 3 indicated earlier planting, which would be suitable for double cropping in the dry season. Meanwhile, the transition in Area 1 was accompanied with delayed planting. Although we do not have a reasonable explanation, flood conditions generally lead to delayed planting. Further field surveys would be required for assessment of the flood condition at the beginning of the dry season.

As mentioned above, the long-term time-series MODIS LAI product indicated that cultivation patterns changed over the study period, and are also influenced by floods. However, the cluster analysis appears to have detected little damage to LAI due to flood except for Cluster 2 in Area 3. The low resolution and accuracy of the data may have reduced the detection power for flood damage, and higher resolution satellite data may be required for evaluation of the flood damage [23]. However, low observation opportunity due to cloud cover restricts the application of visible and near-infrared sensors in the rainy season. Synthetic aperture radar (SAR) sensors are often utilized to observe the ground status in rainy season [24,25]. However, the evaluation of vegetation by SAR generally has low accuracy [26,27], although that of water bodies due to flooding may be acceptable [28]. Accordingly, combining cultivation patterns (assessed by MODIS) and flood area (by SAR) may provide a promising method for the evaluation of flood damage. In order to develop a real-time monitoring system for flood damage in the context of crop production, detection of the start of cultivation is required. Previous studies have reported that SAR can be used to detect the inundation for preparation of rice cultivation in irrigated paddy fields [29,30]. However, natural flooding during rice growing or fallowing may disturb such detection in flood-prone environments. The cultivation patterns evaluated in this study may support such cases.