Evaluation of NDVI Estimation Considering Atmospheric and BRDF Correction through Himawari-8/AHI

Satellite-based vegetation indices are an essential element in understanding the Earth’s surface. In this study, we estimated the normalized difference vegetation index (NDVI) using Himawari-8/Advanced Himawari Imager (AHI) data and analyzed the sensitivity of products to atmospheric and surface correction. We used the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer model for atmospheric correction, and kernel-based semi-empirical bidirectional reflectance distribution function (BRDF) model to remove surface anisotropic effects. From this, top-of-atmosphere, top-of-canopy, and normalized NDVIs were produced. A sensitivity analysis showed that the normalized NDVI had the lowest number of missing values compared with the others and almost no low peaks during the study period. These results were validated by Terra and Aqua/Moderate Resolution Imaging Spectroradiometer (MODIS) and Project for On-Board Autonomy/Vegetation (PROBA) NDVI product, showing the root mean square error (RMSE) and bias of 0.09 and + 0.04 (MODIS) and 0.09 and − 0.04 (PROBA), respectively. These results also satisfied the FP7 Geoland2/BioPar project-defined user requirements (threshold: 0.15; target: 0.10).


Introduction
Vegetation index-based remote sensing is a standard tool of Earth system science and continues to play a central role in global change, carbon cycle, land coverage/use, and terrestrial ecology research. Among various vegetation indices that have been developed, the normalized difference vegetation index (NDVI) is used the most because it provides a good representation of vegetation characteristics over the long term. Moreover, other data products derived from remote sensing use NDVI as their primary input, including the leaf area index (Pontailler et al. 2003), global land cover maps (DeFries et al. 1999), fraction of absorbed photosynthetically active radiation (Ruimy et al. 1994), net primary production (Prince and Goward 1995), air temperature (Prihodko and Goward 1997), and land surface temperature (Jin 2004).
Many studies use three types of reflectance to calculate NDVI: top-of-atmosphere (TOA), top-of-canopy (TOC), and normalized reflectance with bidirectional reflectance distribution function (BRDF) modeling. TOA reflectance does not include atmospheric correction. In this case, atmospheric effects are removed during NDVI calculations (Holben 1986). The NDVI trend is similar, regardless of whether atmospheric correction (Chen et al. 2003). TOC reflectance includes atmospheric correction. Atmospheric elements such as water vapor, total ozone column, and aerosol optical depth (AOD) affect the accuracy of satellite-based land surface products (Nagol et al. 2009). The differences in the NDVI values occurs depending on the atmospheric correction method (Ke et al. 2015). Normalized reflectance with BRDF modeling incurs changes in the NDVI value with changes in the satellite angle and solar angle (Fensholt et al. 2006a). Additionally, the resulting NDVI products must account for the BRDF effect (Vermote et al. 2008); thus, the BRDF must be modeled according to the satellite's characteristics (Yeom and Kim 2013).
NDVI with TOA reflectance are retrieved by several organizations, including the National Oceanic and Atmospheric Administration (NOAA)/Advanced Very High Resolution Radiometer (AVHRR) and the Meteosat Second Generation (MSG)/Spinning Enhanced Visible and Infrared Imager (SEVIRI). The Satellite Pour l'Observation de la Terre (SPOT), Project for On-Board Autonomy/Vegetation (PROBA), Terra and Aqua/Moderate Resolution Imaging Spectroradiometer (MODIS), and the Suomi National Polar-orbiting Partnership (Suomi NPP)/Visible Infrared Imaging Radiometer Suite (VIIRS) have produced NDVI with TOC reflectance. Several of these organizations have geostationary satellites equipped with new sensors, such as the Geostationary Operational Environmental Satellite-16 (GOES-16), Meteosat Third Generation (MTG), and Himawari-8/9. All of these satellites have channels capable of estimating the NDVI.
South Korea operates a geostationary satellite, the Communication, Ocean, Meteorological Satellite (COMS), at the National Meteorological Satellite Center (NMSC); however, this satellite does not measure in the red/near-infrared (NIR) range required for vegetation index calculations. The Korea Meteorological Administration (KMA) launched the Geostationary Korea Multi-Purpose Satellite -2A (GK2A) with Advanced Meteorological Imager (AMI) sensors in December 2018. This satellite has a total of 16 channels that are capable of retrieving vegetation index information. However, the appropriate vegetation index for this satellite remains under investigation.
In this study, we attempted to identify a vegetation index suitable for GK2A by analyzing the sensitivity of the NDVI to atmospheric correction and BRDF effects. Atmospheric corrections were performed through Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer model (RTM) and anisotropic properties of the surface are considered with semi-empirical BRDF modeling. We then constructed three types of NDVIs based on TOA reflectance, TOC reflectance, and normalized reflectance considering BRDF effects. To compare them, we analyzed the temporal and spatial differences between NDVIs by land type.

