Chi Hong Lim · Song Hie Jung · Nam Shin Kim · Chang Seok Lee

Abstract Phenology is a valuable attribute of vegetation to assess the biological impacts from climate change. A challenge of phenological research is to obtain information on both high temporal resolution and fine spatial scale observations.Here,we constructed an air temperature map based on temporal merging and spatial interpolation algorithms to overcome the cloud-related problem from the MODIS LST product. Then, we derived the accumulated growing degree days (AGDD) from the constructed mean air temperature map to use as a meteorological indicator.Further, we verified the indicator with the seasonal mean air temperature and the green-up date of a Quercus mongolica forest determined from the field-based measurements. The AGDD threshold for each Q. mongolica forest when the first leaf has unfolded was detected from the EXG trajectory extracted from digital camera images. A comparison between meteorological and MODIS-derived AGDD showed good agreement between them. There was also high consistency between DoYs extracted from AGDD and EVI based on curvature K for Q. mongolica forests of 30 sampling sites throughout South Korea. The results prove that microclimatic factors such as elevation,waterbody, and land-use intensity were faithfully reflected in the reconstructed images. Therefore, the results of this study could be applied effectively in areas where microclimatic variation is very severe and for monitoring phenology of undergrowth, which is difficult to detect from reflectance imaging.

Keywords Climate change · Digital camera · Growing degree days · MODIS · Phenology · Quercus mongolica

Introduction

Vegetation phenology, referred to as the timing of recurring biological cycles, has been considered a sensitive indicator of climate change at local and global scales(Zhang et al.2004;Cleland et al.2007;Hassan et al.2007;De Beurs and Henebry 2004; Inouye 2008; Richardson et al. 2009a, 2018; Jochner et al. 2013; Fu et al. 2014;Allstadt et al. 2015; Fitchett et al. 2015) because establishment, growth, reproduction, and survival of a plant are directly related to the environmental conditions (Sykes 2009). Among the numerous techniques to observe how phenology has shifted in recent decades are ground-based observations (Reich and Borchert 1984; Walther et al.2002; Menzel et al. 2006; Primack and Higuchi 2007),digital repeat photography (Richardson et al. 2009a; Ide and Oguma 2010; Hufkens et al. 2012; Klosterman et al.2014), satellite remote sensing of ecosystem production(Schwartz and Reiter 2000; Fisher and Mustard 2007;Jeong et al. 2013; Allstadt et al. 2015; Zipper et al. 2016;Richardson et al. 2018), and monitoring atmospheric carbon dioxide (CO2) concentrations (Baldocchi 2003;Richardson et al. 2009b). Satellite remote sensing has the advantage of providing multi-decadal records of vegetation phenology across larger spatial scales than other techniques(Fitchett et al. 2015). In particular, recent advances in temporal and spatial resolution and ease of use have made remotely sensed data derived from MODIS (MODerateresolution imaging spectroradiometry)the most commonly used methods to gather phenological information at regional and global scales (Zhang et al. 2003, 2013).

The data derived from MODIS for the study of vegetation phenology can be divided into two groups: (1) a vegetation index extracted by measuring leaf reflectance and (2) a meteorological index extracted from remotely sensed land surface temperature (LST). Direct observation of changes in leaf traits such as color,area,and vitality are,of course,the most efficient means for monitoring,but has a critical limitation,which is related to the effects of cloud cover and cloud shadow. Cloud cover is the primary disturbance of leaf reflectance measurements, and thus,degrades the reliability of the vegetation index(Park 2013).The LST is similarly affected by cloud contamination(Neteler 2010;Zhang et al.2013).To resolve this problem,researchers often remove the contaminated pixel using quality assessment information in the MODIS product(Kang et al.2003).Then,they spatially interpolate the void map areas using a geographic interpolation system (GIS).However, interpolation of the vegetation index would not valid because the forest is composed of various patches of plant communities that each have intrinsic characteristics.In contrast, the interpolation of air temperature is more accurate than an interpolated vegetation index because the atmospheric condition is fluid. Therefore, the use of a meteorological indicator as an assessment tool in phenology is worth consideration for complex terrain,such as the eastern mountainous area of South Korea.

