A Simple Algorithm for Deriving an NDVI-Based Index Compatible between GEO and LEO Sensors: Capabilities and Limitations in Japan

Geostationary (GEO) satellite sensors provide earth observation data with a high temporal frequency and can complement low earth orbit (LEO) sensors in monitoring terrestrial vegetation. Consistency between GEO and LEO observation data is thus critical to the synergistic use of the sensors; however, mismatch between the sun–target–sensor viewing geometries in the middle-to-high latitude region and the sensor-specific spectral response functions (SRFs) introduce systematic errors into GEO–LEO products such as the Normalized Difference Vegetation Index (NDVI). If one can find a parameter in which the value is less influenced by geometric conditions and SRFs, it would be invaluable for the synergistic use of the multiple sensors. This study attempts to develop an algorithm to obtain such parameters (NDVI-based indices), which are equivalent to fraction of vegetation cover (FVC) computed from NDVI and endmember spectra. The algorithm was based on a linear mixture model (LMM) with automated computation of the parameters, i.e., endmember spectra. The algorithm was evaluated through inter-comparison between NDVI-based indices using off-nadir GEO observation data from the Himawari 8 Advanced Himawari Imager (AHI) and near-nadir LEO observation data from the Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) as a reference over land surfaces in Japan at middle latitudes. Results showed that scene-dependent biases between the NDVI-based indices of sensors were −0.0004±0.018 (mean ± standard deviation). Small biases were observed in areas in which the fractional abundances of vegetation were likely less sensitive to the view zenith angle. Agreement between the NDVI-based indices of the sensors was, in general, better than the agreement between the NDVI values. Importantly, the developed algorithm does not require regression analysis for reducing biases between the indices. The algorithm should assist in the development of algorithms for performing inter-sensor translations of vegetation indices using the NDVI-based index as a parameter.


Introduction
Satellite remote sensing has facilitated biosphere monitoring by, for example, mapping terrestrial vegetation over the past two decades using polar-orbiting, low earth orbit (LEO) satellites [1][2][3]. Geostationary (GEO) satellites have also contributed to terrestrial vegetation monitoring. The Meteosat Second Generation (MSG) program includes four satellites launched over the period 2002-2015 and carrying the Spinning Enhanced Visible and Infrared Imager (SEVIRI). The sensor provides earth observation data across 12 spectral channels, including the visible and near-infrared (NIR) regions, with 3 km spatial sampling at sub-satellite points and a 15 min temporal resolution [4]. The SEVIRI data have been widely used for phenological monitoring [5], drought monitoring [6], and calculating vegetation products involving the leaf area index (LAI) and fraction of vegetation cover (FVC) [7]. of the view angle-dependent endmember spectra can provide a view angle-independent fractional abundance of the endmembers during spectral unmixing [40]. The model used in this study assumed such a situation. Below, the FVC is referred to as an "NDVI-based index" because the endmember spectra used in this study may not necessarily be pure, as described in more detail in Section 3. Our central concern is to demonstrate the sensor-(the geometric condition-and SRFs-) invariant NDVI-based index by utilizing appropriate endmember spectra using GEO and LEO sensor data.
This study was designed with two objectives in mind: (1) to develop an automated algorithm for computing sensor-and scene-dependent endmember spectra that can be used to derive NDVI-based indices, and (2) to inter-compare NDVI-based indices from the AHI reflectance (off-nadir view with backscattering, 40-50 degree view zenith) and those from collocated and coincident Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) reflectances (near-nadir view, ≤10 degree view zenith) over land surfaces in Japan at middle latitude. The TOA reflectances from both sensors that were not corrected for the BRDF effects were used, whereas the bottom-of-atmosphere (BOA) NDVI-based index was approximated. Near-nadir MODIS data were used as a reference. Inter-comparisons of the TOA NDVIs obtained from the AHI and Aqua MODIS were computed simultaneously.
Site and data information are summarized in Section 2, and the algorithm used to compute the NDVI-based index is presented in Section 3. Data processing and comparison methods are presented in Section 4. Section 5 provides the results of the inter-comparisons, followed by a discussion and conclusion in Sections 6 and 7, respectively.  Table 1 summarizes geographic coordinates used to explore the data in this study and AHI view angles for each coordinate in degrees. The locations and sizes of the areas depended on the observation date. A protocol for determining the areas is provided in Section 4.1.