Data
This study was conducted using Himawari-8/Advanced Himawari Imager (AHI) data, because it is difficult to obtain GK2A/AMI data at present. The satellite, developed as a nextgeneration multi-functional transport satellite (MTSAT), was launched on 7 October 2014 and is currently operated by the Japan Meteorological Agency (JMA).
The Himawari-8/AHI is located at longitude 140.7°E as a geostationary satellite. The satellite has a 16 channel multispectral imager to capture visible light and infrared images of the Asia-Pacific region and can provide full disk observations every 10 min and whole target area images every 2.5 min. In this study, we used the red channel (B3; central wavelength: 0.64 μm) and the NIR channel (B4; central wavelength: 0.86 μm). NDVI was calculated using data obtained over a 1-year period from January 2017 to December 2017.
We used various auxiliary data to calculate the NDVI (Table 1). Cloud mask data for aerosol products were provided by the JMA. Snow cover data applied by the algorithm were provided by the NMSC. The solar zenith angle (SZA), solar azimuth angle (SAA), viewing zenith angle (VZA), and viewing azimuth angle (VAA) for the angle calculation module was acquired from the COMS/Meteorological Imager (MI) (Paltridge and Platt 1976). Total ozone column, total precipitable water, and AOD were provided by the European Centre for Medium-Range Weather Forecasts (ECMWF)/Copernicus Atmosphere Monitoring Service (CAMS) and were used as daily product composites (Inness et al. 2019). The MODIS land cover type product (MCD12Q1) was applied to analyze the NDVI characteristics by indicator type.
We used MODIS NDVI and PROBA NDVI to compare the estimated NDVIs in this study. MODIS NDVI were obtained using the 16-day MODIS composite product (MOD13A2) from the Terra MODIS sensor is based on the Terra MODIS level 2 (L2G) daily surface reflectance product (MOD09 series). It provides red and NIR surface reflectance corrected for atmospheric gas, thin cirrus cloud, and aerosol effects. For data processing purposes, MODIS Vegetation Indexes (VIs) are generated in square tile units (~1200 km × 1200 km at the equator) and are mapped as sinusoidal grid projections (equal area projections).
Only tiles containing land features are processed. When mosaicked, all tiles cover the Earth land mass, allowing the global MODIS-VI to be generated every 16 days (Didan et al. 2015).
The data produced through this process have a spatial resolution of 1 km and a temporal resolution of 16 days (Huete et al. 1999).
The PROBA NDVI used in this study was the S10 NDVI 1-km product. PROBA-V was launched by the European Space Agency (ESA) on 7 May 2013. The S10 NDVI product of the satellite was synthesized using the maximum value composite (MVC) method over 10 days based on TOC reflectance data calculated using the simplified method for atmospheric correction (SMAC). The PROBA NDVI product can be rescaled to include quality information (status map) in the main NDVI layer. Pixels wrongly identified as "land" in the S10 NDVI status map are re-classified as "sea". The product was displayed in a regular latitude/longitude grid with the WGS 1984 (terrestrial radius: 6378 km). The resolution of the grid was 1/336°, with the center of the pixel as the reference. The spatial resolution is 1 km, and the temporal resolution was 10 days (Wolters et al. 2014).

Method
We estimated three type NDVIs based on TOA reflectance, TOC reflectance, and normalized reflectance. We compared the three NDVIs to determine the suitable NDVI for Himawari-8. First, we converted from Himawari-8/AHI channel DN data (Red, NIR) to TOA reflectance and TOA radiance using the calibration table. Only cloudy pixels were removed from the TOA reflectance data, as atmospheric correction could not be applied. Based on TOA reflectance data, we retrieved the TOA NDVI. The TOA NDVI was calculated at 1-h intervals to coincide with the 1-h interval of cloud data. The TOC reflectance was produced from TOA reflectance for atmospheric correction. The TOC NDVI was then estimated based on TOC reflectance. Finally, we calculated the NDVI based on normalized reflectance using the BRDF model. The NDVI-based normalized reflectance was calculated at intervals of 1 day.
To compare the data obtained using the three reflectance types, we analyzed temporal and spatial differences between the acquired NDVIs by land type. We chose an NDVI that showed consistent quality, without missing values or low peaks (sudden drops in the time series). From this, we selected the optimal NDVI representation and compared its results with those obtained from MODIS and PROBA NDVIs.
The following is a description of the atmospheric correction and BRDF used for the TOC reflectance and normalized reflectance calculations in this study.

