An integrated method to develop high temporal Landsat NDVI using free and open-source software (FOSS) GIS application

Landsat image series are long-term times series observing the earth’s surface at a relatively fine spatial resolution of 30 m on the visible and infrared channels, which can be used to monitor vegetation dynamics using band ratios such as NDVI. However, due to cloud cover, having Landsat NDVI at a high temporal resolution is impossible. This is especially found in high-humidity regions such as Indonesia, particularly in the rainy season. Whereas some studies, such as long-term phenological observations and detailed land use/land cover change detection, may require Landsat NDVI with a high temporal resolution. One of the higher temporal resolution images with similar spectral characteristics to Landsat is MODIS, which can be used to fill the gap; however, it may have a coarser spatial resolution. Thus, some methods have been developed to blend Landsat and MODIS images to create images with both Landsat spatial detail and MODIS high temporal resolution. Unfortunately, the blending methods are usually not straightforward, where complex steps are required to find and combine several techniques with advanced computational skills. This study, therefore, aims to develop an integrated method to create high temporal Landsat NDVI using free and open-source software (FOSS) GIS application, which is easily accessible and applicable.


Introduction
Observation of biophysical processes such as vegetation phenology and land use/land cover dynamics on the earth's surface can be done using remote sensing technology, such as NDVI [1].Several remote sensing images are free, provided by a range of sensors capable of capturing an object on the earth's surface with variable spatial and temporal resolution.Sensors with a high temporal resolution (i.e., daily acquisition) are very useful for observing the dynamics at the earth's surface to gain insight into its biophysical processes.However, sensors with high temporal resolution mainly work at a relatively coarse spatial resolution that is fit to capture the dynamics of earth surface processes at regional and global scales.Studies at the local scale also require remote sensing images with a fine or relatively detailed spatial resolution.Unfortunately, freely available fine spatial resolution NDVI such as provided by Landsat image [2], usually have a low temporal resolution.In addition, in the case of Indonesia, with very frequent high cloud cover, the temporal resolution is even further decreased [3].This situation may 1291 (2024) 012014 IOP Publishing doi:10.1088/1755-1315/1291/1/012014 2 hamper the local study of earth surface dynamics, requiring observations using Landsat NDVI at a high temporal resolution.To fill the gap, several methods have been introduced to blend remote sending image that has a fine spatial resolution, such as Landsat, with one that has a high temporal resolution, such as MODIS, to generate a synthetic image that has both fine spatial and high temporal resolution at the same time [4].
Available remote sensing image blending methods from the literature mainly need advanced computational procedures.This requires a complex compilation of processes not provided by standard remote sensing and geographic information system (GIS) applications.Therefore, developing an integrated method to blend a fine resolution of Landsat NDVI with a high temporal resolution NDVI such as MODIS NDVI will be helpful for a local study that requires a high temporal Landsat NDVI.To make the integrated method widely available and its code accessible, it should be developed on free and open-source software (FOSS) GIS.At present, some FOSS GIS can be selected to generate high temporal Landsat NDVI [4].However, there is still a limited application based on FOSS GIS that integrates methods to develop such images.Therefore, this study aims to develop an integrated method to generate high temporal Landsat NDVI based on an available image blending approach using FOSS GIS software.

Methods
For this study, we select an area of interest that cover the eastern slope of Mount Merapi and the lower flatland at Klaten and Sleman District in Central Java and Yogyakarta Province, Indonesia, due to being located in the region with very frequent cloud cover [5].

