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

Abstract

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\pm 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.

1. 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].

A new generation of GEO satellites has been launched over the past five years and includes the Himawari 8 [8], the Geostationary Operational Environmental Satellite (GOES) 16/17 [9], Fengyun (FY)-4 [10], and the Geostationary Korea Multi-Purpose Satellite-2A (GEO-KOMPSAT-2A) (https://www.wmo-sat.info/oscar/satellites/view/34). The first satellite of the Meteosat Third Generation program, a follow-on program of MSG, will be launched in 2021 [11]. Sensors onboard these satellites provide enhanced temporal, spatial, and spectral capabilities relative to previous generations of GEO sensors. For example, the temporal resolution of a full disk in the Himawari 8 Advanced Himawari Imager (AHI) is 10 minutes, and its spatial resolution is 0.5–1 km in the visible to NIR bands at sub-satellite points. The AHI was successfully used to improve the characterization of seasonal changes in terrestrial vegetation in the middle latitude region [12,13]. Top-of-atmosphere (TOA)/surface reflectance products for the sensors in the new generation GEO satellites, as well as higher level products, have been generated (e.g., [14,15,16,17,18]).

The high temporal resolution of GEO observations can complement LEO observations, which are influenced by cloud contamination [4]. SEVIRI vegetation products have been fused with data from a polar orbiting sensor, the Satellite Pour 1’Observation de la Terre (SPOT) VEGETATION for climate model analysis in an effort to fill in gaps identified in the SPOT VEGETATION data collected during the rainy season in West Africa [19]. The consistency between GEO and LEO products should be evaluated if data from GEO and LEO sensors are to be used synergistically.

The primary issue associated with combining GEO and LEO sensors is that the GEO and LEO angular conditions (the sun–target–sensor viewing geometries) do not agree except near the equator [15]. Mismatch between the angular conditions and the anisotropic features of the surface influences the consistency between the GEO and LEO sensor reflectances [20] and the downstream products, such as the vegetation indices, including the Normalized Difference Vegetation Index (NDVI). Comparative analyses of the GEO and LEO vegetation indices have been conducted previously (e.g., [18,21,22,23,24,25]). The NDVI values from the MSG SEVIRI and MODIS were compared [21]. Results showed a bias between the NDVIs that relied on geometric and land cover conditions. Numerous studies have compared the NDVI values computed using Bidirectional Reflectance Distribution Function (BRDF)-corrected reflectances calculated by employing Ross-thick/Li-sparse reciprocal (RTLSR) kernels [26,27]. Although the effects of differences between the angular conditions might have been mitigated by use of the BRDF model (e.g., [18,22]), BRDF-corrected reflectances are sensitive to accuracy in cloud mask used and aerosol contamination [22].

In addition, reflectances depend on a sensor’s spectral response function (SRF), which spectrally convolves the signatures of the surface and/or atmosphere, both of which depend on the land cover type. Therefore, reflectances and NDVI values obtained using different SRFs tend to show systematic differences [28]. Such differences have been mitigated using polynomial functions [28,29] and spectral band adjustment techniques [30]. Considerable effort has been devoted to tackling the problems associated with LEO–LEO NDVI intercomparisons/intercalibrations [31], and applications to actual data often require training data for fitting the NDVI transformation functions using polynomials or hyperspectral data for spectral band adjustments.

The NDVI can be transformed into the FVC (an areal measure of the vertical projection of green vegetation [32], which is essentially independent of the sensor characteristics) using the framework of a linear mixture model (LMM) [33]. The FVC may be obtained using an NDVI-isoline-based LMM [34,35] in which a target spectrum is modeled as a linear mixture of pure (vegetation and non-vegetation) spectra. The NDVI is then used to invert the FVC. Although the FVC is a function of the NDVI [36], the FVC should be independent of a sensor’s SRF if the pure spectra in the model “endmember spectra” can be appropriately parametrized.

The FVC cannot be obtained from an off-nadir view, such as the view used by a GEO to make observations in the middle-to-high latitude region because the directional FVC can change along the zenith angle view [37,38,39]. The directional fraction is less sensitive to the view angle for pixels in which vegetation shows a high degree of spatial segregation. In such a case, appropriate preparations 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 Section 6 and Section 7, respectively.

2. Test Site and Data

2.1. Test Site

Figure 1 highlights, in red parallelograms, the regions of interest (ROIs) used in numerical experiments across five regions in Japan, including Hokkaido (10 scenes), Tohoku (12 scenes), Tokai (14 scenes), Shikoku (16 scenes), and Kyushu (12 scenes), from northeast to southwest.

Table 1. Geographic coordinates used to explore the Moderate Resolution Imaging Spectroradiometer (MODIS) data and AHI view angles for each coordinate in degrees.

	Hokkaido	Tohoku	Tokai	Shikoku	Kyushu
Latitude	43	38.29	35.12	33.69	33.28
Longitude	141.38	140.83	137.38	133.49	130.34
AHI view zenith angle	49.6	44.3	40.9	39.9	40.3
AHI view azimuth angle	181.0	180.2	174.2	167.1	161.6

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.

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–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.

2.2. 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.

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. Date and timestamp for MODIS and AHI data used in this study. The timestamp is represented by four digits for hour and minute “hhmm”.

	Hokkaido	Tohoku	Tokai	Shikoku	Kyushu
Date and	16 July 2015, 0340	27 October 2015, 0345	19 December 2015, 0405	28 September 2015, 0415	31 July 2015, 0435
MODIS hhmm	8 May 2016, 0335	5 November 2015, 0340	17 March 2016, 0400	14 October 2015, 0415	3 October 2015, 0435
	12 August 2016, 0335	15 May 2016, 0340	4 May 2016, 0400	21 October 2015, 0420	19 October 2015, 0435
	2 May 2017, 0340	22 May 2016, 0345	7 July 2016, 0355	1 December 2015, 0415	4 November 2015, 0435
	18 May 2017, 0340	7 November 2016, 0340	31 August 2016, 0405	8 December 2015, 0420	16 January 2016, 0430
	12 Jun 2017, 0335	18 May 2017, 0340	12 November 2016, 0355	10 February 2016, 0420	23 May 2016, 0430
	14 July 2017, 0335	10 November 2017, 0340	30 December 2016, 0355	22 March 2016, 0415	30 May 2016, 0435
	21 May 2018, 0340	3 December 2017, 0345	31 January 2017, 0355	19 December 2016, 0415	19 February 2017, 0430
	6 Jun 2018, 0340	20 January 2018, 0345	11 March 2017, 0405	21 February 2017, 0415	2 Jun 2017, 0435
	24 July 2018, 0340	19 April 2018, 0340	28 April 2017, 0405	19 May 2017, 0420	18 Jun 2017, 0435
	-	26 April 2018, 0345	26 February 2018, 0405	4 Jun 2017, 0420	27 December 2017, 0435
	-	4 November 2018, 0345	14 March 2018, 0405	24 February 2018, 0415	27 April 2018, 0430
	-	-	30 March 2018, 0405	13 April 2018, 0415	-
	-	-	8 October 2018, 0405	29 April 2018, 0415	-
	-	-	-	13 October 2018, 0420	-
	-	-	-	29 October 2018, 0420	-
AHI hhmm	0340	0340	0400	0420	0430

. 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).