Atmospheric Correction for TOC Reflectance
We used the 6S-based look-up-table (LUT) to correct for atmospheric effects because it can be applied to various satellites. The 6S RTM calculates the satellite bandwidth; the bandwidth is divided into 2.5-nm intervals over the shortwave area. This RTM has a high accuracy but requires a significant amount of computing time for TOC reflectance estimations of large areas (Zhao et al. 2001). To compensate for this, several studies have used the LUT method (Liang et al. 2001;Nunes et al. 2008). Here, we used the 6S-based LUT referred to by the GOES-16 / Advanced Baseline Imager algorithm for surface albedo (Liang et al. 2010); the contents are shown in Table 2.
If the LUT method is applied directly to atmospheric correction, discontinuities occur in the results. This is due to the discontinuity in the atmospheric correction coefficients that coincide with LUT intervals and is most prominent in SZA and VZA data in particular. Thus, in this study, we interpolated the SZA and VZA at 0.05°intervals.

BRDF for Normalized Reflectance
We used the semi-empirical BRDF model on Ross-Thick/Li-Sparse-Reciprocal kernels for estimating normalized reflectance (Roujean et al. 1992). Polar orbiting satellites observe the reflectance under various VZA conditions during the BRDF modeling synthesis period. In contrast, geostationary satellites such as the Himawari-8 observe the surface reflectance at fixed VZA conditions for each pixel. Therefore, in this study, the normalization method of reflectance was performed considering the characteristics of geostationary satellites. This method was proposed by Yeom and Kim (2013) as a modification of the VZA given by Duchemin et al. (2002). In this approach, the VZA is fixed to the VZA of the pixel, and the SZA and relative azimuth angle (RAA) are modulated by the average of the synthesis period to calculate normalized reflectance. This can be represented as follows: where θ s is the SZA, θ v is the VZA and ϕ is the RAA. where ρ norm is the normalized reflectance, ρ model (θ s = mean, θ v = θ v , ϕ = mean) is the adjusted reflectance at the VZA of the pixel and the mean SZA during the composite period, ρ measured is the TOC reflectance measured by the satellite, and ρ model (θ s , θ v , ϕ) is the calculated TOC reflectance using kernels with yielded empirical coefficients. The applied VZA is fixed to the VZA of the pixel, and the SZA and RAA are modulated by the average of the synthesis period to calculate the normalized surface reflectance. In this study, the synthesis period was set to 5 days. We used the mean value of the normalized reflectance calculated over 5 days to estimate the NDVI.

Himawari-8/AHI NDVI Spatial Distribution
We examined three NDVI types (TOA NDVI, TOC NDVI, and normalized NDVI) for the Himawari-8/AHI over a 1-year period from January to December 2017. The MVC values for the TOA and TOC NDVIs were synthesized over a 1-day period to equalize the temporal resolution of the three NDVIs.  Figure 1a-c shows the spatial distribution of the month mean TOA, TOC, and normalized NDVIs, respectively, for February 2017. Because the cloud mask data used in this study are calculated only up to 60°based on the center of the satellite image, the results of this study were also calculated according to the cloud area. Thus, the regions of India and northern Russia were not produced. The overall distribution of NDVIs was similar, regardless of the reflectance type used to calculate the NDVI. Barren areas with less vegetation showed only a slight difference in their NDVIs. The maximum NDVI difference is 0.07 in barren area. The TOA NDVI generally showed lower values than the other two, and the TOC NDVI was generally higher than the other NDVIs in all areas. The normalized NDVI tended to be distributed between the values of the TOA and TOC NDVIs. The average values are 0.23, 0.35, 0.26 for the TOA NDVI, TOC NDVI and normalized NDVI.
The TOA and TOC NDVIs are relatively large number of missing in several regions, such as South China, Vietnam, and Laos ( Fig. 2a and b). The normalized NDVI, in contrast, showed a smaller number of missing value pixels compared with the other two NDVIs (Fig. 2c). However, a small number of missing values occurred in areas near the equator, such as Papua New Guinea, where clouds continued to persist during the synthesis period (Kim et al. 2016).