Image blending method
The main principle of the image blending method is to generate a synthetic fine spatial resolution image (f1) using a coarse spatial resolution image but has high temporal resolution (c1) based on pair(s) of cloud-free fine (f0) and coarse (c0) spatial resolution images as the reference.Several image blending methods are available from the literature, such as the spatial and temporal adaptive reflectance fusion model (STARFM) [6], the unmixing-based data fusion techniques (UBDF) [7], the flexible spatiotemporal data fusion (FSDAF) [8] and many more.From those available methods, we select one that the code is available for public use, can be adapted to FOSS GIS, and require minimal input data but has acceptable performance.The STARFM is a proven method [10] that fulfills the abovementioned requirement.Therefore, STARFM was used in this study which was set up using single pair reference images of Landsat 8 OLI NDVI and MODIS NDVI and another pair of the same pair image but different in acquisition date for developing synthetic Landsat 8 OLI NDVI image (using MODIS NDVI) and to validate the result (using Landsat 8 OLI NDVI).The STARFM algorithm works based on weighted functions, assuming that the reflectance changes are consistent and comparable for fine and coarse spatial resolution.Further, STARFM also assumes that the condition of the earth's surface captured by the coarse resolution pixels contains a homogeneous land cover type.Regarding those assumptions, changes that occurred on the earth's surface and were detected from the comparison between the coarse pixels reference (c0) and predictor (c1) images are used as input for weighted functions to generate the synthetic fine resolution pixels (f1).

Dataset acquisition and FOSS GIS selection.
The dataset used in this study consists of Landsat  c1) is one with a composite date of 30 September 2013.We use the FOSS application of R programming, GRASS-GIS, and MODIS Reprojection Tools (MRT) for code compilation and image data processing.R is widely used in many programming applications, while GRASS-GIS has also been applied to many studies.The MRT is a standard application that is used to prepare MODIS images ready to use in standard image format [10].In addition, all the software is crossplatform, and GRASS-GIS also MRT can be integrated into R.

Code compilation
Code compilation was done mainly in R, which covers the whole process of data pre-processing and image blending, which are described below.For this study, we use Rstudio software as an R code editor.

Data pre-processing
Data pre-processing aims to prepare image data for the image blending process.Data pre-processing was applied mainly to MODIS NDVI, while Landsat 8 OLI image was obtained from collection 1 with the status of Level-2, which was already corrected geometrically and radiometrically.However, another pre-processing still needs to be done for Landsat 8 OLI image.Therefore, we provide data preprocessing code compilation for Landsat 8 OLI and MODIS NDVI.Pre-proceeding applied to Landsat 8 OLI mainly to re-project the image to adapt to the local coordinate system (WGS 84 UTM zone 49S) and to crop/subset the image to match the area of interest boundaries.In addition, NDVI was derived using near-infrared (Band 5) and red (Band 4) channels of Landsat 8 OLI.All the pre-processing for Landsat can be done using standard remote sensing software using image reprojection and cropping module, or it also can be done in R in case of automation when processing large numbers of data using a combination of GRASS GIS reprojection and cropping module and iteration processing code.
A bit more complex pre-processing needs to be passed by MODIS NDVI image.Firstly, MODIS NDVI image needs to be read from the original file format (*.hdf), then mosaicking, re-projecting, subset/crop, and exporting into a standard image format (*.gtif) that were done by employed MRT, which integrated to R code.Secondly, resulted MODIS NDVI from the first step needs to be resampled and aligned so that the pixel size and position match the pixels of the Landsat NDVI pair.This process was done using GRASS-GIS modules that integrated into R. Thirdly, the resulting MODIS NDVI in the second step is calibrated radiometrically using a linear model [11] based on the correlation between pixels of MODIS (l0) and Landsat (h0) NDVI that used as a reference image (MODIS NDVI at the composite date of 29-August-2013 and Landsat 8 OLI acquired on 12-September-2013).Fourthly, the resulting MODIS NDVI is then exported into ENVI image format (*.hdr).

Image blending
The blending of MODIS NDVI and Landsat NDVI was done using the STARFM algorithm, which was adapted into R code.The adaptation of the STARFM algorithm into R code mainly to set up the process parameters such as input data of base pair (c0 and f0) and predictor images (c1), the output image (f1), and some detailed parameters of the algorithm.In addition, R code also can help to run batch processing that involves a series of MODIS images to provide a series of high temporal resolution images at Landsat spatial resolution.

Evaluation of the synthetic image
To define whether the result of generated synthetic Landsat 8 OLI NDVI using MODIS NDVI is reasonable, the Root Mean Square Error (RMSE) and Pearson's Correlation (r) were calculated between the synthetic image against the actual one.The synthetic Landsat 8 OLI NDVI was generated using MODIS NDVI (c1) with a composite date of 30 September 2013, while the actual one (h1) was acquired on 14 October 2013.Small areas that were covered by cloud were excluded from the evaluation.