2.3. 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].

3. 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.

3.1. 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,

$$\begin{array}{c}\hfill {\widehat{\rho }}_{band}=\omega {\rho }_{v,band}+(1-\omega ){\rho }_{s,band},\end{array}$$

(1)

where ${\widehat{\rho }}_{band}$ is the modeled reflectance of an arbitrary $band$, $\omega$ is the FVC, ${\rho }_{v,band}$ and ${\rho }_{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)),

$$\begin{array}{cc}\hfill v& =\widehat{v},\hfill \end{array}$$

(2)

where

$$\begin{array}{cc}\hfill v& =\frac{{\rho }_n-{\rho }_r}{{\rho }_n+{\rho }_r},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \widehat{v}& =\frac{{\widehat{\rho }}_n-{\widehat{\rho }}_r}{{\widehat{\rho }}_n+{\widehat{\rho }}_r},\hfill \end{array}$$

(4)

and the red and NIR bands for $band$ are represented by r and n, respectively. Substituting Equation (1) into Equation (2) and solving for $\omega$ [35] yielded

$$\begin{array}{cc}\hfill \omega & =\frac{f_1(v,{\mathit{\rho }}_s)}{f_2(v,{\mathit{\rho }}_v,{\mathit{\rho }}_s)},\hfill \end{array}$$

