Phenological Analysis of Sub-Alpine Forest on Jeju Island, South Korea, Using Data Fusion of Landsat and MODIS Products

: Climate change poses a disproportionate risk to alpine ecosystems. Effective monitoring of forest phenological responses to climate change is critical for predicting and managing threats to alpine populations. Remote sensing can be used to monitor forest communities in dynamic landscapes for responses to climate change at the species level. Spatiotemporal fusion technology using remote sensing images is an effective way of detecting gradual phenological changes over time and seasonal responses to climate change. The spatial and temporal adaptive reﬂectance fusion model (STARFM) is a widely used data fusion algorithm for Landsat and MODIS imagery. This study aims to identify forest phenological characteristics and changes at the species–community level by fusing spatiotemporal data from Landsat and MODIS imagery. We fused 18 images from March to November for 2000, 2010, and 2019. (The resulting STARFM-fused images exhibited accuracies of RMSE = 0.0402 and R 2 = 0.795. We found that the normalized difference vegetation index (NDVI) value increased with time, which suggests that increasing temperature due to climate change has affected the start of the growth season in the study region. From this study, we found that increasing temperature affects the phenology of these regions, and forest management strategies like monitoring phenology using remote sensing technique should evaluate the effects of climate change.


Introduction
Rising global temperatures and atmospheric carbon dioxide levels, changes in precipitation frequency, and increasing severity of extreme climatic events are some of the impacts of climate change. The resultant effects on the ecosystems due to these changes will intensify in a negative direction. These effects include changes in the distribution of forest organisms and populations as well as ecosystem function and composition [1][2][3][4]. In particular, forest biodiversity adjusting to new environmental conditions will lead to species migration. Due to their vertical dimensions, mountains have unique climatic and biogeographical features and create gradients of temperature, precipitation, and insolation. Species compositions in forests are therefore highly susceptible to climate change, and as a result vulnerable species and populations may become extinct.
Thus, climate change poses a disproportionate risk to alpine ecosystems because they are determined by relatively strict climatic parameters [5][6][7][8][9]. Forests and trees play a critical role in supporting basic human life by performing many critical ecosystem services, such as soil stabilization and erosion control, air quality control, climate regulation, carbon sequestration, biodiversity, and recreation and tourism [10]. Therefore, changes in climatic parameters have a substantial impact on both the physical environment and biosphere in a sub-alpine ecosystem.
Effective monitoring of forest phenological responses to climate change in forested regions is critical to predicting and managing threats to populations [11]. However, identifying community-level phenological responses of forests through field surveys is a large-scale and time-consuming task [12][13][14]. Because ecosystems exhibit spatiotemporal variation, phenological characteristics should be analyzed for spatiotemporal variability in forest type at the species level [15,16].
Remote sensing can be used in the forest of a dynamic landscape to monitor phenological responses to climate change at the species level [17][18][19]. Several studies have conducted forest phenological monitoring, with a focus on the type of land cover, using multi-temporal data (LandSat ETM+) [20,21]. These studies were limited to daily monitoring due to the poor temporal resolution of the data acquired by these sensors (Quickbird, Ikonos, Landsat, and SPOT5) [22,23]. An example of this study examined forests in the community and their species [24].
Spatiotemporal data fusion is an effective way to solve the trade-off problem between spatial and temporal resolution. Data fusion is used to create synthetic reflectance data with high spatial and temporal resolutions [25]. This method combines the spatial information from images with a high spatial resolution with the temporal information from satellite images with a coarse spatial resolution but a short revisit period to generate images with high spatial and temporal resolution. This method is particularly useful for detecting phenological changes and seasonal responses to climate change [26,27]. The spatial and temporal adaptive reflectance fusion model (STARFM) [26] is perhaps the most widely used data fusion algorithm for combining Landsat and MODIS (Medium Resolution Imaging Spectroradiometer) imagery [28,29], showing high estimation accuracy at a local scale [30,31]. Moreover, it is one of the few data fusion methods that result in synthetic Landsat-like surface reflectance [32]. Previous studies have been focused on large-scale [33,34] and single time-series [35,36] analysis to identify phenology changes at a forest vegetation community or species because of the time-consume and data processing. Thus, the purpose of this research is to use spatiotemporally fused Landsat and MODIS imagery to identify forest phenological characteristics and change at the species-community level. In this study, the normalized difference vegetation index (NDVI), most widely used for phenological analysis [37,38], was used for analyzing the changes of phenology in the sub-alpine ecosystem. The employment of data fusion using Google data processing could help overcome the shortcomings of the previous approaches.

Study Area
The study site (Figure 1) is Mount Halla on Jeju Island (33 • 13 -33 • 36 N, 126 • 12 -126 • 57 E), Republic of Korea. This site was chosen because it exhibits a high biodiversity and latitudinal gradients in species diversity throughout its elevation (300-1950 m) and contains many endemic and endangered species and habitats, typical of a densely populated area. In 2019 (1 January-31 December), the average temperature of Jeju Island was 17.1 • C which is 0.9 • C higher than the average temperature (16.2 • C), the second highest (after a record high of 17.3 • C in 1998) since 1961 [39]. These temperatures led to changes in the composition of habitat, species, forest type, and flora in the sub-alpine area of Mount Halla, South Korea ( Figure 1). The lowest altitude within this study site is 1218 m, the highest is 1945 m, and the average is 1544 m. Based on data from the past 20 years, this area has an annual average temperature of 9.4 • C, and an annual average cumulative monthly precipitation of 288 mm. Mount Halla, South Korea ( Figure 1). The lowest altitude within this study site is 1218 m, the highest is 1945 m, and the average is 1544 m. Based on data from the past 20 years, this area has an annual average temperature of 9.4 °C, and an annual average cumulative monthly precipitation of 288 mm. The conservation area on Jeju Island is characterized by high biodiversity and forms the primary ecological axis of Jeju Island. Within the study area, deciduous broad-leaved trees and Abies koreana are the dominant tree species and formed the target species for this study; they occupy 20.7% and 18.3% of the study area, respectively. However, the forest ecosystems have undergone intensive disturbances from climate change and reclamation for urban development. Even the recent climate change has affected forest types and flora. In particular, Abies koreana is a rare endemic species on Jeju Island that is at threat from and being killed by other species as well as recent climate change [40][41][42].

Research Workflow
In this study, the identification of forest phenological characteristics and creation of a time-series of phenological changes at the species-community level was conducted in three main stages ( Figure 2): (1) medium-spatial and high-temporal resolution products (Landsat OLI, Thematic Mapper (TM), MODIS MOD09GQ, and MOD13Q1) were obtained from the Google Earth Engine (GEE, https://developers.google.com/earth-engine/datasets, accessed on 2 March 2021); (2) the Landsat 30-m images were fused with the MODIS 500-m images using STARFM, to produce synthetic imagery, and an accuracy assessment through the R-square and RMSE was conducted to set the optimal parameters for STARFM; (3) NDVI images, predicted through STARFM, were then classified by referencing the forest community distribution map of the current Mount Halla Absolute Conservation Area, and the time-series was analyzed for phenological changes in growth characteristics. The conservation area on Jeju Island is characterized by high biodiversity and forms the primary ecological axis of Jeju Island. Within the study area, deciduous broad-leaved trees and Abies koreana are the dominant tree species and formed the target species for this study; they occupy 20.7% and 18.3% of the study area, respectively. However, the forest ecosystems have undergone intensive disturbances from climate change and reclamation for urban development. Even the recent climate change has affected forest types and flora. In particular, Abies koreana is a rare endemic species on Jeju Island that is at threat from and being killed by other species as well as recent climate change [40][41][42].

Research Workflow
In this study, the identification of forest phenological characteristics and creation of a time-series of phenological changes at the species-community level was conducted in three main stages ( Figure 2): (1) medium-spatial and high-temporal resolution products (Landsat OLI, Thematic Mapper (TM), MODIS MOD09GQ, and MOD13Q1) were obtained from the Google Earth Engine (GEE, https://developers.google.com/earth-engine/datasets, accessed on 2 March 2021); (2) the Landsat 30-m images were fused with the MODIS 500-m images using STARFM, to produce synthetic imagery, and an accuracy assessment through the R-square and RMSE was conducted to set the optimal parameters for STARFM; (3) NDVI images, predicted through STARFM, were then classified by referencing the forest community distribution map of the current Mount Halla Absolute Conservation Area, and the time-series was analyzed for phenological changes in growth characteristics.

Data Collection
The method used in this study relies on high-temporal and multispectral resolution images for creating a time-series of phenological changes at the species-community level. For this, remote sensing data were obtained from the GEE platform, as it would be too

Data Collection
The method used in this study relies on high-temporal and multispectral resolution images for creating a time-series of phenological changes at the species-community level. For this, remote sensing data were obtained from the GEE platform, as it would be too timeconsuming to collect sufficient data through field surveys. GEE is a cloud-based platform that enables the online processing of satellite exploration data [43]. Fast processing and results can be obtained by linking Google supercomputers with Java or Python code. Additionally, GEE reduces the data processing time and the need for advanced computers as the processing of existing satellite image data is performed online through coding.
For this study, Landsat OLI and TM products, as well as MODIS MOD09GQ and MOD13Q1 products, were obtained from GEE. Landsat 30-m resolution observations provide sufficient spatial detail for surface and change monitoring [44]. However, cloud cover and the 16-day revisit cycle of this product make it less useful for studying rapidly evolving global biophysical processes during the growing season [45]. Therefore, data fusion was attempted between Landsat and MODIS products. MODIS products are the most commonly used daily satellite imagery for observing the seasonal phenology of vegetation [46].
The Landsat TM sensor has seven bands with a spatial resolution of 30 m. MOD09GQ is a MODIS product that has undergone all atmospheric corrections and provides two surface spectrum reflectances (bands 1 and 2) with a spatial resolution of 250 m. MOD13Q1, a vegetation indices product, provides two bands: the NDVI and enhanced vegetation index, each with a spatial resolution of 250 m. NDVI is a vegetation index used to estimate vegetation vitality, productivity, and green cover rate [47]. It has been widely used for phenological detection of terrestrial ecosystems.
We obtained Landsat images with less than 10% cloud cover to minimize the influence of clouds [48,49] during 2000, 2010, and 2019, a period covering the entire growing period, and paired them with corresponding MODIS data, as shown in Table 1. Three NDVI products were constructed using bands 5 and 4 of the Landsat OLI products, bands 4 and 3 of Landsat TM products, and bands 2 and 1 of MOD09GQ. Using NDVI instead of bands allows for high-accuracy phenological analysis [31]. Both the Landsat and MODIS satellite survey data were geometrically corrected with the same coordinates (EPSG Geodetic Parameter Dataset: 32,652). Due to the characteristics of MODIS, the constructed NDVI was multiplied by 10,000, and the result was derived after downscaling to a 30-m likelihood resolution. The entire processing process was carried out on the GEE platform.

Auxiliary Data
A map of forest type in the study area was acquired from the Korean Ministry of Environment. The map depicts vegetation classes (Abies koreana derived from multi-season 2015-2017 aerial imagery mapped at 20 cm resolution). The phenology was then analyzed based on the forest type map. Also, the average monthly temperatures of 2000, 2010, and 2019 were obtained from the automatic weather stations (AWS) provided by the Korean Meteorological Administration (KMA) to find the relationship between NDVI and temperature increase.

Data Fusion
Using STARFM, we fused the Landsat 30 m data with the MODIS 500 m data to generate time-series data with a high spatial resolution [31]. STARFM was chosen as the data fusion algorithm because it is the basis of other fusions and has few restrictions [38].
To increase the fusing accuracy, we minimized the influence of clouds by using 16day composites of the MODIS data covering the prediction period [48]. Images for the same period were fused and predicted using only clean Landsat and MODIS images without clouds. For the 3 years of 2000, 2010, and 2019, a total of 54 images were fused via this method.
Using STARFM to estimate the reflectance over the fusion period assumes that the time-series change in MODIS data is the same as that in the Landsat data. Weighted combining was performed, considering the temporal, spectral, and spatial similarities of MODIS and Landsat data per pixel. We used STARFM to directly blend vegetation indices (VIs) derived from SPOT5 and MOD09Q. Several studies have found that the STARFM method performs better when directly fusing VIs than when the reflectance is fused before the VIs are calculated [50]. The value of the prediction timing for each pixel was calculated using the following equation: where w is the size of the search window and is a parameter that determines the spatial similarity at the predicted position; W is a weight, determined according to the temporal, spectral, and spatial similarity between MODIS and Landsat data. Time similarity at the fused time is calculated using the time difference between the MODIS data and Landsat data. Spatial similarity is given a higher weight as it is closer in terms of the distance between the predicted position (the center of the search window) and the surrounding pixels. Spectral similarity is calculated by assigning a high weight to the pixels having similar spectral difference between surface reflectance of MODIS-Landsat. In this study, we used the Python code of STARFM from [51].

Parameters of STARFM
The image fusion accuracy and the time required for STARFM depend on two parameter values: the size of the search window and the number of classes classifying spectral differences between the surface reflectance of MODIS-Landsat. To improve the accuracy, it is important to determine the ideal size of the search window and the number of classes [29]. The size of the search window, which sets the areas searching for similar pixels, was set to an appropriate value through pilot prediction because, otherwise, the result may be less accurate than a certain value, depending on the target location [25]. In most studies, the size of the search window has a more direct influence on the prediction accuracy rather than the number of classes; therefore, the size of the search window will be set first, followed by the number of classes [27]. For the pilot prediction, the size of the search window and the number of classes were set by referring to the results of previous studies, and the optimal value such as window size and the number of classes was chosen in consideration of the fusion accuracy and the time required.
Two statistical criteria, the coefficient of determination (R 2 ) and the root mean square error (RMSE), were selected to set the optimizing value of window size and the number of classes. These criteria are widely used to evaluate phenology extraction results [29,52]. They are calculated as follows: where M i is the measured value, F i is the linear fitting value, A i is the average value of the measured data, and E i is the estimated value.

Phenological Analysis Based on Forest Types
NDVI is very effective for monitoring vegetation changes over large areas [53], and is easily calculated from satellite imagery with the following equation: Here, NIR refers to the value of the near-infrared band, and Red refers to the value of the red-light band. NDVI ranges from −1 to 1, and values closer to 1 generally indicate higher vegetation vitality. The seasonality of local plants can be analyzed through changes in NDVI values throughout the year. NDVI image data predicted through STARFM were classified according to forest type to analyze changes in growth characteristics. The image data were classified into forest types using the ArcGIS program and based on the community distribution map of the current Mount Halla Absolute Conservation Area.

Results of STARFM
The two parameter values, namely the size of the search window and the number of classes, were set through pilot prediction. In accordance with previous studies [27], the pilot prediction was conducted by setting the size of the search window to 5, 10, 25, 50, 100, 150, 200, 250, and 300 [27]. As a result of the pilot prediction, when the window size was set to 50, the R 2 (0.795) and RMSE (0.0402) values indicated a sufficient accuracy of the model (Figure 3a). The number of classes was set to 10, 20, 40, and 80, and the pilot prediction was carried out. The RMSEs for the different number of classes were similar, but the highest accuracy (R 2 = 0.795) was obtained when the number of classes was set to 10 ( Figure 3b). These STARFM accuracy evaluation results were similar to those (R 2 : 0.69-0.86, RMSE: 0.06-0.11) of [38] and were acceptable results. Based on the pilot prediction results, image fusion and prediction were performed with a search window of 50 and a number of classes of 10. Even though the images from the prediction period were used as 16-day composites, there was one case where a cloud existed in the prediction result. Because clouds were found mainly in the images from the summer period, the images of the summer period were excluded from later seasonality analysis.

Phenological Analysis
After applying STARFM to the sub-alpine area of Mount Halla, the predicted images were classified using the ArcGIS program and a community distribution map of the absolute conservation area of Mount Halla. The classification was made for eight communities, and a total of 54 × 8 = 432 fused imagery were created. Out of the eight communities, the Abies koreana and the deciduous broad-leaved forest were compared because they are the dominant types in this study area.
The results ( Figure 5) showed that the time-series phenology pattern of Abies koreana had opposite characteristics to the deciduous broad-leaved forest. In the case of deciduous broad-leaved forests, the NDVI value showed a tendency to increase gradually over each season (Figure 5b). Similar to deciduous broad-leaved forests, the NDVI of Abies koreana showed a tendency to increase with each season (Figure 5a). In particular, in March and April, the average value of the total NDVI showed that the number of regions with high NDVI values gradually increased, as indicated by the average total NDVI. These results could serve as evidence to support the results of the previous study [52] that spring has higher productivity effects than fall.
Deciduous broad-leaved forests generally have higher NDVI values than Abies koreana since the latter is a coniferous tree. Abies koreana showed the highest NDVI value in September due to its growth characteristics. In contrast, in 2000, the NDVI value of Abies koreana was low, and the start of the season was later than that of deciduous broad-leaved forests. Considering that Abies koreana is a coniferous tree, such an increase in NDVI is a peculiar phenomenon. These results are similar to those of previous studies [17,54] and are likely due to climate change creating a longer growing season. Therefore, for rare species like Abieas koreana, sensitive to climate change, continuous NDVI monitoring is needed to observe climate change phenomena and preserve ecosystems.
September due to its growth characteristics. In contrast, in 2000, the NDVI value of Abies koreana was low, and the start of the season was later than that of deciduous broad-leaved forests. Considering that Abies koreana is a coniferous tree, such an increase in NDVI is a peculiar phenomenon. These results are similar to those of previous studies [17,54] and are likely due to climate change creating a longer growing season. Therefore, for rare species like Abieas koreana, sensitive to climate change, continuous NDVI monitoring is needed to observe climate change phenomena and preserve ecosystems.

Phenology and Temperature
To examine the relationship between NDVI and temperature increase due to climate change, the average monthly temperatures in 2000, 2010, and 2019 were obtained from the AWS provided by the KMA (Figure 6). Considering the observation points in five different AWS locations around the study area, temperature data were obtained from the observatory being able to include the temporal range of this study (Figure 6a).
The average monthly temperature increases slightly over time, except in May, where the NDVI value is relatively high. In particular, from the monthly average temperature

Phenology and Temperature
To examine the relationship between NDVI and temperature increase due to climate change, the average monthly temperatures in 2000, 2010, and 2019 were obtained from the AWS provided by the KMA (Figure 6). Considering the observation points in five different AWS locations around the study area, temperature data were obtained from the observatory being able to include the temporal range of this study (Figure 6a). change, the increase in average temperature occurred at the start of the season rather than at the end of the season (Figure 6b). These changes could be used as evidence to show the rise in NDVI during March and April ( Figure 5). Previous studies have shown similar results [17,54], and it can be speculated that a temperature increase due to climate change can affect sensitive forest communities. In contrast, the target area also receives snow falls as much as the alpine region from late November to February. As in previous studies [55,56], phenomena and sub-alpine ecosystems are affected by the amount of snow cover and the period of cover. Furthermore, a temperature rise due to climate change affects the duration and accumulation of snow and must be observed carefully. Through analyzing phenology by forest type in the target region ( Figures 5 and 6), the seasonal growing characteristics of species sensitive to climate change were confirmed, as in previous studies. In particular, the use of the STARFM technique in target sites with rainy climates sufficiently complements the qualitative limitations of Landsat data, which are heavily influenced by the weather, to increase the analytical possibilities. In addition, surface-based in situ phenological analysis provides detailed leaf development information, however, it has limitations of data availability spatially and temporally [57]. To monitor the seasonal dynamics of vegetation greenness, it is necessary to collect phenological data using remote sensing.

Conclusions
This study uses fusion technology to build on previous studies that have focused on time-series phenological change at a species-community level to identify phenological changes in the sub-alpine area. Our study constructed NDVI on the GEE platform using Landsat OLI and TM products with MODIS MOD09GQ and MOD13Q1 products and variables such as the coordinate system and window size. These approaches significantly re- The average monthly temperature increases slightly over time, except in May, where the NDVI value is relatively high. In particular, from the monthly average temperature change, the increase in average temperature occurred at the start of the season rather than at the end of the season (Figure 6b). These changes could be used as evidence to show the rise in NDVI during March and April ( Figure 5). Previous studies have shown similar results [17,54], and it can be speculated that a temperature increase due to climate change can affect sensitive forest communities. In contrast, the target area also receives snow falls as much as the alpine region from late November to February. As in previous studies [55,56], phenomena and sub-alpine ecosystems are affected by the amount of snow cover and the period of cover. Furthermore, a temperature rise due to climate change affects the duration and accumulation of snow and must be observed carefully.
Through analyzing phenology by forest type in the target region ( Figures 5 and 6), the seasonal growing characteristics of species sensitive to climate change were confirmed, as in previous studies. In particular, the use of the STARFM technique in target sites with rainy climates sufficiently complements the qualitative limitations of Landsat data, which are heavily influenced by the weather, to increase the analytical possibilities. In addition, surface-based in situ phenological analysis provides detailed leaf development information, however, it has limitations of data availability spatially and temporally [57]. To monitor the seasonal dynamics of vegetation greenness, it is necessary to collect phenological data using remote sensing.

Conclusions
This study uses fusion technology to build on previous studies that have focused on time-series phenological change at a species-community level to identify phenological changes in the sub-alpine area. Our study constructed NDVI on the GEE platform using Landsat OLI and TM products with MODIS MOD09GQ and MOD13Q1 products and variables such as the coordinate system and window size. These approaches significantly reduced the time required to collect large amounts of satellite data and the uncertainty of the preprocessing method for each step. In fact, this process took only a day which included the work of coding JavaScript in the GEE platform. This process would have taken over a month if Landsat and MODIS data were downloaded from the United States Geological Survey (USGS) EarthExplorer (https://earthexplorer.usgs.gov/, accessed on 2 March 2021) to calculate NDVI. We used STARFM, which is the most widely used spatiotemporal data fusion model and shows high prediction accuracy at local scale, to fuse Landsat images for dates not covered by data, using fusion of Landsat and MODIS images. For the fusion, Landsat and MODIS data from the same period, in the form of clean images without clouds, and the forecast period data were preprocessed by reducing the effect of clouds (as much as possible) through 16-day composites. Eighteen images were predicted for each of the three years (2000, 2010, and 2019, from March to November).
Through phenological analysis in 2000, 2010, and 2019 using estimated NDVI values obtained from STARFM, we obtained three major results. First, by examining the phenology of forest types, we confirmed that NDVI values increase with time. In particular, in March and April, the average total NDVI showed that the number of regions with high NDVI values gradually increased, which indirectly confirmed that the start of season was shifting earlier. It also suggests that increasing temperature due to climate change have affected the start of season in the study regions. Second, we compared the average monthly temperatures of 2000, 2010, and 2019 to find a link to the temperature rise caused by climate change. We found that the increase in average temperature generally occurred at the start of season, rather than the end. Lastly, the STARFM technique is particularly useful for target sites with rainy climates, where the applications of Landsat data are limited by the weather. Although the climatic characteristics of the site were factored into the STARFM process, and the resulting accuracy (RMSE: 0.0402, R 2 : 0.795) was sufficient, it was not as high as that of other studies [58]. The phenology analysis using remote sensing techniques like STARFM is more useful to monitor the seasonal dynamics of vegetation greenness than surface-based situ phenological analysis in the perspective of data availability spatially and temporally. Additionally, the more clean, cloudless data were used, the better the results were.
In the future, this method can be supplemented by improving the image fusion process using a high-performance machine capable of processing more images. In addition, this study's phenological analysis was only performed for specific years. It seems that it is necessary to determine the trends over multiple years to improve the accuracy and reliability of the method. Furthermore, the average value before and after the relevant year should be used, rather than just selecting a specific year. By improving the disadvantages of the mentioned data period and increasing the level of mechanical defects capable of processing large images, a phenology analysis of a wider area could be attempted as the next step in the study. Additionally, as forests respond to climate change, it would help manage the ecosystem phenology of these regions using this remote sensing technique.