Results and Discussion
This study integrated several methods to generate synthetic Landsat NDVI, covering image preprocessing and blending steps.The integrated method produced the same result as pre-processed image and the synthetic image compared to one from a separated or conventional procedure.However, the integrated method seems faster than the separated or conventional procedure in generating synthetic Landsat NDVI.This acceleration is due to both processes of image pre-processing and blending steps that can be run right away sequentially.In addition, R programming also has a module/package that can run parallel computation that allows exploiting multicore processor ability for faster computational processes.This may cut the computational time spent, particularly when processing a series of images.The complete code for the entire process is available upon request to the author.
STARFM in this study produced a reasonable result of synthetic Landsat 8 OLI NDVI compared to the actual one (Figure 2), indicated by the RMSE value of 0.13 and r value of 0.73.This result seems comparable with the previous study of STARFM performance evaluation [9,5,13].Even though the statistical indices of RMSE and r indicate the STARFM has resulted in a reasonable result of synthetic Landsat 8 OLI NDVI for the case of this study, Figure 2 shows that some details of the feature shape in the actual image are not yet successful in being captured or generated in the synthetic image.For example, the edges of the object in the lower right section of the image, which seems to be paddy field parcels, are less clear to be recognized in the synthetic NDVI (Figure 2-B) than the original one (Figure 2-A).Further, in the center and upper part of the image, which seems to be urban areas, the road corridor and settlement shape are also less clear to be recognized (Figure 2-B) than the actual one (Figure 2-A).This probably occurs because of the time discrepancy between the composite date of the MODIS NDVI as a predictor (30 September 2013) and the acquired date of the original Landsat 8 OLI NDVI (14 October 2013) is quite long (15 days), even though it is still in the range of MODIS NDVI 16 days compositing time.This is in accordance with the previous study that found that time discrepancies may affect the quality of synthetic Landsat images [14,15].However, the synthetic Landsat 8 OLI NDVI, in general, can capture the spatial variance of the vegetation greenness in the study area.This is most probably because the synthetic Landsat NDVI was created directly from both Landsat and MODIS NDVI input data instead of the original reflectance bands and then calculated NDVI, as suggested by the previous study of Jarihani et al. [16].
Another programming application, as well as another blending algorithm, are available widely from many sources; however, this study may be one of the simple, useful examples which can be applied in various operating system platforms following the FOSS GIS and R programming compatibility.

Conclusion
A high spatial and temporal remote sensing image is necessary for the study of earth surface dynamics at the local scale; however, limitations due to cloud cover occurrence may restrict the number.Fortunately, image blending algorithms such as STARFM can fill the gap; however, it is not straightforward because several pre-processing is required, and the algorithm itself needs to be customized for practical use.Simple programming languages such as R may help integrate all the necessary processes, including FOSS GIS inside.This study shows that integration in R is quite successful in producing the synthetic Landsat 8 OLI NDVI as an example; however, the quality of the synthetic image depends on the performance of the blending algorithm used.In this study, the performance of STARFM is quite reasonable according to statistical indices of RMSE and r.This study may help researchers who are seeking a continuous NDVI at Landsat spatial resolution for their study.

Figure 1 .
Figure 1.The study area (indicated using the red box) is at the eastern slope of Mount Merapi in Central Java, Indonesia (base map source from ESRI).

Figure 2 .
Figure 2. Comparison of the original Landsat 8 OLI NDVI captured on 14-October-2013 (A) and the synthetic one generated using STARFM (B).
8 OLI 30m spatial resolution and MODIS NDVI ver.006 16 days composite image 250m spatial resolution (MOD/MYD13Q1) obtained from the NASA's Land Processes Distributed Active Archive Center (LPDAAC) server and the Comprehensive Large Array-Data Stewardship System (CLASS) of the NOAA data server, respectively.For this study, the reference image consists of Landsat 8 OLI NDVI, acquired on 12 September 2013, with the companion pair of MODIS NDVI composite date of 29 August 2013.At the same time, the MODIS NDVI used as a predictor (