Himawari-8/AHI NDVIs Time-Series Distribution According to Land Type
We analyzed the time-series distributions of the average NDVI to ascertain the NDVI trend among land types (Fig. 3). Six land types were included: grass, crop, shrub, mixed forest, evergreen forest, and barren. Figure 3 shows the NDVI value as a function of the estimation period. Three NDVI trends were calculated for the time series distribution; the calculation approach was similar to that used for the spatial distribution calculations. The NDVI values were calculated relatively higher in grass, crop, and mixed forest areas from May to August (Fig. 3a-d); this was attributed to the seasonal vegetation distribution, as there are more land areas in the Northern Hemisphere than in the Southern Hemisphere.
Conversely, NDVI values remained relatively constant over the 1-year period for shrub, evergreen forest and barren areas ( Fig. 3c-f). The shrub region almost represents the Australia's desert in MODIS land cover (MCD12Q1). For this reason, the shrub area also calculated consistent values.
The TOA NDVI showed the lowest values for all land types, with especially low values in areas where some vegetation was present, such as in evergreen forest, mixed forest, grass and crop areas (Fig. 3a-e). Evergreen forest is of particular interest, because it covers most tropical patches of Southeast Asia. The area is healthy throughout the year; nevertheless, the average TOA NDVI was 0.6, which did not reflect the characteristics of healthy vegetation.
In contrast, the TOC NDVI showed the maximum value for most land types and was relatively high in barren and shrub areas ( Fig. 3c and f), which have sparser vegetation distributions. The 6S RTM estimated high coefficient values with high SZA and VZA; as a result, the NDVI values were high. We then synthesized the NDVI using the MVC method; high The normalized NDVI showed tendencies similar to those of the TOC NDVI in mixed forest and evergreen forest areas where vegetation is abundant (Fig. 3d and e). The TOA NDVI and the normalized NDVI were similar in shrub and barren areas ( Fig. 3c and f), where vegetation is sparse; here volumetric scattering rarely occurs. A pattern in which the TOA NDVI and TOC NDVI decreased sharply then increased again (low peak) was apparent for most land types. These trends were determined to be due to an error that had arisen from unremoved clouds, surface moisture caused by rain, and other factors (Kim et al. 2011). In contrast, the normalized NDVI did not significantly change over the time series, regardless of land type. Therefore, normalized NDVI, with few missing values or low-peak occurrences, provided the optimal NDVI model (Gao et al. 2002).

