Detection of Corn Phenological Stages with Landsat Satellite Imagery: A Case Study in Ngawi Regency, Indonesia

Ali, Z.,¹ Jaelani, L. M.^1* and Sumargana, L.^2,3

¹Department of Geomatics Engineering, Faculty of Civil, Planning, and Geo-Engineering, Institut Teknologi Sepuluh Nopember (ITS), Surabaya 60111, Indonesia

²Departement of Civil Engineering, Faculty of Civil, Planning, and Geo-Engineering, Institut Teknologi Sepuluh Nopember (ITS), Surabaya 60111, Indonesia

²National Research and Innovation Agency, Gedung B.J. Habibie, Jl. M.H. Thamrin No.8 Jakarta Pusat 10340, Indonesia

*Corresponding Author

Abstract

This study explores how to use the Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM) to combine Landsat and MODIS images and track the growth stages of corn using the vegetation index (NDVI) and the water index (NDWI). The idea was implemented in Ngawi Regency as it is a major contributor to corn production in East Java, Indonesia. Based on spectral index analysis, interviews with local farmers, and data from the Ngawi Food Security and Agriculture Office, the fused images identified two cycles of corn planting and harvesting per year. The first cycle occurred in October–January with the following corn phenology phases: early vegetative (NDVI ~0.3), late vegetative (NDVI ~0.6), early generative (NDVI ~0.8), and late generative (NDVI ~0.7, harvest). The second cycle occurred in February-May/June with the following corn phenology phases: early vegetative (NDVI ~0.3-0.4), late vegetative (NDVI ~0.6-0.7), early generative (NDVI ~0~0.8), and late generative (NDVI ~0.5-0.7, harvest). The results showed that the corn area in May 2019, June 2022, and May 2023 was mostly in the early generative and late generative phases, respectively, covering about 219.15, 178.74, and 236.43 ha from Landsat imagery and 200.88, 180.18, and 205.92 ha from ESTARFM fusion imagery. Despite the differences between Landsat images and ESTARFM fusion images, in detecting the total area of corns, the ESTARFM model successfully predict the corn phenology pattern in Ngawi.

Keywords: Corn, ESTARFM, Ngawi Regency, NDVI, Phenology, Sustainability

1. Introduction

Corn is a versatile and seasonal crop that can complete its life cycle in 80-150 days [1]. The social economy relies on it for human consumption, animal feed, and industrial materials. Along with rice and wheat, corn is also an important source of carbohydrates in Indonesia [2]. Ngawi Regency is one of the largest rice producers in East Java, with a total production of 818.62 thousand tons of Miled Dry Grain (MDG) [3]. The area of Ngawi Regency is 1,394.74 km², of which about 40%, or about 506.6 km², are rice fields. Furthermore, Ngawi Regency contributes to the national corn plantation industry. According to data from the Central Bureau of Statistics (BPS), in 2021, 259,339 tons of corn were harvested from a plantation area of 33,948 ha (hectare) with a productivity of 7.64 tons/ha [4].

Remote sensing is a convenient and inexpensive method to determine the phenological characteristics of cropland, which describes the different stages of a plant's life cycle [5]. These characteristics are often used to assess phenological variation and predict crop yields [6]. Satellite-derived vegetation indices [7] have been used to study different phenological transition phases. A common method used in previous studies is using Normalized Difference Vegetation Index (NDVI) spectral indices from time-series satellite imagery, such as MODIS, to identify crop growth [8]. However, MODIS imagery has a low spatial resolution (≥ 250 m), making it unsuitable for detecting small farms (< 2 ha). Satellite imagery with medium or high spatial resolution (≤ 30 m) is also used to solve this issue.