The relationship between temperature and phenophases such as flowering and leaf unfolding is well known and has been widely reviewed (Cleland et al. 2007). In particular,the concept of a heat unit, measured in growing degree days (GDD, °C d), has vastly improved description and prediction of phenological events compared with other approaches such as day of year or number of days(McMaster and Wilhelm 1997). Although AGDD has hardly been considered as a link between changes in climate and phenology in ecology(Cayton et al.2015),it has a long history of use in predicting plant and insect phenology in agriculture (Parry and Carter 1985; Bonhomme 2000). Moreover, temperature is known as a major driver of phenology (Miller-Rushing et al. 2010).

The main objectives of this study include (1) the reconstruction of an air temperature map based on the algorithm to overcome the cloud-related problem from the MODIS LST product and(2)evaluation of the applicability of a meteorological indicator as a tool to monitor phenology. To achieve these aims, first, we reconstructed a mean air temperature map at 8-d intervals using the MODIS Terra/Aqua LST product based on temporal merging and spatial interpolation algorithms. Then, the accumulated growing degree days (AGDD) as a meteorological indicator were deducted from the reconstructed mean air temperature map. The indicator was verified with the seasonal mean air temperature and the green-up date of Quercus mongolica forest determined from field measurements.

Materials and methods

Study area

The study region, South Ko rea (33°–38°N, 125°–131°E)(Fig. 1), covers 100,033 km2and has a population of 51.0 million people (Statistics Korea 2016). The annual mean air temperature is between 10 and 15 °C, and annual precipitation is between 1000 and 1900 mm (Korea Meteorological Administration 2016). The landform is mountainous,and most of the land is not arable.Lowlands,located primarily in the west and southeast,constitute only 30% of the total land area. The total forest area in South Korea covers 63,346 km2, and most forests are deciduous broadleaved and evergreen needle-leaved (Korea Forest Service 2015).The temperate deciduous broadleaved forest zone covers most of the Korean Peninsula, except for the southern and western coastal areas and the high mountainous area. It is usually dominated by various broadleaved deciduous forests, but coniferous plantation forests also appear frequently in lowland areas, which were intensively used. In this study, Mongolian oak (Q. mongolica), a dominant species in the temperate deciduous broadleaved forest zone on the Korean Peninsula, was selected as the target species.

Data acquisition and preprocessing

The details for the available data sources with remotely sensed,station-based,and other auxiliary geographical data are listed in Table 1. The two types of MODIS data,including LST/emissivity products (MOD11A2.006,MYD11A2.006) and a surface reflectance product(MOD09A1.005) from 1 January 2015 to 26 June 2015were downloaded. In the case of LST/emissivity products,the overpass times for the day image were 10:30(MOD11A2) and 13:30 (MYD11A2) local solar time. The times for night images were 22:30 (MOD11A2) and 01:30(MYD11A2).

Fig. 1 Map of the study area and spatial distribution of the primary sampling points (sites where digital cameras were installed, meteorological stations, and verificationation sites in Q. mongolica forests).Meteo-station: meteorological station

Table 1 Summary of available data in the study

The MODIS re-projection tool software (MRT V4.1)was used to re-project the original Sinusoidal (SIN) to a TIFF file using the Universal Transverse Mercator projection (UTM Zone 52N, WGS82 ellipsoid) (Neteler 2010).For calculating a raster at the pixel level,all MODIS images were resampled to 1000 m with the same spatial resolution of the LST data using nearest neighbor resampling (Zhang et al. 2013).

MODIS 8-d mean LST images (MOD11A2/MYD11A2) were used in this study for two reasons: (1)They can minimize the loss of data due to clouds to a large extent, and (2) there is a high consistency between the estimations of heat accumulation from meteorological daily air temperature and 8-d composited mean air temperature (Zhang et al. 2013).

The enhanced vegetation index(EVI;Huete et al.2002)was calculated from the MODIS surface reflectance product based on Eq. 1:

where ρNIR,ρRED,and ρBLUEare near-infrared, red, and blue bands in MOD09A1, respectively, L is the canopy background adjustment, C1and C2are coefficients, and G is the gain factor. Then, a bit-pattern analysis was implemented to remove cloudy pixels identified by the information in the State flags (16-bit unsigned integer).