Test Site
Regions selected in this study mostly belonged to the humid subtropical climate and humid continental climate in the Köppen-Geiger climate classification system [41]. Each region was selected to include various land cover types. The characteristics of the land cover types in these areas were investigated using the Terra and Aqua Collection 6 Land Cover Type Yearly L3 Global 5 km SIN Grid (MCD12Q1C) for 2015 [42,43], which was downloaded from the Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/). The 17-class International Geosphere-Biosphere Programme (IGBP) legend was available in the product. The class of "Water bodies" (class 17) was excluded prior to investigating the land covers. The percentage of each land cover was calculated on each date in each region. The "Open shrublands" (class 7), "Permanent snow and ice" (class 15) and "Barren" (class 16) classes were excluded from the analysis because their percentages were 0.0 in all regions. The land cover types for evergreen forests (evergreen needleleaf and broadleaf forests, classes 1 and 2 in the IGBP classification system) were categorized as "Evergreen forests". The types of deciduous forests (classes 3 and 4) were categorized as "Deciduous forests". The type for mixed forests (class 5) was left as is. Other types with vegetation (classes 6,[8][9][10][11][12]14) were categorized as "Non-forest vegetation", and urban and built-up areas (class 13) were denoted "Built-up". Figure 2 shows the overall average of each category using individually averaged percentages of land cover type of the ROIs for each region. Overall, the average percentages of evergreen forests, deciduous forests, and mixed forests were, respectively, 23.1%, 12.1%, and 26.4%, and non-forest vegetation and built-up were 26.0% and 12.4%, respectively. Kyushu regions in Japan, demarcated by red parallelograms, from northeast to southwest over the SRTM30 Plus from the NASA WorldWind WMS (https://data.worldwind.arc.nasa.gov/elev?) [44]. Yellow crosses indicate the points at which the MODIS scene was extracted in the first data extraction step. Their coordinates and view angles of Advanced Himawari Imager (AHI) for the points are summarized in Table 1.

Satellite Data
The polar-orbiting satellite, Aqua, has flown on an afternoon orbit since 2002, was originally known as NASA's Earth Observing System (EOS) PM-1, and carries MODIS. The Aqua MODIS Collection 6.1 Calibrated Radiances, the Daily L1B Swath 1 km (MYD021KM), and the Geolocation Fields Daily L1A Swath 1 km (MYD03) [45,46] were downloaded from the LAADS DAAC. TOA reflectances over bands 1 and 2 were calculated using calibration coefficients defined by the sun-to-Earth distance, in addition to the cosine of the solar incidence angle on the Earth target.
The Terra MODIS provided an alternative dataset; however, the magnitude of the relative azimuthal angle between the sun and AHI viewing angle on the Aqua overpass time, averaged over a year, was about 41-52 degrees for the target regions whereas the corresponding average angle on the Terra MODIS with a morning orbit was about 26-31 degrees. In this case, the AHI observations over the Aqua overpass time were less influenced by the backscattering effects [47][48][49] because the AHI viewing angle was more perpendicular to the principal plane. For this reason, the Aqua data was more suitable for our comparison. Overall average of each land cover category using individually averaged percentages of land cover type of the ROIs for each region based on the MCD12Q1C. The land cover type corresponding to evergreen forests (classes 1 and 2 in the IGBP classification system) was categorized as "Evergreen forests". The types for deciduous forest (classes 3 and 4) were categorized as "Deciduous forests". The type for mixed forests (class 5) was left as is. Other types with vegetation (classes 6, [8][9][10][11][12]14) were categorized as "Non-forest vegetation", and urban and built-up areas (class 13) were categorized as "Built-up".
Himawari 8, launched in October 2014, is located at 140.7 • E on the equator and covers Eastern Asia and the Pacific ocean to provide full-disk observation [8]. The full-disk data of AHI were downloaded from the National Institute of Information and Communications Technology (NICT) Science Cloud (https://sc-web.nict.go.jp/himawari/himawari-data-archive.html). The temporal resolution of the product was 10 min. Bands 3 and 4, corresponding to red and NIR bands, respectively, were used to compute the TOA reflectances using calibration information, the solar incidence angle, and the sun-to-Earth distance. The spatial resolutions of bands 3 and 4 were 0.5 and 1 km, respectively. The band 3 reflectances were then downscaled by the arithmetic average across 2-by-2 pixels to obtain 1 km resolution data.
The date and timestamp information of the scenes used in this study are summarized in Table 2. The data were selected by following the protocol in Section 4.1. The timestamp is represented by four digits for hour and minute hhmm and corresponds to the beginning of the full disk scan for AHI and acquisition time for the MODIS 5-min granule. The geographic coordinates of all areas extracted (red parallelograms in Figure 1) are shown in Supplementary Material (Tables S1-S5).

Geographic Coordinates and Illumination and View Angles
The MYD03 provides MODIS geolocation, illumination and view angle information. The geographic coordinates for AHI were computed using projection information in the AHI data and functions for the Normalized Projection Information [50]. The Global Multi-Resolution Terrain Elevation Data (GMTED2010) with spatial resolution of 30 arc-seconds were downloaded via EarthExplorer (https://earthexplorer.usgs.gov/) and used to obtain elevation data [51], which were then used to derive solar zenith and azimuth angles based on the library for the Full Vectorization of Solar Azimuth and Elevation Estimation [52] and AHI view angles based on the ellipsoidal model of the Earth [53].

Algorithm for Computing the NDVI-Based Index
The FVC using TOA reflectances were first derived in Section 3.1. We consider the FVC as the NDVI-based index in this study, because endmember spectra used to derive FVC were not necessarily pure.
All spectral variables in the next subsection depended on the sun-target-sensor viewing geometry. The angular variables were omitted here for brevity.

NDVI-Isoline-Based LMM Including Atmospheric Effects
A reflectance spectrum for the BOA with red and NIR bands was modeled by the two-endmember LMM,ρ whereρ band is the modeled reflectance of an arbitrary band, ω is the FVC, ρ v,band and ρ s,band are vegetation and non-vegetation endmember reflectances for band, and the subscripts v and s correspond to vegetation and non-vegetation. Again, the FVC is considered the NDVI-based index in this study. The NDVI from the target spectrum was assumed to be equal to the NDVI calculated from the modeled spectrum (Equation (1) where and the red and NIR bands for band are represented by r and n, respectively. Substituting Equation (1) into Equation (2) and solving for ω [35] yielded where ρ ρ ρ em = [ρ em,r , ρ em,n ].
The subscript em corresponds to either vegetation (v) or non-vegetation (s). The FVC calculated using Equation (5) was identical to the NDVI mixture model [32]. In addition, the FVC was identical to that calculated from the linear NDVI model [32], in which the NDVI could be modeled as a linear mixture of endmember NDVI values corresponding to vegetation and non-vegetation, if the 1-norms of the vegetation and non-vegetation endmember spectra were identical [54]. Unique values of ρ ρ ρ em should be computed using actual data for each scene or location.
Surface reflectances should be derived prior to computing the FVC in the application using the satellite data, although precise corrections to the effects of the atmosphere required ancillary data, such as the aerosol optical thickness. We then attempted to approximate the FVC using the TOA reflectances by simplifying the atmosphere-surface light interaction.
In this study, the FVC calculated using the TOA reflectances was derived by adding the atmospheric layer to the reflectances modeled by the LMM and by applying certain approximations.
The TOA reflectances were modeled as the sum of the path radiances and the first-and higher-order interaction terms between the atmosphere and the surface. The derivation of the following equation is described in Appendix A. The FVC could be approximated using the same form with Equation (5) but using the TOA variables, as follows: where v atm = ρ atm n − ρ atm ρ ρ ρ atm em = [ρ atm em,r , ρ atm em,n ], (11) ρ atm r and ρ atm n are the TOA red and NIR reflectances, and ρ ρ ρ atm em is an endmember spectrum from the TOA reflectances of the endmember category em.
Equation (9) indicated that the FVC for the BOA, ω, could be approximated using only the TOA measurement, v atm , and the endmember spectra, ρ ρ ρ atm v and ρ ρ ρ atm s . The advantage of the approximated equation was that no spectral parameters for the atmospheric layer (e.g., path radiances and atmospheric transmittances) were required to estimate the FVC if the assumption, that the atmospheric effects among all target spectra were equivalent, held. Note that the endmember spectra will be simulated using target spectra. The effects of approximating the atmosphere as homogeneous were expected to be small, as identified in the derivation presented in Appendix A, although this assumption may break down if there is atmospheric contamination by, for example, aerosol loading and water vapor, which is spatially heterogeneous.

Automated Computation of the Pseudo-Endmember Spectra
Endmember spectra of vegetation and non-vegetation for the TOA were computed using the characteristic distributions of the TOA reflectance spectra in the red-NIR reflectance space for each partial scene of each sensor. Below, the word "pseudo" is placed in front of each endmember spectra because the endmember spectra were not necessarily pure (100% vegetation or zero vegetation), especially within our ROIs. Finding pure spectra for a 1 km pixel size is expected to be challenging. For example, the pseudo-non-vegetation endmember may contain small amounts of vegetation. Use of pseudo-endmember spectra results in exceeding possible range of values of FVC, 0-1, and this is one of the reasons for using the term NDVI-based index.
The Soil Adjusted Vegetation Index (SAVI) [55] was used to compute the pseudo-vegetation endmember spectrum (green circle in Figure 3a). First, the reflectance spectra that computed the SAVI value p 1 ± p 2 % (p 1 = 90, p 2 = 1 in this study) from the bottom were extracted. The extracted spectra were sorted using the red reflectance values in ascending order. The p 3 % (p 3 = 5) of reflectance spectra from the smallest red reflectances were averaged to produce the pseudo-vegetation endmember spectrum (ρ ρ ρ atm v = [ρ atm v,r ,ρ atm v,n ]). The pseudo-non-vegetation endmember spectrum was computed using a soil line-like equation. Soil lines are defined by plotting the bare soil reflectances over the red and NIR reflectance space [56]. Many factors determined the slope and intercept, including the soil physical and chemical properties, soil moisture, soil mineral composition, roughness, etc. [57]. The equation is called a "soil line-like equation" due to the fact that any soil line is likely to be influenced by vegetation in the 1 km resolution data over our target region.
We included water bodies in the target area, and reflectance spectra rotation was implemented to obtain a stable soil line. The low reflectances of the water bodies provided a unique cluster of points in the reflectance space that was used to describe lines corresponding to soil lines. Water bodies were present because Japan is comprised of numerous islands with inland seas, and has extensive coastline on the Pacific Ocean. Spectral rotation resulted in stable lines corresponding to soil lines, especially over areas including pixels with NDVI < 0, corresponding to an unstable line if the rotation was not implemented. Reflectance spectra over a target area (ρ ρ ρ = [ρ r , ρ n ]) were then rotated θ = p 4 × π/180 (p 4 = −30) in the reflectance subspace to provide ρ ρ ρ = [ρ r , ρ n ], where In the present study, the soil line-like equation was automatically computed using a quantile regression [58,59]. The quantile regression fit the linear regression over the selected quantile, τ (0-1), to provide the response variable distribution [60]. This process was robust in the presence of extreme response variables compared to the mean-based classical regression [61]. This process was also useful for computing boundary lines in two-dimensional reflectance space [60]. The quantile regression has been widely used in soil line estimation [57,60]. The quantile level, τ was set to p 5 (=0.04) in the quantile regression. The function in MATLAB R f minsearch was used to compute the slope and intercept in the quantile regression [61]. The function f minsearch identified the minimum without constraints using a derivative-free method. The resultant slope (β 1 ) and intercept (β 0 ) in the rotated reflectance space (Figure 3b) were restored to the original reflectance space, where β 1 and β 0 are the slope and intercept of the soil line-like equation in the original reflectance space. The mean of the red reflectances (ρ r ) was computed after masking the pixels corresponding to water bodies. Likewise, the mean of the NIR reflectances (ρ n ) was extracted. A spectrum consisting of the computed mean values (ρ ρ ρ = [ρ r , ρ n ]) was then prepared (the yellow circle in Figure 3c). The intersection of the line passing through computed spectra (ρ ρ ρ atm v and ρ ρ ρ) and the soil line was derived. The intersection was considered to be a pseudo-non-vegetation endmember spectrum (ρ ρ ρ atm s = [ρ atm s,r ,ρ atm s,n ]). This spectrum is indicated by the brown circle in Figure 3c and conforms to the relation:ρ where We used the pseudo-endmember spectra,ρ ρ ρ atm v andρ ρ ρ atm s for ρ ρ ρ atm v and ρ ρ ρ atm s in Equation (9), which are geometric condition-, SRF-, and scene-dependent, for computing the NDVI-based index instead of the FVC. The values of the parameters, p 1 -p 5 , used to compute pseudo-endmember spectra in this subsection are listed in Table 3; these values were fixed for all AHI and MODIS scenes for our ROIs. The use of these constant values was important from a practical standpoint to avoid the need to adjust parameters scene-by-scene. Table 3. Sensor-and scene-independent parameter settings for the endmember computation algorithm over ROIs in Japan.

Processing Steps Used to Prepare the AHI and MODIS Scene Pairs
Pairs of nearly coincident and collocated AHI and MODIS partial scenes were prepared for the comparisons. Steps (1-4) were repeated across each region for the period between 7 July 2015 and 31 December 2018 to create pairs of MODIS and AHI data. The process flow is illustrated in Figure 4.

1.
Search near-nadir MODIS data: MODIS scenes that included a point P(λ i , φ i ), where λ and φ are the latitude and longitude and i identifies the region, designated in Table 1, were explored. The scenes that satisfied the following two conditions were obtained: View zenith angle of the point (θ v ) was equal to or less than 10 degrees, and the area around the target region was not influenced by clouds, cirrus clouds, or snow observed by visual inspection.

2.
Extract an area of at most 300-by-300 km of MODIS data: If the MODIS data satisfying the above condition were found, a 300-by-300 pixel area around the center point P(λ i , φ i ) or largest possible rectangular area that did not exceed 300 pixels on any side around P(λ i , φ i ) was extracted. Subsequently, pixels in the data that exceeded 10 degrees of the viewing zenith angle were masked. 3.
Extract a target area from the MODIS: A target area in the MODIS, selected to be as extensive as possible, was manually extracted from the 300-by-300 km data or the largest possible area of data. The shape of the area was a parallelogram that included more than 1000 non-water body pixels and did not contain cloud or snow. This was accomplished by visual inspection rather than by using cloud and snow flags, as this process almost completely avoided the effects of cloud contamination, including cirrus clouds and snow cover. The target area was selected to include water bodies, as discussed in Section 3.2.

4.
Extract the AHI target area and create a scene pair: The AHI scene observed at the time closest to the MODIS observation time was selected, and the AHI partial scene corresponding to the MODIS target area was extracted to create a MODIS and AHI partial scene pair for comparison.   Figure 5). The NDVI-based index, computed using the algorithm introduced in Section 3 using the MODIS data on the date, is shown in the center panel in Figure 5. The NDVI-based index computed using the nearly coincident AHI data collocated to the MODIS data is shown in the right panel in Figure 5. Pixels with no data in the NDVI-based index correspond to water bodies identified by the land/sea mask in the MYD03 product.

Comparison Method
In this study, the NDVI differences between AHI and MODIS using TOA reflectances were evaluated in addition to the NDVI-based index differences. Pixels corresponding to water bodies were excluded using the land/sea mask in the MYD03 product after pseudo-endmember computation for deriving the NDVI and NDVI-based index.
Differences between the NDVI/NDVI-based indices from two sensors were evaluated using a scene-by-scene comparison by describing their relative frequency distributions, computing the differences between the average NDVI/NDVI-based indices of the AHI and MODIS, and conducting a non-parametric statistical test (Wilcoxon rank sum test) of the frequency distribution to statistically identify whether the outputs from two sensors were extracted from distributions with the same median, against the alternative that they did not, using the ranksum of MATLAB R .
Pixel-by-pixel comparisons were conducted using the AHI and MODIS pairs. The effects of the lower spatial resolution of the off-nadir AHI data may have influenced the pixel-by-pixel comparison using the MODIS near-nadir data with a 1 km resolution. Uncertainties in the geometric calibration (relative geolocation errors) of the AHI and MODIS may also have impacted the comparison. The AHI and MODIS data were downscaled using the four-by-four pixel area average to mitigate the effects of the discrepancy on the actual spatial resolution and the relative geolocation errors. A nearest neighbor algorithm was then implemented to create a pixel-by-pixel pair for each AHI and MODIS pair. Coefficients of determination were computed using the results of the linear regression between the AHI (independent variable) and the MODIS (response variable) NDVI/NDVI-based indices (R 2 v and R 2 ω for NDVI and NDVI-based indices, respectively). The variability of the differences between NDVIs can be compared with that between NDVI-based indices using the coefficients of determination. Scatter plots of the NDVI/NDVI-based index for AHI and MODIS were also shown for visually investigating the variability. Figure 6 shows the illumination conditions and view angles of the Tokai region for 7 July 2016 and 19 December 2015, respectively, over the polar plot. Each angle is the average value of the angles calculated on that date. Points corresponding to the solar angle for other dates are near or between the points shown in Figure 6a,b. Other regions in Japan showed similar angular conditions. The MODIS and AHI solar zenith angles (circles) differed slightly due to small differences in the observation times. During the winter, the AHI view angle was closer to the principal plane ( Figure 6b) and the reflectances were sensitive to the backscattering effects that increase the reflectances. The backscattering effects were small if the solar zenith angle was small because the shadow effects decreased, especially during summer, as shown in Figure 6a.

Scene-by-Scene Comparison
The relative frequency distributions of the TOA NDVI and NDVI-based indices over the Hokkaido region on 10 dates are shown in Figure 7. Data from only the spring and summer seasons could be analyzed due to the presence of snow and cloud cover in other seasons. The NDVI distribution from the AHI data differed from that of the MODIS data (Figure 7a). The MODIS NDVIs were higher than the AHI NDVIs, especially at higher NDVI values (>0.3). The distributions of the NDVI-based indices for the two sensors were more similar than the NDVI distributions (Figure 7b). The NDVI-based index could assume negative values or values exceeding 1 because the NDVI of the pseudo-vegetation endmember was lower than the maximum NDVI, or the NDVI of the pseudo-non-vegetation endmember was higher than the minimum NDVI over the area. Table 4 presents a statistical analysis of the NDVI/NDVI-based index differences. The averages of the AHI NDVI (v a ) were 0.31-0.66 and the differences between the average NDVI (v a minus the average MODIS NDVI), δ v , ranged from −0.036 to −0.014. The average AHI NDVI-based index (ω a ) was 0.44-0.65, and the magnitude of the ω a minus the average MODIS NDVI-based index, δ ω , ranged from 0.001 to 0.095.   The spectra corresponding to water bodies were excluded prior to plotting. The magnitudes of the AHI reflectances were slightly higher than those of the MODIS. Differences between the view zenith angles of the AHI and MODIS did not significantly affect the results, unlike the results observed over the Tokai region, discussed further below. The small differences effect was attributed to the fact that all dates used over the Hokkaido region fell between May and August, during which the AHI solar zenith angle average within each area was small (27-34 degrees) at 1:30 local time (the equator crossing time of Aqua). For this reason, the effects of the view zenith angle on the density scatter plots were small.
The absolute values of δ ω for 12 August 2016 and 14 July 2017 in the Hokkaido region were higher (0.0417 and 0.0945) than those for other dates (Table 4). This result was attributed to differences in reflectance density distribution between sensors. For example, Figure 8g reveals two clusters of high density in the plot of the MODIS analysis. One of the clusters was superimposed on the vegetation endmember and another assumed smaller values of the NIR reflectances (marked by an arrow in Figure 8g). Figure 8c shows that the NIR reflectances of the cluster with smaller NIR reflectances were much higher (denoted by an arrow), and the NIR reflectances of the two clusters were nearly identical. The different characteristics of the reflectance density distributions resulted in inconsistent NDVI-based index values across sensors. For example, 4.7% of the AHI NDVI-based index for 14 July 2017 was less than −0.5 while 0.2% of MODIS NDVI-based index for the date was less than −0.5, leading to large biases in δ ω (−0.0945).
The rejection of the null hypothesis (H0), which is that the outputs from the two sensors were extracted from a distribution with the same median, corresponded to 1, whereas H0 returned 0 if the null hypothesis were not rejected. H0 v and H0 ω were the null hypotheses for the NDVI and NDVI-based indices. The eighth and ninth columns in Table 4 listed H0 v and H0 ω , respectively. H0 v was 1 for all dates whereas H0 ω was 0 for all dates except for 12 August 2016 and 14 July 2017.
The relative frequency distributions of the NDVI and NDVI-based indices over the Tokai region on 14 dates are shown in Figure 9. The AHI NDVI distribution was close to the MODIS distribution (Figure 9a), although slight differences were observed, especially at higher NDVI values (>0.5). The distributions of the NDVI-based indices of the two sensors were more similar than the NDVI distributions ( Figure 9b). The average AHI NDVI (v a ) was 0.41-0.59, and differences between the average NDVIs (v a minus the average MODIS NDVI) ranged from −0.03 to −0.011, as summarized in Table 5. The average AHI NDVI-based index (ω a ) was 0.53-0.72 and the magnitude of the difference between the average AHI NDVI-based indices (ω a minus the averaged MODIS NDVI-based index) were less than 0.017. The     (Figure 10e). This result was attributed to the anisotropy of the surfaces that is backscattering effects.  Statistical analyses of the differences between the AHI and MODIS NDVI as well as between the AHI and MODIS NDVI-based indices are summarized for the Tohoku, Shikoku, and Kyushu regions in Tables 6-8, respectively. Trends in the statistics were similar to those observed over the Hokkaido and Tokai regions (Tables 4 and 5). The null hypothesis for the NDVI-based index (H0 ω ) over the three regions was rejected at a higher rate than in the Hokkaido and Tokai regions. Nevertheless, the number of NDVI-based index rejections was always smaller than the corresponding value for the NDVI. The δ ω computed from all ROIs was −0.0004 ± 0.018 (mean ± standard deviation) while δ v was −0.021 ± 0.007.
The density scatter plots of Hokkaido region except for plots of four dates shown in Figure 8 are provided in Supplementary Materials (Figures S1 and S2 for AHI and MODIS), and the plots of the Tokai region except for plots of four dates shown in Figure 10 are provided in Supplementary Materials ( Figures S3 and S4 for AHI and MODIS). In addition, histograms of the relative frequency distributions of the NDVI and NDVI-based indices and the density scatter plots of the AHI and MODIS for Tohoku, Shikoku, and Kyushu regions are provided in Supplementary Materials (Figures S5-S16). Table 6. Summary of the statistical analyses showing agreement between the AHI and MODIS NDVIs as well as between the AHI and MODIS NDVI-based indices over the Tohoku region. N is the number of pixels sampled on each date. θ a is the AHI solar zenith angle averaged over the area for each date.

Pixel-by-Pixel Comparison
The variability of NDVI/NDVI-based index differences was evaluated by the pixel-by-pixel comparison using downscaled AHI and MODIS reflectance pairs. The NDVI between-sensor coefficients of determination, R 2 v , and the corresponding NDVI-based coefficients, R 2 ω , were computed on each date and within each region, as summarized in the right-most two columns in Tables 4-8. The values of R 2 v and R 2 ω were comparable on each date within each region. The values of the coefficients depended on the region, scene, and number of pixels included in each area. Figure 11 shows the scatter plots of the MODIS vs. AHI NDVIs and MODIS vs. AHI NDVI-based indices for the Hokkaido region (Figure 11a) and the Tokai region (Figure 11b). The results indicate that the variability of NDVI-based index differences was comparable to that of the NDVI differences while biases between NDVI-based indices of sensors were mitigated by using the developed algorithm, as shown in previous subsection. 6. Discussion

Differences between the Sensor-Specific NDVI/NDVI-Based Indices
The AHI NDVI calculated using the TOA reflectances (40-50 degree view zenith) displayed a negative bias (about 0.02) relative to the near-nadir MODIS TOA NDVI. The differences between the NDVIs were characterized using a scatter plot of the AHI and MODIS reflectances. Figure 12a shows the AHI and MODIS reflectances (blue and red dots) used in pixel-by-pixel comparisons over the Tokai region on 8 October 2018. The gray line revealed the trajectory of the reflectance change between the AHI and MODIS. The change of spectrum resulted primarily from surface anisotropy, but the SRF effects also contributed. The slope of the trajectory line was smaller when the red and NIR reflectances were similar; that is, the NDVI values were small. The slope gradually increased as the red reflectances decreased (increasing NDVI). The trajectory line slope resembled the NDVI isoline described by the dotted lines in Figure 12a, suggesting that in this case, the NDVI mitigated the effects of the sensor dependencies. The magnitude of δ v on this date was not large (0.024); however, the trajectories did not necessarily match the NDVI isoline on other dates or in other regions. Figure 12b shows an additional example over the Hokkaido region for 12 Jun 2017. The trajectory of the reflectance change differed from the NDVI isoline, for example, in areas with medium vegetation corresponding to the two ellipsoids in Figure 12b. In this area, the MODIS spectra (magenta ellipsoid) shifted almost vertically into the AHI spectra (cyan ellipsoid). Effects of relative geolocation errors manifested as anomalous trajectories. The errors in these trajectories may have been caused by geolocation errors in AHI data from the NICT showing temporal variations [62]. Geolocation errors of Aqua MODIS, by contrast, have been reported to be less than 60 m [63]. Although the NDVI might have mitigated the anisotropy effects, the SRF effects biased the NDVIs calculated by the two sets of sensors. The NDVI-based index computed by the algorithm developed here further mitigated the effects of the sensor dependencies. The value of the NDVI-based index is determined at the intersection of the NDVI isoline and the subspace of pseudo-endmember spectra described by yellow lines in Figure 12a,b [36]; hence, the geometric relationship between the target spectra and the simulated endmember spectra for AHI should be similar to those calculated for the MODIS data. Figure 12a shows that the geometric relationship between target and pseudo-endmember spectra for AHI is indeed similar to those for MODIS. In addition, the gradual increase in the slope of the trajectory lines (gray lines in Figure 12a) with the proper endmember spectra indicates that the NDVI-based index is similar across sensors even though the view zenith angles are different. This implies that vegetation in each pixel was likely spatially segregated. The magnitude of δ ω on a particular date was small (0.009). Caution should be taken when interpreting areas with medium vegetation (the actual FVC likely spanned 0.2-0.8), within which dramatic increases in the NIR bands were identified at larger viewing zenith angles (Figure 12b), as discussed above. The distance between the AHI target spectra (cyan ellipsoid) and the AHI pseudo-non-vegetation endmember (red filled circle) was larger than the corresponding MODIS distance (magenta ellipsoid and brown filled circles), which increased the AHI NDVI-based index relative to the MODIS index. Similar trends in the reflectance shift were identified at other ROIs in the Hokkaido region (e.g., Figure 8c,g), resulting in lower AHI NDVI-based index values, especially in pixels with a low NDVI-based index (<−0.5), as discussed in the previous section.
The relatively large increase in the NIR reflectances likely stemmed from changes in the directional fraction of the vegetation cover. Intuitively, visible soil surface areas between the canopies are expected to decrease at large view zenith angles [38], especially in areas that include, for example, scattered trees over artificial materials with 50% tree cover. The increases in the directional fraction of vegetation cover should increase the top-of-canopy NIR reflectances. This situation may have been present in several pixels. This work did not seek to correct this effect.
The characteristic differences between the NDVI-based sensor indices provided important information about where the LMM assumption could be applied and where the NDVI-based index, i.e., the FVC-like index from AHI was sufficiently consistent with the index based on the near-nadir MODIS data, within a certain uncertainty.

Comparison with Previous Studies
To the best of our knowledge, only a few studies have discussed differences in FVC data obtained using GEO and LEO sensors. A key work on this topic was a series of studies on the MSG SEVIRI FVC data products [7], which determine FVC using an algorithm based on the stochastic spectral mixture model. In that series of studies, the MSG SEVIRI FVC data were compared with those of LEO's FVC products, with SPOT VEGETATION data as a reference [7]. The SPOT FVC was based on a neural network that was trained to generate the best estimates of LAI, fraction of absorbed photosynthetically active radiation (FAPAR), and FVC from the fused and scaled MODIS and Carbon Cycle and Change in Land Observational Products from an Ensemble of Satellites (CYCLOPES) products [64]. Although the MSG SEVIRI and SPOT VEGETATION approaches yielded consistent results of FVC for a site in Gabon (Figure 11 in [7]), the algorithms implemented to retrieve FVC were not identical. Our intention in the present work was to compare FVCs (NDVI-based indices in this study) obtained using GEO and LEO sensors with an identical algorithm. Thus, the results of our study shed light on possible discrepancies in FVCs between the GEO and LEO sensors caused by the differences in geometric conditions and sensor specifications.
Most comparisons of data from GEO and LEO sensors have examined differences in reflectances and spectral vegetation indices such as the NDVI as representative parameters. Some of them performed corrections for atmospheric influences and BRDF effects [18,[22][23][24]65]. Those studies reported biases in NDVI. For example, the AHI NDVI computed from reflectances after atmospheric and BRDF correction was higher than the MODIS NDVI (the 16-day MODIS composite product (MOD13A2)) [18]. By contrast, the TOA MODIS NDVI was somewhat higher than the TOA AHI NDVI in this study. This discrepancy may be due to differences in correction levels and the test site. Specifically, uncorrected TOA reflectances were used in this study whereas corrected reflectances were used in previous works. Further studies are needed to investigate this issue.
In another study comparing NDVI data in the GEO-LEO framework, SEVIRI NDVI (atmospherically and BRDF-corrected) and MODIS NDVI based on the MCD43A4 16-day Nadir BRDF-Adjusted Reflectance (NBAR) product were compared [22]. The study found differences in NDVI between the two systems, which were attributed to sensor-dependent variations in the accuracies of the BRDF corrections as well as differences in SRF. In another study, NDVI data obtained from atmospherically corrected and BRDF-adjusted reflectances recorded by the Geostationary Ocean Colour Imager (GOCI) on the Korean Communication, Ocean and Meteorological Satellite (COMS) were compared with MODIS [23,65]. It should be noted that GOCI data used in the comparisons included the reflectance standardized to GOCI's own fixed view angle with the mean solar zenith angle [65]. Unlike the NBAR products, the view zenith angle was not adjusted to the nadir direction. The GOCI NDVI was compared with the MODIS NDVI computed from the nadir BRDF-adjusted reflectances provided by the BRDF model based on the RTLSR kernels that had been calibrated using a surface reflectance product (MOD09GA and MYD09GA) [65]. The study showed that the GOCI NDVI was lower than MODIS NDVI, with differences in the SRF being the suspected major cause of systematic errors in NDVI [23,65].
In general, inconsistencies between GEO and LEO NDVI values derive from differences in SRFs. In addition, uncertainties in the BRDF model (and its inversion process) also influence the accuracy of the BRDF correction. Moreover, such uncertainties are strongly influenced by the quality of the data, which can be degraded by cloud contamination, atmospheric effects, and/or poor scene geometry [22]. The algorithm developed in this study does not employ a BRDF model. For this reason, use of the proposed algorithm may reduce systematic errors between GEO and LEO vegetation products caused by BRDF and SRF effects.

Characteristics in the Developed Algorithm
The NDVI-based index for the BOA as a function of the TOA reflectances was approximated using the radiative transfer theory. The atmospheric model was homogeneous over the entire area of the single partial scene. Actual data were not homogeneous. The Rayleigh scattering depended on the atmospheric pressure specific to the elevation, and the aerosol loading and water vapor content were heterogeneous in the spatial domain. These effects may have added errors to the NDVI-based index [66]. The atmospheric effects required correction to mitigate the effects of the atmospheric heterogeneity in the subsequent comparison steps.
Significantly, training data, ancillary data, or statistical functions may not be required to obtain the view zenith angle-and SRF (sensor)-independent NDVI-based indices for areas where the LMM assumption can be applied in the context of the synergistic use of GEO and LEO sensors. It means that regression analysis is not required to obtain indices compatible between sensors. It is also important to note that the NDVI-based index was less sensitive to calibration uncertainties because the endmember spectra simulated using the algorithm developed here and the target spectra were both influenced by uncertainties, meaning that systematic components of the uncertainties in the endmember and target spectra canceled one another. This point should be further investigated.
The NDVI-based indices independently estimated from different locations could not be directly compared because the simulated endmember spectra were scene-dependent; however, small differences between the NDVI-based indices without external data would be beneficial for the development of algorithms for the inter-sensor translation of the NDVI. The NDVI-based index may be a parameter for relating vegetation indices of a reference sensor and a destination senor. Previously, soil brightness was used as a parameter to derive the inter-sensor relationships between vegetation indices [67], and the soil isoline equations were used to relate the vegetation indices [68]. Small differences between the NDVI-based AHI and MODIS indices thus demonstrated the applicability of the developed algorithm and suggested that the vegetation indices might be translated between sensors using the NDVI-based index as a parameter for GEO-LEO and/or LEO-LEO data fusion.

Conclusions
This study developed an algorithm for deriving the NDVI-based index that was designed to be compatible between GEO and LEO sensors. In the inter-comparison results, the NDVI-based indices calculated from AHI off-nadir observations (40-50 degree viewing zenith angles) were, in general, compatible with the MODIS near-nadir observations made in mid-latitude (Japanese) areas in which the vegetation in a pixel might have been spatially segregated, and the fractional abundances of each component were less sensitive to the view zenith angle. The algorithm developed here to automatically compute the pseudo-vegetation and non-vegetation endmember spectra mitigated the sensor dependencies, including the effects of differences between the viewing geometries and the sensor's SRFs over the areas. By contrast, the NDVI-based indices showed large biases in some areas. This effect can be attributed to the effects of sparsely scattered vegetation (e.g., scattered trees) within individual pixels. The directional fractions of vegetation canopy in such areas increased with increasing view zenith angle.
The NDVI-based index computed using the algorithm developed here showed better consistency in the relative frequencies relative to the TOA NDVI. In this study, the AHI TOA NDVIs displayed a negative bias relative to the MODIS TOA NDVIs (about 0.02) due to sensor dependencies.
The inter-comparison and integration of GEO and LEO data have been hindered by BRDF and SRF effects. As shown in several studies, one way to minimize these biases is to perform BRDF correction. However, this approach has the drawback that it requires prior knowledge of the BRDF model for each target. The accuracy with which the output of one sensor can be transformed to that of another sensor depends largely on the accuracy of the BRDF model. Moreover, minimization of the effects of differences in SRFs between the two sensors often requires statistical methods such as regression analysis. In general, the set of coefficients determined by the regression varies significantly according to land cover type. Furthermore, the coefficients are specific to each sensor pair. On the contrary, the algorithm developed in this study can derive comparable NDVI-based indices for GEO and LEO sensors without using regression analysis (training data). Because this novel approach does not require regressions between the two sensors, it has advantages over existing methods in which the set of coefficients determined by the regression depends on land cover type and the sensor pair.
The characteristics of the NDVI-based index differences over other areas, in the presence of atmospherically corrected reflectances, should be further investigated. Such studies will explore the applicability of fusing NDVI-based indices from GEO and LEO sensors. At least small biases in the NDVI-based indices calculated from different sensors over certain areas encourage the use of GEO sensors to complement LEO observations and the inter-sensor translation of vegetation indices for monitoring terrestrial vegetation.  Acknowledgments: AHI data were downloaded from the NICT science cloud. MODIS land cover data (MCD12Q1C) and Aqua MODIS calibrated radiances (MYD021KM) and their geolocation data (MYD03) were downloaded from the LAADS DAAC. The GMTED2010 data with spatial resolution of 30 arc-seconds were downloaded from EarthExplorer. Map services and data for the SRTM30 Plus are available from the NASA WorldWind WMS.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study, in the collection, analysis, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: where o em,band is the reflectance of the second-and higher-order interaction terms between the atmosphere and surface. Equation (A10) describing the vegetation endmember was equal to Equation (A2) if ω = 1, and Equation (A10) for the non-vegetation endmember was equal to Equation (A2) if ω = 0. The values of i (i = 1, 2) in Equation (A6) were expected to be small because the terms of each band (e.g., o n and −o s,n in Equation (A8)) canceled out their values if the values of the second-and higher-order interaction terms were comparable. Assuming that i in Equation (A6) was equal to zero, ω could be represented as follows: .
Numerical simulations to evaluate and validate these approximations may be presented in a separate study.