Comparison of Other Satellite Products
The final calculated normalized NDVI was compared with the widely used MODIS and PROBA NDVIs. We performed a normalized NDVI 16-day MVC synthesis for a comparison with the MODIS NDVI and a 10-day synthesis for a comparison with PROBA NDVI data. Both datasets were tested using homogeneous pixels. Based on a 3 × 3 array of validation pixels, the standard deviation was 0.03 or less.   Figure 4 shows the difference in NDVI values according to land type. The differences are expressed in quartiles. The line in the black box represents the median. The difference between MODIS and PROBA was confirmed before their comparison with the normalized NDVI (Fig. 4a). The difference between each NDVI values occurred inevitably due to differences in atmospheric correction models, spectral response functions and geometric features (Fensholt et al. 2006b). MODIS NDVI showed lower values than PROBA NDVI for all land types. The difference was relatively small in barren and shrub areas. However, the difference was more significant in evergreen forest and mixed forest areas. Median values were the most similar for grass areas. The normalized NDVI was high for all land types compared with the MODIS NDVI (Fig. 4b). The smallest difference appeared in the shrub area; this area pertains mostly to Australia. During the study period, this region had the least cloud-induced impact compared with other regions; thus, the data quality was better for this region, as indicated by the small difference between indices. Similar trends were observed in barren areas. However, low NDVI values in desert regions also have an impact. The biggest difference appeared in evergreen forest. Most of the region resides in Southeast Asia near the equator. The area has significant cloud cover; thus, there is a high probability that the pixel quality is reduced, due to an incomplete cloud mask. However, the median value was smallest for the mixed forest area. The overall difference in all indicator areas was less than the difference between the validation data set. The normalized NDVI was lower for all land types compared with PROBA NDVI data (Fig. 4c). Overall, the trend was similar to the interquartile range between the two validation data sets, i.e., a small interquartile range for the shrub area (0.04) and large interquartile range for evergreen (0.14~0.18) and mixed forest (0.12~0.16) areas. Figure 5 shows the difference between the normalized NDVI and the validation data set according to month. We did not perform comparisons between validation data because the temporal resolution was not match. There was no significant difference between the normalized NDVI data and the validation data sets in a comparison of monthly data. Compared with MODIS and PROBA NDVI data, the normalized NDVI data resided in the middle, being relatively high with respect to MODIS data and relatively lower than PROBA data over all months. This trend was the same over all land types. Table 3 shows the calculated root mean square error (RMSE) and bias of the validation data sets. The average RMSE was 0.10, and the bias was 0.06. The black circles in Fig. 5 show the results of the comparison between the normalized NDVI and MODIS NDVI; 23 days in total were considered for the comparison. The average RMSE was 0.09, and the bias was 0.04. Overall, similar tendencies were observed; however, the RMSE and bias were smaller than the values given in Table 3.
The white circles in Fig. 6 show the results of a comparison between the normalized NDVI and PROBA NDVI. The average RMSE for all dates was 0.09, and the bias was −0.04. The RMSE was similar to other comparison results. However, the absolute value of the bias decreased and showed a negative trend.
NDVI is not an essential climate variable, and there are no Global Climate Observation System specifications on required accuracy. The FP7 Geoland2/BioPar project defined user requirements for NDVI in terms of acceptable differences with existing satellite-derived products. The requirements specify a threshold RMSE of 0.15 and a target RMSE of 0.10 (Swinnen and Toté 2017). Regardless of the type of validation data set, the RMSE was <0.10 in this study. Other studies related to NDVI estimation have calculated the RMSE to be between 0.10 and 0.15 (Yeom and Kim 2013;Ke et al. 2015).
We compared the NDVI over a 1-year period for the areaaverage representing each of the six land types (barren, mixed forest, crop, and grass), as shown in Fig. 7. In this graph, the red circle indicates the normalized NDVI, the blue triangle represents the MODIS NDVI, and the green rhombus shows the PROBA NDVI. The overall trend of each NDVI was similar over the time series. Unlike other NDVIs, the normalized NDVI do not perform temporal synthesis. PROBA NDVI had relatively high values over all seasons and land types. The MODIS NDVI and the normalized NDVI exhibited similar values. September and October correspond to autumn and winter seasons in the Northern Hemisphere. During this time, the NDVI decreased sharply by 0.2~0.3, as shown in Fig. 7d. MODIS NDVI and PROBA NDVI data based on polar orbit satellites also showed a similar trend in some of the data (two or three); however, the normalized NDVI acquired more than 30 data points during this time period (Fensholt et al. 2011;Yeom and Kim 2015).

Concluding Remarks
In this study, we analyzed the sensitivity of atmospheric correction and BRDF effect to determine an NDVI suitable for the GK2A/AMI using Himawari-8/AHI data. A TOA NDVI, TOC NDVI, and normalized NDVI were compared to determine the optimal NDVI for GK2A/AMI operations. The results of our spatiotemporal analysis revealed the normalized NDVI to have the smallest number of missing values compared with the other two NDVIs, and there were almost no low peaks during the time series. The normalized NDVI calculated in this study was compared with the widely used MODIS and PROBA NDVIs. We calculated an RMSE of ≤0.1 for the normalized NDVI, regardless of the validation data set. This result satisfies FP7 Geoland2/BioPar projectdefined user requirements and demonstrates the high accuracy of this NDVI estimation approach.
However, the normalized NDVI showed missing values in some areas. In the present study, the NDVI was calculated using 1-h interval data. The problem of missing values may be resolved if cloud data become available at 10-min intervals in the future. In addition, the reanalysis data of ECMWF was used to perform the current atmospheric correction. We expect that high accuracy can be achieved using satellite data with high spatial and temporal resolution in our study area.
The normalized NDVI calculated in this study can be applied to the future GK2A for geostationary-based NDVI in Asia and Australia; here, the data have a high temporal resolution, as opposed to that available from polar orbit satellites. We expect these findings to be useful for monitoring and mitigating changes in the Earth's surface in the near future.