Daily air temperature data (minimum and maximum temperature) for 1 January through 2 July 2015 from meteorological stations were acquired from the Korea Meteorological Administration (KMA). Corresponding to the temporal resolutions of MODIS LST products, daily data were averaged at 8-d intervals (Zhang et al. 2013).From daily air temperature data,the AGDD at each station was calculated based on Eq. 5.

Additional auxiliary geographical data(Elevation above sea-level,Distance from sealine,Coordinate X and Y)were used to refine the accuracy during the model fitting and spatial interpolation.All geographical data were resampled to 1000 m pixel resolution.

Commercial digital cameras (model Ltl-6210M, Little Acorn Outdoors, Denmark, WI, USA) were installed in nine Q.mongolica forests in South Korea to determine the AGDD threshold at which the Q. mongolica forest starts spring green-up (Table 2). For camera installation and environmental setting such as orientation, field of view,and distance from canopy, we applied internationally accepted protocol (Tierney et al. 2013). Each camera was set to record high-quality JPEG images at 12:30 local time using intervalometers. From the collected daily images,which captured canopy layer of Q. mongolica community during growing period in 2015 and 2016, we extracted the vegetation index (excess green index; ExG) indicating phenological signal (Richardson et al. 2009a; Klosterman et al. 2014; Lim et al. 2018). The ExG was defined as:

where ρRED,ρGREEN,and ρBLUEare values in the red,green, and blue bands in the digital camera image,respectively. Then, the ExG was smoothed to the 80th percentile using an exponentially weighted moving average(Lim et al. 2018). To derive a green-up start date from the phenological time-series data, we used a sigmoid-based equation (Zhang et al. 2003; Fisher and Mustard 2007;Richardson et al. 2009a; Ide and Oguma 2010; Hufkens et al. 2012; Klosterman et al. 2014; Lim et al. 2018).Finally, the date for each site was attained by obtaining local extrema

i

n the rate of change of curvature K (hereafter, simply curvature K) (Zhang et al. 2003; Klosterman et al. 2014; Lim et al. 2018). In addition, we installed meteorological sensors (model HOBO pro v2 U23-001,Onset Co., Bourne, MA, USA) 1 m above the ground at each site to measure the air temperature during the growing season.Each sensor was protected from direct sunlight by a protective case.Measurements were sampled every minute and logged every 30 min.From daily air temperature data,the AGDD at each study site was calculated based on Eq. 5.

Table 2 Summary of environmental factors and DoY of spring green-up in the sites in which digital camera were installed. AGDD indicates the accumulated growing degree days.We assumed the minimum temperature threshold for the growing of in Q. mongolica as 5 °C based on studies of Kira(1945) and Schenker et al.(2014)

Reconstruction of air temperature map

The overall procedure of the algorithm to filter and convert MODIS LST data into several indices are shown in Fig. 2.First, low-quality pixels due to clouds or other processing failures were filtered using QA bitmap layers (8-bit unsigned integer) for each LST image (Neteler 2010;Zhang et al. 2013). The criteria to determine the pixel quality followed Neteler(2010).The calculated LST values were converted from Kelvin to Celsius.

After the filtering and conversion steps,the 8-d synthetic day/night LST temperature maps were merged from Terra and Aqua LST images. We followed the algorithm of temporal merging proposed by Zhang et al. (2013): (1) At the pixel level, if both input maps had values, the values were averaged. (2) If only one map had values, the single value was used.(3)If neither map had values,it was left as a null pixel to be filled by spatial interpolation in the following step. In the study by Zhang et al. (2013), this temporal merging step was done after the air temperature conversion step.However,in the present study,the number of valid pixels for a regression analysis was often insufficient in the winter due to adverse weather. Therefore, we merged the Terra and Aqua images before the regression analysis.

Fig. 2 Flowchart to filter and convert MODIS LST into several indices that affect vegetation phenology

To convert LST values into air temperature values, the LST values, EVI, and elevation were linked with air temperature values from meteorological stations by multiple stepwise regression analysis (Benali et al. 2012; Zhang et al. 2013).