However, medium- or high-spatial resolution sensors such as Landsat are unsuitable for detecting plant growth and development stages due to their low temporal resolution (16 days) [9]. A previous study attempted to overcome this problem by combining medium/high spatial resolution and high temporal resolution images with different acquisition dates using the Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM) method. ESTARFM is an evolution of Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) that aims to improve prediction accuracy in heterogeneous environments. It often uses images with high temporal resolution, such as MODIS, and images with medium or high spatial resolution, such as Landsat and Sentinel [10]. The temporal resolution produced by ESTARFM, which follows the daily resolution of MODIS, is superior to the combined resolution of Landsat 8 and Landsat 9, which achieves a tandem temporal resolution of 8 days (https://www.usgs.gov/index.php/landsat-missions/ news/landsat-update-special-issue-landsat-9).

This study aims to identify the growth stage of corn plants based on NDVI and Normalized Difference Water Index (NDWI) spectral indices, which will be achieved through data fusion of MODIS and Landsat 8-9 satellite images, as well as Landsat images in Ngawi Regency using the ESTARFM method. To support this study, data validation was conducted using spectral indices obtained from spectrometer measurements, coordinate points of corn samples, and photos of corn in the field, and correlation tests of accuracy using spectral indices from Landsat image processing, MODIS, and ESTARFM image fusion results. In addition, WorldView-2 satellite imagery was used to classify corn and non-corn land to centralize the study area. The results of this study are expected to provide information on corn planting patterns and the distribution of corn growth stages in Rejuno Village, Ngawi Regency. Hopefully, this research can be used to predict corn production based on crown growth stages and planting patterns.

2. Methodology

2.1 Study Area

Rejuno Village is one of the villages in Karangjati District, Ngawi Regency, as shown in Figure 1. Geographically, Rejuno Village is located at 111°38'11.976"E and 7°26'9.46"S. Rejuno Village is located on the east-north side of Ngawi Regency and is characterized by hills. Rejuno Village is directly bordered by several areas, with the north side bordering Bojonegoro Regency, the west side with Plosolor Village, the east side, and the south side directly bordering Rejomulyo Village. Rejuno Village has an area of 1066.60 ha.

2.2 Data

Following the workflow that can be seen in Figure 2, this study used MODIS Nadir BRDF-Adjusted Reflectance (NBAR; MCD43A4; https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MCD43A4) satellite image data at 500 m resolution and Landsat 8-9 Surface Reflectance Collection 2 Tier 1 satellite images to analyze the growth of corn plants in the study area. Image data were obtained from the Google Earth Engine platform ( https://earthengine.google.com/") at acquisition times from April 2013 to June 2023. MODIS-MCD43A4 was selected due to its daily acquisition frequency and the nadir and atmospheric corrections performed [11], despite its low resolution (500 meters). Meanwhile, Landsat 8-9 Surface Reflectance Collection 2 Tier 1 was chosen because of its good spatial resolution (30 meters) and atmospheric correction [12] but low temporal resolution (16 days). The ESTARFM method was used to perform image fusion to produce better spatiotemporal images to utilize its advantages. In addition, Worldview-2 satellite imagery was used for land classification of corn, sugarcane, and others. The data was acquired on May 25, 2023. This study also used vector data of administrative boundaries of Ngawi Regency obtained from the Geospatial Information Agency (BIG), while hamlet boundaries in Rejuno Village, Karangjati Subdistrict, Ngawi Regency were obtained from the Central Bureau of Statistics (BPS) to determine the study area. Spectral data of corn samples of three varieties (NK-Perkasa, NK-212, NK-Sableng, and NK-Sumo) were obtained through spectrometer measurements in Rejuno Village, Ngawi Regency, and processed to obtain validation values of spectral indices (NDVI and NDWI). Corn planting realization data from Ngawi Regency Food Security and Agriculture Office validated corn planting patterns in Ngawi Regency. This study used Google Earth Engine (GEE), Spyder, and ArcGIS software for data analysis.

2.2.1 Determination of the cornfield using the support vector machine classification method

WorldView-2 satellite image data was used to classify corn, sugarcane, and neither field using the support vector machine method. The class was determined to avoid potential overfitting between corn and sugarcane fields. The classification used training data that was manually digitized based on the location of the corn samples and based on elements of image interpretation.

Figure 1: Overview of the study area in Rejuno Village, Karangjati District, Ngawi Regency

Table 1: Cohen's Kappa value

The result was 51 data for each class, for a total of 153 data. The data was then divided into two compositions, training data and validation data, with a ratio of 70:30. This resulted in 105 training data and 48 validation data. Once the land classification is determined, the results must be accurate. In this research, Cohen's Kappa was used to calculate accuracy. Landis and Koch Cohen's Kappa value was categorized into several classes [13], as shown in Table 1.

Figure 2: Workflow of research

2.2.2 Corn growth stage determination using NDVI and NDWI spectral indices

The growth stage of corn in this study was determined using the NDVI and NDWI spectral indices. These indices were obtained by processing Landsat 8 and Landsat 9 imagery. These images were used because corn fields are relatively narrower. Therefore, MODIS images are less appropriate when used to determine the growth phase of corn in narrow fields due to their pixel size of 250-1000 meters. The Landsat image data used were acquired from April 2013 to June 2023 with path 119 and row 65 from the Google Earth Engine platform. NDVI is a common and standard vegetation index used in remote sensing for vegetation monitoring [14]. NDVI is defined in equation 1.

Equation 1

Where:

ρNIR = Near Infrared band reflectance

ρRED = Red band reflectance

NDVI produces good results in chlorophyll processing because NDVI is based on the value between the red band's maximum absorption and the infrared wave's maximum reflectance, which is influenced by the leaf structure [15] NDVI has a range of values from -1 to +1, with values less than 0.2 representing less-vegetation, clouds, and water objects. While the value range of 0.2 to 0.8 is an area of vegetation that is not dense to dense [16]. NDWI is an index used to estimate water potential in vegetation canopies using remote sensing spectral data [17], one of which is the corn plant [18]. NDWI can measure plant physiological stress from leaf changes or water loss [19]. NDWI uses satellite imagery's green and near-infrared (NIR) reflectance bands. Two versions of NDWI are commonly used, the McFeeters NDWI and the Gao NDWI. This study uses the McFeeters version of NDWI, which combines the green and near-infrared (NIR) bands [20].

The NDWI index in equation 2 represents water bodies and measures their turbidity levels [16].

Equation 2

Where:

ρGREEN = Green band reflectance

ρNIR = Near Infrared band reflectance

The determination of the growth stage of corn in this study was obtained by extracting NDVI and NDWI spectral index values using 3 ROIs with a size of 90 x 90 meters or equivalent to the size of 3x3 pixels of Landsat satellite images as sample points in corn fields of NK-Perkasa, NK-Sableng, and NK-Sumo varieties. This sampling was done to obtain values in areas processed by cloud masking on Landsat images. The extracted spectral index values were calculated as median values and determined in a time series graph from April 2013 to June 2023 to identify the growth stage of corn and determine the threshold value of the corn growth stage.

2.3 ESTARFM Framework

The ESTARFM algorithm was developed bya href="#[5]"> [5] to overcome the shortcomings of the STARFM method, namely the inaccurate prediction of surface reflectance in heterogeneous landscapes. ESTARFM was developed further by a href="#[21]">[21] to operate in the Google Earth Engine environment based on the Python API. As a result, algorithms running in Google Earth Engine can operate on a more extensive time series and do not require large local memory. ESTARFM uses two pairs of MODIS-Landsat images acquired on different dates (t1 and t2) to produce images of medium spatiotemporal resolution and is effective for heterogeneous areas. The ESTARFM algorithm is divided into four steps, selecting similar pixels with a neighborhood approach, calculating the weights of similar pixels, calculating the conversion coefficient, and calculating the reflectance at the center of the pixel a href="#[22]">[22].

2.4 Accuracy Assessment

The accuracy assessment in this study used regression analysis and Pearson correlation to evaluate the accuracy and interrelationship between spectral indices from Landsat and MODIS satellite image processing, as well as image fusion results using the ESTARFM method. The correlation coefficient between variables X₁and Y and X₂ and Y is determined based on the Carl Pearson correlation as defined in equation 3 [23]:

Equation 3

Where:

r_Xi,Y = Correlation coefficient

n = Sample number

X = X variable score

Y = Y variable score

∑X = Sum of X variable scores

∑Y = Sum of Y variable scores

∑X²= Sum of squares of X variable scores

∑Y² = Sum of squares of Y variable Scores

The Pearson correlation (r) has a range of values from -1 to 1. A positive correlation indicates a directly proportional relationship, while a negative value indicates an inversely proportional relationship. Evan [24] has provided guidelines for interpreting the Pearson correlation coefficient in several classes, as shown in Table 2.

3. Results and Discussion

3.1 Corn Spectral Library Results

In this study, direct field measurements were made using the Ocean Optics Spectrometer to obtain reflectance values of corn plants. This value used to reference and validate the spectral index processing results of Landsat satellite imagery, MODIS, and ESTARFM image fusion. Corn varieties NK-Perkasa, NK-Sableng, and NK-Sumo were obtained in the late generative phase on the measurement date of February 2-5, 2023 in Rejuno Village.

Table 2: Pearson's correlation coefficient

Figure 3: Reflectance of late generative corn varieties

Table 3: Spectral indices at late generative phase from spectrometer measurement

Table 4: Confusion matrix of corn, sugarcane, and others land cover classification

The NDVI and NDWI spectral index values were calculated from the average measured reflectance by spectrometer (Figure 3) following selected Landsat band range (Red, Green, and NIR) using Equations 1 and 2. The result for three varieties is presented in Table 3.

3.2 Classification of Corn Crop fields

The results obtained from the land classification process using the Support Vector Machine (SVM) method (Figure 4) have a kappa accuracy of 81.18% and an overall accuracy of 87.50%. These results fall into the "almost perfect" category [13]. Table 4 describes the Confusion Matrix value for each corn, sugarcane, and others class. Once good accuracy was achieved, the satellite imagery was clipped according to the administrative boundary vector data of Rejuno Village. The Region of Interest (ROI) was then created with 3x3 Landsat pixels or 90⨯90 meters for NK-Perkasa, NK-Sableng, and NK-Sumo corn varieties as presented in Figure 5. This aimed to obtain the spectral index value when cloud cover and shadows mask the satellite image.

3.3 Identification of Corn Growth Stages over Ten Years

Detection of corn growth phase using two spectral indices (NDVI and NDWI). The temporal changes were identified using Landsat and MODIS images from April 2013 to June 2023. However, MODIS has too large a spatial resolution (500 meters), so the corn growth phase was detected using Landsat-8 and Landsat-9 imagery, respectively. Figure 5 shows the temporal changes of the two spectral indices from April 2013 to June 2023.

Figure 4: Landcover classification

Figure 5: ROI of NK-Perkasa, NK-Sableng, and NK-Sumo corn fields with a size of 90&#215;90 meters

Figure 6: Graphs of NDVI and NDWI using Landsat April 2013 – June 2023 in the ROI of corn varieties
(a) NK-Perkasa (b) NK-Sableng (c) NK-Sumo

Figure 6 shows the correspondence of spectral index patterns between the non-corn column (before October 2018) and the corn column (after October 2018). Thus, the detection of the growth phase of corn requires supporting data from historical satellite imagery obtained from Google Earth Pro, Google Street Map, and data from the Ngawi Food Security and Agriculture Office to determine historical land cover. Despite the similarity of the pattern of increase and decrease in spectral index between corn and non-corn plots, some differences can be seen in the figure. One notable difference is the consistency of the lowest NDVI values in the corn column, which had the lowest values of ~0.2 - ~0.3 in October, while the non-corn column had the highest minimum NDVI values of ~0.5 - ~0.6 in October 2016, respectively. Thus, the difference compared to the corn column is quite large.

Given the inversely proportional spectral index patterns of NDVI and NDWI, one of the indices, such as NDVI, can be used to identify the growth stage of corn. When the NDVI from the October 2018 ESTARFM fusion image is compared to each corn variety, the results are shown in Figure 7. Figure 6 shows that the NDVI values produced by the ESTARFM image fusion have a data gap. This results in the loss of NDVI values used to analyze the growth phase of corn. According to local farmers, the values obtained are from June to October when corn is no longer planted. Therefore, the NDVI values of the Landsat satellite images were used to determine the threshold for the growth stage of corn. The threshold values are shown in Table 5.

3.4 Corn Crop Phenology Model

The NDVI and NDWI values from Landsat-8 and Landsat-9 images were used to determine the phenology of corn plants per variety each month throughout the year. The phenology of corn plants was analyzed based on the phase transitions experienced by corn during each month. These transitions were determined based on the threshold values for each phase, as shown in Table 5. From these values, the phenology of each corn variety NK-Perkasa, NK-Sableng, and NK-Sumo was modeled as shown in Tables 6 to 8. Based on Tables 6 to 8, a similar growth pattern was found for each corn variety. In 2018, the planting phase for varieties NK-Perkasa and NK-Sableng started in October, which was characterized by the lowest NDVI values of 0.300 and 0.302, respectively. These values represent the early vegetative stage. On the 18th day after planting, the plants entered the late vegetative phase, indicated by an increase in NDVI values to 0.474 and 0.473. This phase lasted for approximately 32 days. Due to significant cloud cover in December 2018, NDVI values could not be collected. During this month, the corn plants had entered the early generative phase, which indicates the highest NDVI values during the growing season. This phase lasted for 25 days, after which the plants entered the late generative phase in January 2019, characterized by a decrease in NDVI values for 34 days for NK-Perkasa and 25 days for NK-Sableng. Therefore, by the end of January or early February 2019, the corn plants had reached the physiological maturity stage and were ready to be harvested. In 2019, the planting cycle started again with a pattern known as Rainy Season 2: planting from February to June. After a post-harvest rest period, corn planting resumed mid to late February. The NDVI values in the early vegetative phase in this pattern were not as low as in the Rainy Season 1 planting pattern, namely 0.744 for NK-Perkasa and 0.844 for NK-Sableng, because the rainy season conditions were still ongoing. This phase lasted for 18 days after planting. Then, the entry into the late vegetative phase, indicated by increasing NDVI values, lasted for about 32 days. Thus, in April 2019, the corn had entered the early generative phase, indicated by the highest NDVI values of 0.884 and 0.789. The early generative phase lasted about 25 days. Then, the NDVI values decreased, namely 0.293 and 0.260, indicating that the plants entered the late generative phase until harvest time. From 2019 onwards, the corn planting pattern on the sample plots of corn varieties NK-Perkasa and NK-Sableng did not change. However, extensive cloud cover during the rainy months resulted in some NDVI values that could not be identified. The monthly columns on the phenology calendar that indicate this condition are colored light gray. Corn fields in Rejuno Village are planted twice a year during the rainy season (Mbah Wo, personal communication). During the dry season, the land is left fallow and not planted. In addition, the new corn variety NK-Sumo has been planted since October 2022, and the new variety has been successfully tested.

Figure 7: Landsat NDVI (April 2013 – June 2023) in the ROI of corn varieties NK-Perkasa, NK-Sableng, and NK-Sumo

Table 5Threshold value of corn phenology

Table 6: Monthly median values of NDVI and corn phenology in the NK-Perkasa's ROI for 2018-2023 using Landsat 8-9 imagery

Table 7: Monthly median values of NDVI and corn phenology in the NK-Sableng's ROI for 2018-2023 using Landsat 8-9 imagery

Table 8: Monthly median values of NDVI and corn phenology in the NK-Sumo's ROI for 2018-2023 using Landsat 8-9 imagery

Figure 8: Distribution of corn phenology in May 2019
(a) Landsat imagery (b) ESTARFM fusion

Figure 9: Distribution of corn phenology in June 2022
(a) Landsat imagery (b) ESTARFM

3.6 Corn Growth Stage Distribution Identification

Based on the NDVI and NDWI thresholds, as well as the corn cropping pattern in the Rejuno Village obtained from Landsat imagery and ESTARFM fusion results, the distribution of corn growth phases in the Rejuno Village, especially in the late generative phase, is predicted to be around May to June. The late generative phase was chosen based on spectral validation of corn plant samples obtained during this phase. Visualizations of the distribution of the late generative phase of corn plant growth in May 2019, June 2022, and May 2023 are shown on the maps in Figures 8 to 10, respectively.

Figure 10: Distribution of corn phenology in May 2023
(a) Landsat imagery (b) ESTARFM fusion

Table 9: Corn growth phase distribution area

These years are images that can be fully formed when corn is in the late generative phase from ESTARFM image fusion processing from 2018 to 2023. Based on the corn growth stage distribution results in Figures 8 to 10. Table 9 illustrate the growth stage on a monthly basis and the area information.

In May 2019, the most dominant phase was the early generative phase with an area of 219.15 ha based on Landsat image processing and 200.88 ha based on ESTAR-FM fusion image processing. This phase was most widespread in the Karangjati forest area and least in the Gurdo Barat Hamlet.

Other phases that were obtained were the late vegetative phase, with an area of 0.36 ha based on Landsat imagery processing and 0.27 ha based on ESTARFM fusion imagery processing, and the late generative phase, with an area of 70.83 ha based on Landsat imagery processing and 85.23 ha based on ESTARFM fusion imagery processing. In June 2022, the most dominant phase was also the early generative phase with a total area of 178.74 ha based on Landsat imagery processing and 180.18 ha based on ESTARFM fusion imagery processing. This phase is spread in Karangjati Forest, Kalisoko Hamlet, Karanganyar Hamlet, Gurdo Timur Hamlet, and Pakulan Hamlet. Other phases that were obtained were the early vegetative phase with an area of 2.43 ha based on Landsat imagery processing and 1.35 ha based on ESTARFM fusion imagery processing, the late vegetative phase with an area of 3.24 ha based on both Landsat imagery processing and ESTARFM fusion imagery processing, and the late generative phase with an area of 128.34 ha based on Landsat imagery processing and 119.52 ha based on ESTARFM fusion imagery processing.

In May 2023, the most dominant phase was the late generative phase with an area of 236.43 ha based on Landsat imagery processing and 205.92 ha based on ESTARFM fusion imagery processing. This phase occurs throughout the study area and is most common in the Karangjati Forest. Other phases that were obtained were the early vegetative phase with an area of 1.17 ha based on Landsat imagery processing and 0.99 ha based on ESTARFM fusion imagery processing, the late vegetative phase with an area of 3.96 ha based on Landsat imagery processing and 2.97 ha based on ESTARFM fusion imagery processing, and the early generative phase with an area of 101.88 ha based on Landsat imagery processing and 79.38 ha based on ESTARFM fusion imagery processing. The difference between ESTARFM and Landsat in the term of area come from the flag (usually cloud flag) from Landsat and/or MODIS. The ESTARFM only produces clear pixels when no flag at the exact same pixel in MODIS or Landsat data.

4. Conclusion

Based on the processing results, most of the corn in May 2019 is early generative. It covers about 219.15 ha from Landsat imagery and 200.88 ha from ESTARFM fusion imagery. In June 2022, the corn phase remained mostly early generative. It covers about 178.74 ha from Landsat imagery and 180.18 ha from ESTARFM fusion imagery. In May 2023, the most abundant corn phase was late generative. It covers about 236.43 ha from Landsat imagery and 205.92 ha from ESTARFM fusion imagery these phases are the most widespread in Karangjati Forest each year. The growth cycle pattern of corn varieties NK-Perkasa, NK-Sableng, and NK-Sumo in Rejuno Village is twice a year, namely in the first rainy season (October-January) and in the second rainy season (February-June). NDVI values average about 0.3 in the early vegetative phase, increase to about 0.6 in the late vegetative phase, peak at about 0.8 in the early generative phase, and decrease to about 0.7 or 0.5 in the late generative phase and harvest period. Corn is not planted during the dry season (July-September).

Future studies should consider utilizing Sentinel-2 satellite imagery for corn crop phenology modelling due to its superior spatial, temporal, and spectral resolution compared to Landsat-8 and Landsat-9. The selection of satellite imagery should inherently ensure clarity and accuracy of data. Additionally, the use of hyperspectral satellites could be beneficial as they offer detailed and comprehensive reflectance spectra at various wavelengths, thereby enabling a more detailed determination of crop characteristics.

Acknowledgments

The authors gratefully acknowledge the National Research and Innovation Agency (BRIN); PT. Aria Agri Indonesia and the Office of Food Security and Agriculture in Ngawi Regency, which assisted in providing data for this research. The authors thank Ty C. Nietupski for making the ESTARFM script code based on Google Earth Engine available to the community as an open source.

References

[1] Cheesman, E. E., (1973). Tropical Crops: Monocotyledons Vols. 1 and 2 By J. W. Purseglove London: Longman (1972), 334 and 273. Experimental Agriculture , Vol. 9(3), 287–287. < a href="https://doi.org/10.1017/S0014479700005822"> https://doi.org/10.1017/S0014479700005822.

[2] Kemendag, (2021). Profil Komoditas Jagung. [Online]. Available: https://ews.kemendag.go.id/sp2kplanding/assets/pdf/120116_ANK_PKM_DSK_Jagung.pdf. [Accessed Oct. 28, 2022].

[3] Bappeda Jatim, (2023). Jawa Timur Kembali Sebagai Penghasil Padi Terbesar Di Indonesia. [Online]. Available: https://bappeda.jatimprov.go.id/2022/01/26/jawa-timur-kembali-sebagai-penghasil-padi-terbesar-di-indonesia/". [Accessed Jan. 05, 2023].

[4] BPS, (2022). Kabupaten Ngawi Dalam Angka. Ngawi: BPS.

[5] Zhu, X., Chen, J., Gao, F., Chen, X. and Masek, J. G., (2010). An Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model for Complex Heterogeneous Regions. Remote Sensing Environment, Vol. 114(11), 2610–2623. https://doi.org/10.1016/j.rse.2010.05.032.

[6] Tang, J., Zeng, J., Zhang, Q., Zhang, R., Leng, S., Zeng, Y., Shui, W., Xu, Z. and Wang, Q., 2020. Self-Adapting Extraction of Cropland Phenological Transitions of Rotation Agroecosystems Using Dynamically Fused NDVI Images. International Journal of Biometeorology, Vol. 64(8). 1273–1283. https://doi.org/10.1007/s00484-020-01904-1.

[7] Helman, D., (2018). Land Surface Phenology: What Do We Really ‘See’ from Space?, Science of the Total Environment, Vol. 618, 665–673. https://doi.org/10.1016/j.scitotenv.2017.07.237.

[8] Pan, Z., Huang, J., Zhou, Q., Wang, L., Cheng, Y., Zhang, H. K., Blackburn, G. A., Yan, J. and Liu, J. (2015). Mapping Crop Phenology using NDVI Time-Series Derived from HJ-1 A/B Data. International Journal of Applied Earth Observation and Geoinformation , Vol. 34(1). 188–197. https://doi.org/10.1016/j.jag.2014.08.011.

[9] Zhang, X., Wang, J., Henebry, G. M. and Gao, F., (2020). Development and Evaluation of a New Algorithm for Detecting 30 m Land Surface Phenology from VIIRS and HLS Time Series, ISPRS Journal of Photogrammetry and Remote Sensing , Vol. 161, 37–51. https://doi.org/10.1016/j.isprsjprs.2020.01.012.

[10] Zhou, X., Wang, P., Tansey, K., Zhang, S., Li, H. and Tian, H., (2020). Reconstruction of Time Series Leaf Area Index for Improving Wheat Yield Estimates at Field Scales by Fusion of Sentinel-2, -3 and MODIS Imagery. Comput Electron Agric, Vol. 177. https://doi.org/10.1016/j.compag.2020.105692.

[11] Wang, Z. and Schaaf, C. B., (2023). MODIS/Terra+Aqua BRDF/Albedo Nadir BRDF Adjusted Ref Daily L3 Global - 500m V061 [Dataset]. NASA EOSDIS Land Processes DAAC . [Online]. Available: https://lpdaac.usgs.gov/products/mcd43a4v061/". [Accessed: May 01, 2023].

[12] Earth Resources Observation and Science (EROS) Center, (2023). USGS EROS Archive - Landsat Archives - Landsat 8-9 OLI/TIRS Collection 2 Level-2 Science Products,” SCIENCE. [Online]. Available: https://www.usgs.gov/centers/eros/science/usgs-eros-archive-landsat-archives-landsat-8-9-olitirs-collection-2-level-2#overview. [Accessed May 01, 2023].

[13] Landis, J. R. and Koch, G. G., (1977). The Measurement of Observer Agreement for Categorical Data. Biometrics, Vol. 33(1). https://doi.org/10.2307/252931.

[14] Bhandari, A. K., Kumar, A. and Singh, G. K., (2012). Feature Extraction using Normalized Difference Vegetation Index (NDVI): A Case Study of Jabalpur City. Procedia Technology, Vol. 6, 612–621. https://doi.org/10.1016/j.protcy.2012.10.074.

[15] Alharbi, S. W. R., Arnall, B. and Zhang, H., (2019). Prediction of Maize (Zea mays L.) Population Using Normalized-Difference Vegetative Index (NDVI) and Coefficient of Variation (CV). Journal of Plant Nutrition , Vol. 42(6), 673–679. https://doi.org/10.1080/01904167.2019.1568465.

[16] Hardianto, A., Dewi, P. U., Feriansyah, T., Sari, N. F. S. and Rifiana, N. S., (2021). Pemanfaatan Citra Landsat 8 Dalam Mengidentifikasi Nilai Indeks Kerapatan Vegetasi (NDVI) Tahun 2013 dan 2019 (Area Studi: Kota Bandar Lampung). Jurnal Geosains dan Remote Sensing, Vol. 2(1), 8–15. https://doi.org/10.23960/jgrs.2021.v2i1.38.

[17] Gao, B., (1996). NDWI-A Normalized Difference Water Index for Remote Sensing of Vegetation Liquid Water from Space. Remote Sensing of Environment , Vol. 58(3), 257–266, https://doi.org/10.1016/S0034-4257(96)00067-3.

[18] Kerner, H. R., Sahajpal, R., Pai, D. B., Skakun, S., Puricelli, E., Hosseini, M., Meyer, S. and Becker-Reshef, I., (2022). Phenological Normalization Can Improve In-Season Classification of Maize and Soybean: A Case Study in the Central US Corn Belt. Science of Remote Sensing, Vol. 6. https://doi.org/10.1016/j.srs.2022.100059.

[19] Alvino, F. C. G., Aleman, C. C., Filgueiras, R., Althoff, D. and da Cunha, F. F., (2020). Vegetation Indices for Irrigated Corn Monitoring. Engenharia Agrícola, Vol. 40(3), 322–333. https://doi.org/10.1590/1809-4430-eng.agric.v40n3p322-333/2020.

[20] McFeeters, S. K., (1996). The Use of the Normalized Difference Water Index (NDWI) in the Delineation of Open Water Features. International Journal of Remote Sensing , Vol. 17(7), 1425–1432. https://doi.org/10.1080/01431169608948714.

[21] Nietupski, T., (2021). GoogleEarthEngine _ImageFusion. Mendeley Data, Corvallis: 1. https://doi.org/10.17632/bcbptkrbsg.1.

[22] Liu, M., Liu, X., Dong, X., Zhao, B., Zou, X., Wu, L., and Wei, H. (2020). An Improved Spatiotemporal Data Fusion Method Using Surface Heterogeneity Information Based on ESTARFM. Remote Sensing, Vol. 12(21), 3673. https://doi.org/10.3390/rs12213673.

[23] Abdullah, Yusrizal, Z. and Syamsuddin, S., (2018). Analisis suhu permukaan laut dan klorofil-a untuk mengetahui potensi ikan di dalam fishing ground purse seine menggunakan data citra satelit aqua modis level 3. Jurnal Kelautan dan Perikanan Terapan JKPT, Vol. 1(1), 38–45.

[24] Evans, J. D., (1996). Straightforward Statistics for The Behavioral Sciences. Thomson Brooks/Cole Publishing Co.