(5)

where

$$\begin{array}{cc}\hfill f_1(v,{\mathit{\rho }}_s)& ={\rho }_{s,n}-{\rho }_{s,r}-v({\rho }_{s,n}+{\rho }_{s,r}),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill f_2(v,{\mathit{\rho }}_v,{\mathit{\rho }}_s)& =v({\rho }_{v,n}+{\rho }_{v,r}-{\rho }_{s,n}-{\rho }_{s,r})-{\rho }_{v,n}+{\rho }_{v,r}+{\rho }_{s,n}-{\rho }_{s,r},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\mathit{\rho }}_{em}& =[{\rho }_{em,r},{\rho }_{em,n}].\hfill \end{array}$$

(8)

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 ${\mathit{\rho }}_{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:

$$\begin{array}{c}\hfill \omega \approx \frac{f_1(v^{atm},{\mathit{\rho }}_s^{atm})}{f_2(v^{atm},{\mathit{\rho }}_v^{atm},{\mathit{\rho }}_s^{atm})},\end{array}$$

(9)

where

$$\begin{array}{cc}\hfill v^{atm}& =\frac{{\rho }_n^{atm}-{\rho }_r^{atm}}{{\rho }_n^{atm}+{\rho }_r^{atm}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathit{\rho }}_{em}^{atm}& =[{\rho }_{em,r}^{atm},{\rho }_{em,n}^{atm}],\hfill \end{array}$$

(11)

${\rho }_r^{atm}$ and ${\rho }_n^{atm}$ are the TOA red and NIR reflectances, and ${\mathit{\rho }}_{em}^{atm}$ is an endmember spectrum from the TOA reflectances of the endmember category $em$.

Equation (9) indicated that the FVC for the BOA, $\omega$, could be approximated using only the TOA measurement, $v^{atm}$, and the endmember spectra, ${\mathit{\rho }}_v^{atm}$ and ${\mathit{\rho }}_s^{atm}$. 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.

3.2. 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\pm p_2\%$ ($p_1=95,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 (${\widehat{\mathit{\rho }}}_v^{atm}=[{\widehat{\rho }}_{v,r}^{atm},{\widehat{\rho }}_{v,n}^{atm}]$).

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 ($\mathit{\rho }=[{\rho }_r,{\rho }_n]$) were then rotated $\theta =p_4\times \pi /180$ ($p_4=-30$) in the reflectance subspace to provide ${\mathit{\rho }}^{\prime }=[{\rho }_r^{\prime },{\rho }_n^{\prime }]$,

$$\begin{array}{cc}\hfill {\mathit{\rho }}^{\prime }& =\mathit{T}\left(\theta \right)\mathit{\rho },\hfill \end{array}$$

(12)

where

$$\begin{array}{cc}\hfill \mathit{T}\left(\theta \right)& =\left[\begin{array}{cc}\hfill cos\theta & \hfill -sin\theta \\ \hfill sin\theta & \hfill cos\theta \end{array}\right].\hfill \end{array}$$

(13)

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, $\tau$ (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, $\tau$ was set to $p_5$ (=0.04) in the quantile regression. The function in MATLAB® $fminsearch$ was used to compute the slope and intercept in the quantile regression [61]. The function $fminsearch$ identified the minimum without constraints using a derivative-free method. The resultant slope (${\beta }_1^{\prime }$) and intercept (${\beta }_0^{\prime }$) in the rotated reflectance space (Figure 3b) were restored to the original reflectance space,

$$\begin{array}{cc}\hfill {\beta }_1& =\frac{{\beta }_1^{\prime }-tan\theta }{1+{\beta }_1^{\prime }tan\theta },\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\beta }_0& =\frac{{\beta }_0^{\prime }}{cos\theta +{\beta }_1^{\prime }sin\theta },\hfill \end{array}$$

(15)

where ${\beta }_1$ and ${\beta }_0$ are the slope and intercept of the soil line-like equation in the original reflectance space.

The mean of the red reflectances (${\overline{\rho }}_r$) was computed after masking the pixels corresponding to water bodies. Likewise, the mean of the NIR reflectances (${\overline{\rho }}_n$) was extracted. A spectrum consisting of the computed mean values ($\overline{\mathit{\rho }}=[{\overline{\rho }}_r,{\overline{\rho }}_n]$) was then prepared (the yellow circle in Figure 3c). The intersection of the line passing through computed spectra (${\widehat{\mathit{\rho }}}_v^{atm}$ and $\overline{\mathit{\rho }}$) and the soil line was derived. The intersection was considered to be a pseudo-non-vegetation endmember spectrum (${\widehat{\mathit{\rho }}}_s^{atm}=[{\widehat{\rho }}_{s,r}^{atm},{\widehat{\rho }}_{s,n}^{atm}]$). This spectrum is indicated by the brown circle in Figure 3c and conforms to the relation:

$$\begin{array}{c}\hfill {\widehat{\mathit{\rho }}}_s^{atm}=\left[\frac{{\gamma }_0-{\beta }_0}{{\beta }_1-{\gamma }_1},\frac{{\gamma }_1({\gamma }_0-{\beta }_0)}{{\beta }_1-{\gamma }_1}+{\gamma }_0\right],\end{array}$$

(16)

where

$$\begin{array}{cc}\hfill {\gamma }_1& =\frac{{\widehat{\rho }}_{v,n}^{atm}-{\overline{\rho }}_n}{{\widehat{\rho }}_{v,r}^{atm}-{\overline{\rho }}_r},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\gamma }_0& ={\overline{\rho }}_n-{\gamma }_1{\overline{\rho }}_r.\hfill \end{array}$$

(18)

We used the pseudo-endmember spectra, ${\widehat{\mathit{\rho }}}_v^{atm}$ and ${\widehat{\mathit{\rho }}}_s^{atm}$ for ${\mathit{\rho }}_v^{atm}$ and ${\mathit{\rho }}_s^{atm}$ 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. Sensor- and scene-independent parameter settings for the endmember computation algorithm over ROIs in Japan.

${\mathit{p}}_1$ (%)	${\mathit{p}}_2$ (%)	${\mathit{p}}_3$ (%)	${\mathit{p}}_4$ (Degree)	${\mathit{p}}_5$ ($\mathit{\tau }$)
95	1	5	$-30$	0.04

; 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.

4. Data Processing and Comparison Methods

4.1. 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.

• Search near-nadir MODIS data: MODIS scenes that included a point $P({\lambda }_i,{\varphi }_i)$, where $\lambda$ and $\varphi$ 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 (${\theta }_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.
• 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({\lambda }_i,{\varphi }_i)$ or largest possible rectangular area that did not exceed 300 pixels on any side around $P({\lambda }_i,{\varphi }_i)$ was extracted. Subsequently, pixels in the data that exceeded 10 degrees of the viewing zenith angle were masked.
• 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.
• 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 illustrates the target area extraction and NDVI-index computation for the Tokai region for 14 March 2018. The black area corresponds to the view for which the zenith angle exceeded 10 degrees (left panel in 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.

4.2. 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®.

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_v^2$ and $R_{\omega }^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.

5. Results

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.

5.1. 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. Summary of the statistical analyses comparing the AHI and MODIS NDVIs as well as the AHI and MODIS NDVI-based indices over the Hokkaido region. N is the sample number of pixels used on each date. ${\theta }_a$ is the solar zenith angle for the AHI average value over the area on each date.

Date	N	${\mathit{\theta }}_{\mathit{a}}$	${\overline{\mathit{v}}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{v}}$	${\overline{\mathit{\omega }}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{\omega }}$	$\mathit{H}{0}_{\mathit{v}}$	$\mathit{H}{0}_{\mathit{\omega }}$	${\mathit{R}}_{\mathit{v}}^2$	${\mathit{R}}_{\mathit{\omega }}^2$
16 July 2015	1033	28	0.581	−0.0232	0.608	0.00569	1	0	0.772	0.779
8 May 2016	3601	32.5	0.378	−0.0205	0.454	0.00467	1	0	0.746	0.768
12 August 2016	1582	33.8	0.609	−0.0228	0.464	−0.0417	1	1	0.728	0.719
2 May 2017	1645	34	0.312	−0.0236	0.435	0.00914	1	0	0.854	0.862
18 May 2017	1433	30.5	0.473	−0.0141	0.486	−0.00294	1	0	0.71	0.696
12 Jun 2017	5873	27	0.664	−0.0253	0.652	−0.00767	1	0	0.865	0.866
14 July 2017	4034	27.9	0.589	−0.0359	0.474	−0.0945	1	1	0.842	0.82
21 May 2018	2240	30	0.497	−0.0189	0.581	0.00289	1	0	0.714	0.719
6 Jun 2018	1678	27.8	0.565	−0.0227	0.559	−0.0013	1	0	0.779	0.778
24 July 2018	1105	29.3	0.548	−0.0276	0.556	−0.00488	1	0	0.657	0.659

presents a statistical analysis of the NDVI/NDVI-based index differences. The averages of the AHI NDVI (${\overline{v}}_a$) were 0.31–0.66 and the differences between the average NDVI (${\overline{v}}_a$ minus the average MODIS NDVI), ${\delta }_v$, ranged from $-0.036$ to $-0.014$. The average AHI NDVI-based index (${\overline{\omega }}_a$) was 0.44–0.65, and the magnitude of the ${\overline{\omega }}_a$ minus the average MODIS NDVI-based index, ${\delta }_{\omega }$, ranged from 0.001 to 0.095.

Figure 8 shows density scatter plots of the red and NIR reflectances for the Hokkaido region on 16 July 2015, 12 August 2016, 14 July 2017, and 21 May 2018 as examples. 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 ${\delta }_{\omega }$ 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 ${\delta }_{\omega }$ ($-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. $H{0}_v$ and $H{0}_{\omega }$ were the null hypotheses for the NDVI and NDVI-based indices. The eighth and ninth columns in Table 4 listed $H{0}_v$ and $H{0}_{\omega }$, respectively. $H{0}_v$ was 1 for all dates whereas $H{0}_{\omega }$ 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 (${\overline{v}}_a$) was 0.41–0.59, and differences between the average NDVIs (${\overline{v}}_a$ minus the average MODIS NDVI) ranged from $-0.03$ to $-0.011$, as summarized in

Table 5. 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 Tokai region. N is the sample number of pixels analyzed on each date. ${\theta }_a$ is the solar zenith angle for the AHI data, averaged over the area on each date.

Date	N	${\mathit{\theta }}_{\mathit{a}}$	${\overline{\mathit{v}}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{v}}$	${\overline{\mathit{\omega }}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{\omega }}$	$\mathit{H}{0}_{\mathit{v}}$	$\mathit{H}{0}_{\mathit{\omega }}$	${\mathit{R}}_{\mathit{v}}^2$	${\mathit{R}}_{\mathit{\omega }}^2$
19 December 2015	4533	62.1	0.505	−0.0303	0.676	0.00934	1	1	0.918	0.916
17 March 2016	5890	40.6	0.442	−0.0234	0.613	−0.00708	1	0	0.952	0.95
4 May 2016	8115	27.6	0.584	−0.0214	0.722	−0.0111	1	1	0.903	0.905
7 July 2016	2100	22.3	0.588	−0.0107	0.689	0.00401	1	0	0.886	0.882
31 August 2016	1466	33	0.543	−0.015	0.633	−0.00812	1	0	0.927	0.93
12 November 2016	9066	57.8	0.519	−0.0278	0.686	0.0113	1	1	0.932	0.932
30 December 2016	5275	61.5	0.52	−0.0229	0.687	0.0172	1	1	0.92	0.918
31 January 2017	3918	55.3	0.405	−0.0249	0.557	−0.00347	1	0	0.943	0.942
11 March 2017	2478	42.8	0.525	−0.019	0.706	−0.000809	1	0	0.873	0.874
28 April 2017	1166	28.8	0.458	−0.0122	0.525	−0.0121	1	0	0.908	0.917
26 February 2018	6728	47.3	0.414	−0.0174	0.563	−0.00262	1	0	0.952	0.953
14 March 2018	7712	41.8	0.415	−0.0128	0.577	−0.00396	1	0	0.961	0.961
30 March 2018	13749	36.3	0.452	−0.0117	0.613	−0.00394	1	0	0.928	0.929
8 October 2018	4298	46.4	0.497	−0.024	0.592	−0.00904	1	0	0.946	0.948

. The average AHI NDVI-based index (${\overline{\omega }}_a$) was 0.53–0.72 and the magnitude of the difference between the average AHI NDVI-based indices (${\overline{\omega }}_a$ minus the averaged MODIS NDVI-based index) were less than 0.017. The $H{0}_v$ values were equal to 1 on all dates, whereas the $H{0}_{\omega }$ values were 0 except on four dates (19 December 2015, 4 May 2016, 12 November 2016, and 30 December 2016).

Figure 10 shows density scatter plots of the red and NIR reflectances for the Tokai region on 19 December 2015, 4 May 2016, 28 April 2017, and 8 October 2018 as examples, excluding the spectra of water bodies, over the Tokai region. The AHI reflectances were higher than the MODIS reflectances except between the months of April and August. In particular, the results for 19 December 2015 (Figure 10a,e), a date close to the winter solstice, showed large differences. The AHI NIR reflectances were 0.2–0.5 (Figure 10a) whereas the MODIS NIR reflectances were 0.1–0.3 (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

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. ${\theta }_a$ is the AHI solar zenith angle averaged over the area for each date.

Date	N	${\mathit{\theta }}_{\mathit{a}}$	${\overline{\mathit{v}}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{v}}$	${\overline{\mathit{\omega }}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{\omega }}$	$\mathit{H}{0}_{\mathit{v}}$	$\mathit{H}{0}_{\mathit{\omega }}$	${\mathit{R}}_{\mathit{v}}^2$	${\mathit{R}}_{\mathit{\omega }}^2$
27 October 2015	13802	55.3	0.508	−0.0364	0.67	0.0134	1	1	0.843	0.843
5 November 2015	17776	58.6	0.447	−0.0285	0.667	0.0313	1	1	0.823	0.826
15 May 2016	6691	26.7	0.547	−0.0248	0.602	−0.0187	1	1	0.853	0.855
22 May 2016	8180	25.7	0.542	−0.0328	0.591	−0.00602	1	1	0.896	0.897
7 November 2016	9021	59.2	0.433	−0.0238	0.664	0.0113	1	1	0.809	0.812
18 May 2017	3363	26.1	0.59	−0.021	0.652	0.0207	1	1	0.866	0.87
10 November 2017	6563	59.9	0.43	−0.0288	0.645	0.00181	1	1	0.856	0.844
3 December 2017	3422	63.9	0.377	−0.0247	0.626	0.00669	1	1	0.882	0.84
20 January 2018	2339	60.7	0.394	−0.0203	0.579	−0.0047	1	0	0.889	0.896
19 April 2018	7848	33	0.411	−0.0152	0.61	−0.00875	1	0	0.898	0.886
26 April 2018	5287	31.2	0.456	−0.0184	0.574	0.00475	1	0	0.893	0.89
4 November 2018	4067	57.8	0.49	−0.0314	0.674	0.0084	1	1	0.752	0.75

,

Table 7. 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 Shikoku region. N is the number of pixels sampled on each date. ${\theta }_a$ is the AHI solar zenith angle averaged over the area for each date.

Date	N	${\mathit{\theta }}_{\mathit{a}}$	${\overline{\mathit{v}}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{v}}$	${\overline{\mathit{\omega }}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{\omega }}$	$\mathit{H}{0}_{\mathit{v}}$	$\mathit{H}{0}_{\mathit{\omega }}$	${\mathit{R}}_{\mathit{\omega }}^2$
28 September 2015	4030	41.9	0.638	−0.0272	0.777	0.00706	1	1	0.715	0.709
14 October 2015	6969	47.6	0.619	−0.0227	0.765	0.0168	1	1	0.796	0.797
21 October 2015	6034	49.7	0.591	−0.0287	0.801	0.02	1	1	0.766	0.765
1 December 2015	1898	59.5	0.582	−0.0206	0.777	0.0232	1	1	0.76	0.761
8 December 2015	2886	60	0.612	−0.0265	0.656	−0.0053	1	0	0.634	0.648
10 February 2016	7814	50.9	0.586	−0.0174	0.725	0.0163	1	1	0.712	0.731
22 March 2016	13658	37.6	0.558	−0.0189	0.722	0.00808	1	1	0.818	0.822
19 December 2016	13150	60.7	0.536	−0.0145	0.707	0.0245	1	1	0.809	0.818
21 February 2017	2940	47.2	0.594	−0.0207	0.796	0.0169	1	1	0.793	0.804
19 May 2017	6053	23.5	0.688	−0.0194	0.79	−0.0366	1	1	0.683	0.687
4 Jun 2017	1744	21.4	0.679	−0.0156	0.796	−0.0227	1	1	0.766	0.757
24 February 2018	2264	46.6	0.57	−0.0151	0.787	0.0137	1	1	0.855	0.852
13 April 2018	9956	31.3	0.562	−0.0181	0.749	−0.00479	1	0	0.83	0.83
29 April 2018	7586	27.5	0.625	−0.019	0.786	0.00392	1	0	0.777	0.778
13 October 2018	1559	47.1	0.668	−0.0266	0.738	−0.00446	1	0	0.659	0.689
29 October 2018	1211	52.1	0.671	−0.0277	0.698	−0.0226	1	1	0.515	0.516

and

Table 8. 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 Kyushu region. N is the number of pixels sampled on each date. ${\theta }_a$ is the AHI solar zenith angle averaged over the area for each date.

Date	N	${\mathit{\theta }}_{\mathit{a}}$	${\overline{\mathit{v}}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{v}}$	${\overline{\mathit{\omega }}}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{\omega }}$	$\mathit{H}{0}_{\mathit{v}}$	$\mathit{H}{0}_{\mathit{\omega }}$	${\mathit{R}}_{\mathit{v}}^2$	${\mathit{R}}_{\mathit{\omega }}^2$
31 July 2015	1172	23.8	0.634	−0.00304	0.675	0.00426	1	0	0.938	0.937
3 October 2015	2883	43.3	0.621	−0.0221	0.711	−0.0107	1	0	0.916	0.911
19 October 2015	4611	48.9	0.579	−0.0248	0.707	−0.00463	1	0	0.856	0.865
4 November 2015	3613	53.9	0.527	−0.0251	0.671	0.00924	1	1	0.836	0.847
16 January 2016	8706	57.5	0.451	−0.0125	0.646	0.00295	1	1	0.841	0.841
23 May 2016	16654	23.9	0.625	−0.02	0.716	−0.0087	1	1	0.861	0.862
30 May 2016	4808	22.4	0.579	−0.0244	0.671	−0.0136	1	1	0.806	0.801
19 February 2017	16264	48	0.494	−0.0116	0.649	0.00363	1	1	0.846	0.847
2 Jun 2017	5111	22.2	0.593	−0.0147	0.699	−0.00344	1	1	0.856	0.859
18 Jun 2017	1008	21.3	0.419	0.000215	0.42	0.023	0	0	0.816	0.817
27 December 2017	4373	59.9	0.491	−0.0224	0.661	−0.00439	1	0	0.911	0.913
27 April 2018	1483	28	0.546	−0.00759	0.574	0.0108	1	0	0.835	0.834

, respectively. Trends in the statistics were similar to those observed over the Hokkaido and Tokai regions (Table 4 and Table 5). The null hypothesis for the NDVI-based index ($H{0}_{\omega }$) 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 ${\delta }_{\omega }$ computed from all ROIs was $-0.0004\pm 0.018$ (mean ± standard deviation) while ${\delta }_v$ was $-0.021\pm 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).

5.2. 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_v^2$, and the corresponding NDVI-based coefficients, $R_{\omega }^2$, were computed on each date and within each region, as summarized in the right-most two columns in Table 4, Table 5, Table 6, Table 7 and Table 8. The values of $R_v^2$ and $R_{\omega }^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

6.1. 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 ${\delta }_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 ${\delta }_{\omega }$ 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.

6.2. 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.

6.3. 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.

7. 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.