The final step in the reconstruction process was the spatial interpolation for filling residual empty pixels. We converted valid raster pixels to point vectors for each map.Then, a ordinary cokriging algorithm was applied to the vectors with covariates(elevation above sea-level,distance from sealine, coordinate x and y) (Rossiter 2007). The consistency between the measured and interpolated air temperature was evaluated using the root mean squared error (RMSE):

where xm,iand xi,iare measured values and interpolated values, respectively.

Indices calculation

A series of indices that affect vegetation phenology were derived from the completely reconstructed 8-d air temperature maps.First,GDD(°C d)were calculated using the equation from McMaster and Wilhelm (1997):

where Tmax·tand Tmin·tare the maximum and minimum air temperatures at DoYt(day of year at time t;the day of year is the number of days starting with 1 January),respectively,and Tbaseis the temperature below which plant growth is zero. In this study, we set the base temperature as 5 °C(Gordon and Bootsma 1993). In addition, the GDD based on meteorological station data was estimated to validate the GDD derived from MODIS.

We accumulated 8-d-interval GDDs(AGDD)by simple summation when the GDD exceeded the base temperature(De Beurs and Henebry 2004; Zhang et al. 2013):

where GDDtis the 8-d mean GDD at DoYtand i is the time interval coefficient(GDD from MODIS:8;GDD from field data: 1). AGDDtis the accumulated GDD from the beginning of the time period until DoYt+7. Based on the average of AGDD threshold values determined from field measurements when the spring green-up started in Q.mongolica forests, a MODIS-derived AGDD threshold map was produced to assess the timing of green-up of Q.mongolica in all of South Korea.

Results

Air temperature reconstruction and precision validation

A large number of spatial pixels were still missing in the filtered products after the reduction of low-quality pixels by quality assessment information. Therefore, we introduced an algorithm to fill those missing pixels (Neteler 2010; Zhang et al. 2013). The performance of the reconstruction is illustrated in Fig. 3.For the nighttime Terra and Aqua images at DoY 65, 2015, when most of the pixels were contaminated by clouds, a small number of valid pixels could be used in the reconstruction process(Fig. 3a,b). Among the total pixels (98,837), only 20.17% (19,933)in the Terra image and 21.58%(21,332)in the Aqua image were valid. In addition, the spatial arrangement of the pixels was seriously biased.After merging the Terra image and Aqua image into the synthetic image,the proportion of valid pixels increased to 35.19% (34,779), and the bias in the spatial arrangement was moderated (Fig. 3c). Because of the difference in heating mechanisms between the LST temperature and the air temperature(Zeng et al. 2015), we modeled the air temperature as a function of EVI, elevation, and the LST (Zhang et al. 2013) (Fig. 3d). After the conversion step,the remaining invalid pixels were spatially interpolated with the geographic covariables, and then a complete air temperature image was reconstructed(Fig. 3e). The reconstruction process was applied to all filtered images (a total of 96 MODIS 8-d LST images).

The measured and the interpolated air temperatures were highly consistent. The slope of the linear regression gradients and coefficients of determination (R2) were above 0.9 during the whole season except at the minimum temperature in winter (Appendix S1). The highest consistency was in the spring (maximum: 0.99, minimum: 0.99). The conformity between measured and interpolated data likely depends on sample size, which is the number of pixels available in the interpolation process (Fig. 4).

Estimation of the phenological transition date from AGDD indices

Fig. 3 Differences between original filtered minimum land surface temperature maps for Terra (a), Aqua (b), merged minimum land surface temperature map (c), converted minimum air temperature map (d) and reconstructed complete minimum air temperature map(e) for day of years (DoY) 65 in 2015

Fig. 4 Scatter plot of the relationship between root mean squared error (RMSE) and sample size (scaled by natural logarithm) for the spatial interpolation process

The AGDD threshold for each Q. mongolica forest when the first leaf had unfolded was detected in the ExG trajectory extracted from digital camera images (Table 1).The averages of the AGDD in 2015 and 2016 were 157.75(SD = 3.64) and 158.80 (SD = 4.16), respectively. Based on the result, we set up the AGDD threshold as 159 based on the result extracted from field measurements. Further,we assumed that the first day that the value in MODIS daily AGDD images increased above the AGDD threshold is the start date of green-up in Q. mongolica forests.

In 2015, the areas where the AGDD exceeded the AGDD threshold (159 °C d) appeared first in March(Fig. 5b). At this time, the areas were usually distributed around the lowland areas of a southern island and the southern and southeastern coastal areas of South Korea. In addition, the AGDD image shows the local hotspots centered on the urban areas. In April, the areas expanded on a nationwide scale, except for high-altitude mountainous areas (Fig. 5c). Finally, the entire area of South Korea exceeded the AGDD threshold in May (Fig. 5d).

From the MODIS-derived AGDD images, the threshold map was generated by counting pixelwise the number of DoY (hereafter, DoYAGDD) needed to reach the AGDD threshold (Neteler 2010). In the map (Fig. 6), the obvious gradients of DoY in relation to the geographic condition(e.g., altitude, continentality) can be seen. For example, a variation in the AGDD depending on altitude appeared clearly in the case of Jeju Island (southernmost island in South Korea). Aside from the gradient, a number of areas were distinguished,each appearing as a point pattern,were shown on the map. The most obvious example of the area can be seen in the Seoul Capital Area (the metropolitan area of Seoul in northwestern South Korea).

Fig. 5 Spatial distribution of monthly accumulated growing degree days (AGDD) on the South Korea at days of year(DoY) 1–184, 2015

The meteorological and MODIS-derived AGDD from DoY 1 to DoY 184 were compared to display the discrepancy at a daily scale. Figure 7 shows that there was quite a good agreement for most sites, but the discrepancy at the Mt. Jeombong site, located in a high-altitude area,was higher than that at the other sites. It is estimated that the result reflected the spatial heterogeneity caused by micrometeorological factors such as terrain aspect, differences in elevation, solar radiation, cold air drainage flow,and so on, which were not accommodated in the MODIS-derived product due to the spatial resolution (1 km).

To verify the applicability of the AGDD threshold at the national scale,we compared the DoYAGDDwith the greenup DoY, which was extracted from the EVI based on curvature K for 30 sampling sites in Q. mongolica forests in South Korea (Appendix S2). The results of the comparison are shown in Fig. 8a. There was high consistency between DoYs extracted from both criteria (AGDD and curvature K). The coefficient of determination was 0.87,and RMSE was 3.9. The absolute differences between the two types increased with increasing altitude (Fig. 8b).

Discussion

Reconstruction process of satellite images to estimate air temperature

Traditionally, the estimation of air temperature has depended on ground measurements at point levels such as meteorological stations and ground surveys. The groundmeasured air temperatures are spatially interpolated using a GIS with conventional spatial statistics techniques such as kriging, inverse distance weighting (IDW), spline, and so on for expanding the estimation to the polygon level.Although the development of GIS and spatial statistics have led to a drastic refinement in the interpolated result,the limited number of points is a severe weakness.

Fig. 6 Threshold map of the number of days from 1 January 2015(days of year, DoY) needed to reach 159 AGDD (accumulated growing degree days) in 2015

Remote sensing provides an alternative data source since remotely sensed imagery is intrinsically spatialized(Neteler 2010).In particular,MODIS provides an abundant series of LST products with different spatial and temporal resolutions from both the Terra and Aqua platforms(Zhang et al.2013).Previous studies have shown that the LST data measured by MODIS can be successfully used for linear regression estimates of air temperatures at a regional scale(Mostovoy et al. 2006; Hassan et al. 2007). However, part of this work neglects the problem that cloud cover can obscure land and sea surfaces(Zhang et al.2013;Krehbiel and Henebry 2016). In the case of LST images, land/sea surface temperatures in cloud-covered areas are unavailable since cloud-top temperatures are measured instead(Neteler 2010). There are many ‘missing’ pixels in the filtered images after cloud-contaminated pixels have been removed (Fig. 3).

To solve the missing-pixel problem, we applied the algorithm by Neteler (2010) and improved by Zhang et al.(2013). At first, we temporally merged MODIS LST images from four acquisition times from the MODIS Terra/Aqua LST products.The relative close temporal proximity( ～3 h) of the Terra and Aqua overpasses provides an opportunity to combine information from the two data sources, which can reduce the data loss related to clouds(Crosson et al. 2012). The merged LST images were converted to air temperature images, then the images were spatially interpolated to the reconstructed air temperature images with auxiliary covariates by an ordinary cokriging algorithm. Microclimatic factors such as elevation, waterbody,and land-use intensity were faithfully reflected in the reconstructed images (Fig. 3).

Applicability of AGDD as an indicator of phenology

In general, temperatures affect plant functioning by influencing enzymatic activities (Bonhomme 2000). When the temperature is too low, the enzyme is not flexible enough and therefore unable to change conformation, but at high temperature, the protein can be denatured (Bonhomme 2000). In this study, we assumed a minimum temperature threshold of 5 °C for enzymatic activity in Q. mongolica based on studies of Kira(1945)and Schenker et al.(2014).Previous studies set up the maximum temperature threshold as above 40 °C, but we did not set up the threshold in this study because days that exceed 40 °C in daily mean air temperature during the green-up season of Q. mongolica are uncommon in South Korea.

In this study, we determined AGDD threshold of Q.mongolica forest based on ExG extracted from digital camera images and daily air temperature data measured from meteorological sensors from 2015 to 2016. The threshold agreed well with the green-up DoY determined from the EVI based on curvature K(Fig. 8).Therefore,the AGDD threshold can be regarded as a robust indicator of the green-up of Q. mongolica forest.

Another advantage when AGDD is used as a phenological indicator is that it is not influenced by species composition. At the satellite level, most phenological studies have depended on the detection of the phenological signal extracted from several vegetation indices, such as the normalized vegetation difference index (NDVI) or the EVI, calculated by the ratio of the visible light and nearinfrared light. However, the detection is often restricted when the target species is not dominant in the forest due to the limitation of spatial resolution. In this respect, the application of meteorological indicators can be an alternative. Because the AGDD threshold is simply a decision criterion for the phenological condition of a given plant species determined from the field study, it is not affected from factors such as the forest area, species composition,and presence or absence of dominant species that influence a vegetation index. Therefore, the AGDD threshold is still suitable, even when shrub or herbaceous species are selected as the target species to monitor phenology. Most importantly, however, the AGDD threshold applied in the

diagnosis and prediction of phenology must be determined from field experience with a number of repetitions.

Fig. 7 Comparisons of accumulated growing degree days (AGDD) derived from MODIS data and meteorological stations (meteostation) at a daily scale from DoY 1 to DoY 184. DoY:number of days from 1 January 2015

Fig. 8 Scatter plots for verifying DoYAGDD (predicted green-up day of year from accumulated growing degree days) with green-up DoY determined from field measurements (a) and the relationship between the absolute difference and elevation (b)

Conclusions

As interest in phenology increases as a diagnostic tool for climate change,studies of phenology on various scales are being carried out. With the advances of image processing technology, GIS and spatial statistics support these studies by enabling acquisition and analysis of information from various methods. However, since there are limitations inherent in each method, it is necessary to explore a solution to use them properly.In addition,the development or selection of appropriate phenological indices applicable to each method should precede.

In this study, we attempted to develop an algorithm to minimize contaminated data included MODIS LST image.Although the algorithm is based on previously developed methodology, we improved this method by combining it with a near-remote sensing(digital camera)technique.The results showed that this improved method reflects field information well. In addition, we verified that AGDD can be used as an efficient index for phenology monitoring.The AGDD threshold of plant phenology, confirmed through field-based research, can incorporate spatiotemporal variations in temperature and thus can be used as a strong and integrative tool for assessing climate change.

The results of this study demonstrate that the method can be an effective method for phenology study. Since MODIS images have been accumulating for about 20 years, we expect that the data can be used to track phenology patterns and trends in the past. However, there is a limitation in the requirement for an AGDD threshold of each target species for monitoring. Therefore, for increasing usability of the method, field experiments should be carried out to verify the AGDD information for various plant species.