Development and Demonstration of a Method for GEO-to-LEO NDVI Transformation

: This study presents a new method that mitigates biases between the normalized difference vegetation index (NDVI) from geostationary (GEO) and low Earth orbit (LEO) satellites for Earth observation. The method geometrically and spectrally transforms GEO NDVI into LEO-compatible GEO NDVI, in which GEO’s off-nadir view is adjusted to a near-nadir view. First, a GEO-to-LEO NDVI transformation equation is derived using a linear mixture model of anisotropic vegetation and nonvegetation endmember spectra. The coefﬁcients of the derived equation are a function of the endmember spectra of two sensors. The resultant equation is used to develop an NDVI transformation method in which endmember spectra are automatically computed from each sensor’s data independently and are combined to compute the coefﬁcients. Importantly, this method does not require regression analysis using two-sensor NDVI data. The method is demonstrated using Himawari 8 Advanced Himawari Imager (AHI) data at off-nadir view and Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) data at near-nadir view in middle latitude. The results show that the magnitudes of the averaged NDVI biases between AHI and MODIS for ﬁve test sites (0.016–0.026) were reduced after the transformation ( < 0.01). These ﬁndings indicate that the proposed method facilitates the combination of GEO and LEO NDVIs to provide NDVIs with smaller differences, except for cases in which the fraction of vegetation cover (FVC) depends on the view angle. Further investigations should be conducted to reduce the remaining errors in the transformation and to explore the feasibility of using the proposed method to predict near-real-time and near-nadir LEO vegetation index time series using GEO data.


Introduction
The normalized difference vegetation index (NDVI) has been used for monitoring terrestrial vegetation from the regional to the global scale [1], where the NDVI functions as a proxy indicator of biophysical parameters involving the leaf area index (LAI) and the fraction of absorbed photosynthetically active radiation (FAPAR) [2,3], among other parameters. The global NDVI time series of the Advanced Very High Resolution Radiometer (AVHRR) onboard the National Oceanic and Atmospheric Administration (NOAA) polarorbiting, low Earth orbit (LEO) satellite series has been provided since the early 1980s [4,5]. The AVHRR dataset can be integrated into or continued by LEO satellite sensors with advanced sensor attributes, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua platforms and the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi National Polar-Orbiting Partnership (NPP) spacecraft, for long-term vegetation monitoring [6,7].
Next-generation geostationary (GEO) satellites for Earth observation began operation in the 2010s and have provided highly calibrated data from optical sensors with improved temporal, spatial, spectral, and radiometric resolutions compared with the previous generation of GEO sensors [8]. For example, the Advanced Himawari Imager (AHI) onboard Himawari 8 launched in 2014 and the Advanced Baseline Imager (ABI) onboard Geostationary Operational Environmental Satellites (GEOSs) 16 and 17 launched in 2016 and 2018, respectively, can provide full-disk observation data in 10 min intervals and at 0.5-2 km spatial resolution [9,10]. An important aspect of the GEO data for land surface monitoring is its high temporal resolution, which enables observation of diurnal variations of reflectances and radiances as well as the mitigating effects of cloud contamination. The GEO NDVI time series has therefore been used for improved monitoring of land surface phenologies in middle latitude forests [11,12] and for more accurate characterization of the seasonality of the Amazon evergreen forest [13].
Combining the GEO NDVI dataset with the historical, global LEO NDVI dataset increases the spatiotemporal coverage of satellite data, introducing the possibility of advanced vegetation monitoring. Before the GEO and LEO NDVI are combined, their consistency should be evaluated. GEO-LEO NDVI inter-comparisons have thus been reported in numerous studies [11][12][13][14][15][16]. In the inter-comparisons, the reflectances between GEO and LEO sensors showed biases due to the sensor-dependent characteristics of the observations (e.g., [17]). The primary cause of the biases is the coupled effects of surface anisotropy and differences between the illumination and viewing geometric conditions of the sensors. A GEO sensor observes land surfaces from a fixed, oblique direction with a large variation of the solar angle, whereas LEO sensors such as MODIS and VIIRS observe the Earth from various viewing angles along the scan direction, with a certain variation of the illumination condition specific to the location and the satellite overpass time [14]. Reflectance biases between GEO and LEO sensors caused by the geometric condition result in biases between NDVIs.
The reflectance biases can be reasonably avoided by selecting only relative-azimuthalangle-matched reflectances [17], which have been used in radiometric inter-comparisons at forest and urban areas. For vegetation analysis using GEO and LEO reflectance time series, bidirectional reflectance distribution function (BRDF) models were fitted, and the fitted models were used to standardize reflectances to identical geometric conditions [15,16,18]. However, reflectance biases between sensors still remained after the BRDF corrections, possibly as a result of the effects of calibration uncertainties and discrepancies in the band configuration and spectral response function (SRF) [15][16][17][18]. These effects can be mitigated using regression analysis [19] or spectral band adjustments and radiometric intercalibration [20], which require additional data such as hyperspectral surface reflectances and/or atmospheric parameters [21,22].
In another relevant study involving inter-comparisons of GEO and LEO NDVIs [23], the effects of the viewing zenith angle and SRF differences between Himawari 8 AHI (off-nadir view) and Aqua MODIS (near-nadir view) were mitigated at regional scale at middle latitude. More specifically, AHI and MODIS NDVIs were converted to the fraction of vegetation cover (FVC) using endmember spectra automatically computed from each sensor's data. The FVC values of AHI and MODIS were nearly equivalent except for a few cases involving a viewing angle dependency. The results imply that the FVC can be used to relate the NDVIs of GEO at off-nadir view and LEO at near-nadir view and that a GEO-to-LEO NDVI transformation is feasible [23], which should facilitate the combination of GEO and LEO NDVIs with smaller biases.
The objective of the present study is to develop a new method that geometrically and spectrally transforms GEO NDVI directly to LEO-compatible GEO NDVI at near-nadir view by following the insights reported in [23]. First, a GEO-to-LEO NDVI transformation equation is derived based on a linear mixture model of anisotropic vegetation and nonvegetation endmember spectra. Using the derived equation, an NDVI transformation method is developed. Notably, this transformation does not require regression analysis using two sensors' NDVI because endmember spectra, which are used to determine the coefficients of transformation equations, are computed from each sensor's data indepen-dently. The method is numerically demonstrated using Himawari 8 AHI and Aqua MODIS top-of-atmosphere (TOA) reflectance data at middle latitude sites.

Background
A parameter necessary to derive the NDVI transformation equation, the directional FVC, is first defined. A general definition of the FVC, which differs from the directional FVC, is the vertically projected area of vegetation canopy to a unit area. It can be derived using scaled vegetation indices (pixel dichotomy model) [24], spectral mixture analysis [25], machine learning approaches using a radiative transfer model (RTM) simulation [26], or existing FVC products [27].
The directional FVC was assumed to be one minus a viewable gap fraction (a fraction of land surface that is visible from above) [28][29][30][31]. The directional FVC changes with the viewing angle (Figure 1), and the directional FVC at nadir view corresponds to the general definition of the FVC. In the present study, reflectances are expressed by a sum of the vegetation and nonvegetation reflectances showing surface anisotropy with the directional FVC for a fraction of vegetation reflectances: whereρ band(i,j) is a modeled reflectance for band; i corresponds to the direction of the solar angle specific to the solar zenith and azimuth angles; j corresponds to the direction of the sensor (viewing) angle specific to the sensor zenith and azimuth angle; ω (j) is the directional FVC; ρ v,band(i,j) is the vegetation reflectance; ρ s,band(i,j) is the nonvegetation reflectance (subscripts v and s denote vegetation and nonvegetation (soil), respectively); ρ v,band(i,j) and ρ s,band(i,j) vary with the viewing angle and with the sun angle condition that changes sunlit and shaded parts over vegetation canopy and nonvegetation surfaces. A spectrum that consists of multiple bands (two bands in the present study) of vegetation or nonvegetation reflectances is referred to as an endmember spectrum. The endmember spectra and the directional FVC are sensitive to the structure of the canopy and to the density and spatial distributions of canopy stands when the angular condition varies [32][33][34], which is, however, not explicitly expressed in Equation (1). The directional FVC for a pixel can be derived by minimizing the Euclidean distance between modeled and target spectra [25], where endmember spectra for a scene are extracted prior to the minimization. Similarly, the directional FVC can be derived by minimizing differences between NDVIs from modeled and target spectra with endmember spectra for the scene [35,36]. In a previous study, the latter approach was used to derive the FVC for spatially uniform atmospheric conditions [23]. The directional FVC specific to viewing directions can be expressed by and v atm where ρ atm r(i,j) and ρ atm n(i,j) are the TOA, partially atmosphere-corrected or fully atmospherecorrected reflectances of the target spectrum and subscripts r and n correspond to red and NIR bands, respectively; ρ ρ ρ atm em(i,j) is the TOA, partially atmosphere-corrected or fully atmosphere-corrected endmember spectra for vegetation and nonvegetation (em corresponds to v or s). The terms 1(i,j) and 2(i,j) are functions of the target NDVI and the secondand higher-order interaction terms between the atmosphere and surfaces. The details of Equation (2) are shown in Appendix A.
The values of 1(i,j) and 2(i,j) are equal to zero when the reflectances are atmospherically corrected. The values of 1(i,j) and 2(i,j) for TOA-level or partially atmospherecorrected reflectances are expected to be close to zero because the individual terms of the red and NIR bands in 1(i,j) and 2(i,j) tend to cancel out [23]. Thus, the expression form of the equation for directional FVC as a function of TOA or partially atmosphere-corrected reflectances is identical to that as a function of fully atmosphere-corrected reflectances when 1(i,j) = 2(i,j) = 0. Note that, in general, the value of the directional FVC tends to be sensitive to atmospheric correction levels unless correct endmember spectra are chosen for a pixel such that the Euclidean distance between the modeled and target spectra becomes zero.

Derivation of the NDVI Transformation Equation
Directional FVCs of two sensors are defined, and then the FVCs are used to relate NDVIs of different sensors. The resultant NDVI transformation equation is a function of the endmember spectra for two sensors.
Equation (2) is reformulated to simplify the derivation as follows, where Two sensor types are defined: an original sensor that is a sensor-to-be-transformed, and a destination sensor. The directional FVC for each sensor is expressed by ω sen(j sen ) = g v atm sen(i sen ,j sen ) , ρ ρ ρ atm sen,v(i sen ,j sen ) , ρ ρ ρ atm sen,s(i sen ,j sen ) , sen,1(i sen ,j sen ) , sen,2(i sen ,j sen ) , where the subscript sen in Equation (10) and in the following equations corresponds to a (sensor-a, original sensor) or b (sensor-b, destination sensor). The FVC is used as a parameter to relate the NDVI of two sensors (v atm a(i a ,j a ) and v atm b(i b ,j b ) ). The following equation with an error term ( ω(j a ,j b ) ) is assumed: where ω(j a ,j b ) represents an error included in ω a(j a ) with reference to ω b(j b ) and is sensitive to the structure of the canopy and to the density and spatial distribution of the canopy stands. An NDVI transformation equation describing the inter-sensor NDVI relationship is derived using Equations (5)- (11), in which geometric variables (i sen and j sen ) are omitted for brevity: where A flow diagram of the derivation is illustrated in Figure 2.

Overview of the Algorithm
A basic NDVI transformation algorithm is illustrated in Figure 3. First, endmember spectra of the original GEO sensor (ρ ρ ρ atm a,v and ρ ρ ρ atm a,s ) are computed using the automated algorithm of endmember computation which is described in Section 4.2; the NDVI of the original sensor (v atm a ) can also be computed (step 1). Second, endmember spectra of the destination LEO sensor (ρ ρ ρ atm b,v and ρ ρ ρ atm b,s ) are computed using the same endmember algorithm (step 2). Finally, endmember spectra for each sensor and the GEO NDVI are used to estimate the LEO-compatible NDVI (v atm b ) (step 3). The NDVI computed using LEO data (v atm b ), which is not shown in Figure 3, can be used to evaluate the transformed NDVI here denoted byv atm b to distinguish it from the NDVI denoted by v atm b . In our numerical demonstration described in Section 6, this evaluation is carried out with the values of ζ 1 and ζ 2 set to zero. Real applications of the time-series analysis using the method are discussed in Section 8.2.

Automated Endmember Computation
Vegetation and nonvegetation endmember spectra for original and destination sensors are necessary to transform the NDVI using Equation (12). For this transformation, the present study requires a single set of endmember spectra for each sensor scene. A representative approach to endmember extraction is to use vertices of scatter plots for the target area over the reflectance space [37] or feature space [38]. In a previous study, a similar idea was used for computing vegetation and nonvegetation endmember spectra using red-NIR reflectance scatter-plot information [23]; this approach is also used in the present study. Instead of finding vertices, the algorithm automatically finds a point along the left side boundary of the red and NIR reflectance scatter plot for computing the vegetation endmember spectrum and a point along the right side boundary (soil line) for computing the nonvegetation endmember spectrum. Note that the endmember spectra for the NDVI transformation equations are merely parameter vectors and need not be pure [23], which is discussed in Section 8.2.

Test Sites
Land surfaces at middle latitude where the viewing zenith angle of AHI shows middle values-specifically, approximately 40 • -50 • -were explored. Among these land surfaces, that of Japan (Figure 4a) was identified as suitable for our demonstration for the following reasons: (1) Several climatic zones are included because Japan is located in a transitional zone between sub-boreal to sub-tropical [39], and (2) large variations of reflectances along seasons are expected because of Japan's mountainous topography and higher solar zenith angle at LEO (Aqua) overpass time in winter, which is a consequence of Japan being located at true north of the sub-satellite point of AHI. The solar zenith angle variation in our regions of interest (ROIs) was approximately 19-64 • (area-averaged value), which likely caused variations of the areas for shaded parts. These characteristics make our demonstrations more stringent. Descriptions of the ROIs used in the present study are introduced in Section 5.3.
The maps in Figure 4b-f show the five test sites: Hokkaido, Tohoku, Tokai, Shikoku, and Kyushu (from north to south, 31-45 • N, 129-143 • E). Dotted rectangles or parallelograms in the maps are polygons for the ROI, which vary with the date of observation. The geographic coordinates for each site used for exploring data (the center of each site) and the viewing angles of AHI for the coordinates are summarized in Table 1. The Hokkaido and Tohoku ROIs are located primarily in a cool temperate zone; the other three ROI sites are dominated by a warm temperate zone. The compositions of the land cover in the ROIs differ across sites, as shown in Figure 4g, which was created using the MODIS land cover product (MCD12Q1C) [40]. Evergreen forests in Figure 4g include, for example, coniferous species such as Japanese cedar (Cryptomeria japonica) and Hinoki cypress (Chamaecyparis obtusa) and broadleaf species such as Japanese chinquapin (Castanopsis sieboldii) and evergreen oaks (e.g., Quercus glauca). Deciduous forests include deciduous oaks (e.g., Quercus serrata and Quercus crispula) and Japanese beech (e.g., Fagus crenata), among others [39,41,42]. Nonforest vegetation in Figure 4g includes rice paddy, crop, and grasslands [43].

Satellite Data
The Himawari 8 is located at 140.7 • E, and the AHI onboard the Himawari 8 observes 16 bands from the visible to the infrared region. The AHI performs a full disk observation with 10 min intervals. The Center for Environmental Remote Sensing (CEReS) at Chiba University distributes Himawari 8 AHI data, which are gridded over the geographic coordinates [44,45]. Red-band data (band 3) are gridded with a 0.005 • interval, and blue, green, and near-infrared (NIR)-band data (bands 1, 2, and 4, respectively) are gridded with a 0.01 • interval. These data were downloaded from the CEReS ftp site (ftp://hmwr829gr. cr.chiba-u.ac.jp/, (accessed on 27 August 2021)). The AHI TOA reflectances for each band were computed using the apparent reflectance values, solar zenith angle, and the sun-to-Earth distance [23]. The AHI red band (band 3) reflectances were aggregated arithmetically to obtain 0.01 • data.
The MODIS onboard Aqua flies an afternoon orbit, and its Equator crossing time is 1:30 p.m. The Collection 6.1 Calibrated Radiances, the Daily L1B Swath 1 km (MYD021KM) [46], and the corresponding Geolocation (MYD03) [47] were downloaded from the Level 1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) website (https://ladsweb.modaps.eosdis.nasa.gov/, (accessed on 27 August 2021)). Note that the MODIS data were limited to data including the test sites observed from the near-nadir direction. The MODIS TOA reflectances for each band were computed using scaled integer values, calibration coefficients, and the solar zenith angle.

Partial Scene Pair for ROIs
The ROI data for the five sites are the same as those used in [23] except that the AHI data were converted from data from the National Institute of Information and Communications Technology (NICT) Science Cloud into CEReS gridded data because the gridded data contain fewer geolocation errors [45]. Nearly coincident and collocated AHI and MODIS data were used, and the date range of the data spanned between July 2015 and December 2018. The number of pairs of AHI and MODIS partial scenes for the Hokkaido, Tohoku, Tokai, Shikoku, and Kyushu sites were 10, 12, 14, 16, and 12, respectively.
The ROIs were extracted to have as many pixels as possible from 300 × 300 km 2 areas centered on the coordinates for each site (Table 1); during the extraction, pixels corresponding to a MODIS zenith viewing angle larger than 10 • or cloudy pixels were excluded. The number of pixels used for analysis was greater than 1000 for every partial scene. The coordinate information for all of the ROIs is provided in the supplementary materials in [23].

Numerical Demonstration Procedure
To numerically demonstrate the NDVI transformation method, we evaluated differences between MODIS NDVI at near-nadir view and MODIS-compatible AHI NDVI at nearnadir view (transformed AHI NDVI), where the MODIS NDVI was used as a reference.
For the five test sites shown in Section 5, the variables ζ 1 and ζ 2 in Equation (12) were assumed to be zero in the numerical experiments. The reasoning for this assumption is explained as follows. Similar frequency distributions of the directional FVCs (ω sen(j sen ) in Equation (10)) from the TOA AHI (off-nadir view) and TOA Aqua MODIS (near-nadir view) data were observed at our test sites [23]. Although a few scenes exist in which the FVCs show a strong viewing angle dependency, the assumption that ω in ζ 1 and ζ 2 are zero during the NDVI transformation is mostly reasonable. In addition, sen,1 and sen,2 in ζ 1 and ζ 2 were set to be zero because 1(i,j) and 2(i,j) for each sensor would be close to zero. These assumptions enable us to omit ζ 1 and ζ 2 from the NDVI transformation equation.
Scene-level and pixel-by-pixel comparisons of NDVIs were conducted. MODIS and AHI pixels that include water bodies identified by the land/sea mask in the MYD03 product were excluded from our analysis of the NDVI differences after endmember computations [23]. In the scene-level comparisons, differences between area-averaged values of AHI NDVI or transformed AHI NDVI and MODIS NDVI were computed for all partial scenes. Their averages for each site, denoted by m, were computed as follows: and v atm a,m,k = where d orig,m is the average of the differences between the area-averaged AHI and MODIS NDVIs for site m; d trans,m is the average of the differences between the area-averaged transformed AHI and MODIS NDVIs; v atm a,m,k,n , v atm b,m,k,n , andv atm b,m,k,n are the AHI, MODIS, and transformed AHI NDVIs for the n-th pixel of the k-th date at site m, respectively; K m is the number of the date relying on the site; and N a,m,k and N b,m,k are the number of pixels relying on the date and site for the AHI and MODIS, respectively.
In the pixel-by-pixel comparison, NDVIs obtained at 0.01 • or 1 km spatial resolution were aggregated at 4 × 4 pixels using the arithmetic mean for mitigating the effects of geolocation errors and pixel size differences. Each pixel of an aggregated MODIS NDVI was paired with an aggregated original or transformed AHI NDVI using a nearest-neighbor algorithm. Linear regression analysis was conducted to evaluate the results obtained using these NDVIs for each partial scene. In addition, for each site, the aggregated MODIS NDVI for all dates were combined to generate a single MODIS NDVI dataset. Similarly, a single NDVI dataset of AHI or transformed AHI at each site was generated. Note that each pixel in the AHI or transformed AHI dataset is paired with one of the pixels in the MODIS NDVI dataset. The AHI, MODIS, and transformed AHI NDVI of the -th pixel for the dataset of site m were denoted by v atm * a,m, , v atm * b,m, , andv atm * b,m, , respectively. The differences between v atm * a,m, and v atm * b,m, (∆v orig,m, ) and betweenv atm * b,m, and v atm * b,m, (∆v trans,m, ) were computed for each pixel: The mean of the absolute values for Equations (26) and (27) in each site, referred to as the mean absolute difference (MAD) of NDVI, and the standard deviation (SD) of Equations (26) and (27) in each site were also computed. In the following equations, the subscript state corresponds to either the original AHI (orig) or the transformed AHI (trans), and we have and where µ state,m is the MAD for site m and σ state,m is the SD for site m. Parameter S m is the total number of the pixel pair dependent on site m.

Scene-Level Evaluations of NDVI Transformation
Results for relative frequency distributions of TOA NDVIs for AHI and MODIS are shown in Figures 5 and 6. Figure 5 corresponds to the "best date" for each site, where the absolute value of the averaged AHI (transformed) minus the MODIS NDVIs was minimal among all dates. By contrast, Figure 6 corresponds to the "worst date", which provides the maximal absolute value of the averaged AHI minus the MODIS NDVIs. The left column in Figure 5 shows similar distributions of the AHI and MODIS NDVI, although slight differences are observed, especially at higher NDVI values. Thus, the MODIS NDVI was slightly higher than the AHI NDVI. The right column in Figure 5 shows a much more similar distribution of the AHI (transformed) and MODIS NDVI, indicating successful transformations of the NDVI. The left column in Figure 6 shows the same trend as the left column in Figure 5; however, the right column in Figure 6 shows slight under-and over-corrections (e.g., Figure 6d,f), which are not observed in the right columns in Figure 5.  The mean of the averaged AHI or transformed AHI minus MODIS NDVIs for each site (d orig,m and d trans,m ) is summarized in Table 2. The original AHI NDVIs were approximately 0.02 smaller than the MODIS in NDVI units. Biases between the transformed AHI and MODIS NDVIs were much smaller (the absolute value of the mean was smaller than 0.01).  This coefficient reflects nonlinear effects of the NDVI transformation equation and simply η b minus η a . The η b for Shikoku is higher than the corresponding η a in almost all cases, indicating that inter-sensor relationships between the NDVIs for Shikoku could be slightly more nonlinear than those for the other sites.  Figure 7. Averages for the four coefficients in the NDVI transformation equation (p 0 , p 1 , q 0 , and q 1 ) for the five test sites. (a-d) correspond to p 0 , p 1 , q 0 , and q 1 , respectively. The error bars correspond to the standard deviation.

Pixel-by-Pixel Evaluations of the NDVI Transformation
Aggregated NDVIs were used for pixel-by-pixel evaluations. Figure 8 shows the spatial distribution of the NDVI differences between sensors for the partial scene of the Shikoku site on 13 April 2018. Figure 8a shows the spatial distribution of the aggregated MODIS NDVI. Negative biases could be identified in numerous parts in the original NDVI differences (Figure 8b). The large errors are attributable to sensor noises and to random errors caused by remaining relative geolocation differences, atmospheric and cloud contamination differences due to differences in the observation times (within approximately 5 min), and unidentified error sources. Overall negative errors were reduced or turned into small positive errors after the transformations (Figure 8c). Density scatter plots of AHI NDVI vs. MODIS NDVI for the data (Shikoku, 13 April 2018) are shown in Figure 9a, and plots of the transformed AHI NDVI vs. MODIS NDVI are shown in Figure 9b. These plots indicate that the AHI NDVI approached the one-to-one line after the transformations, indicating that the biases were mitigated. The results for other data for Shikoku and other sites typically showed the same trend, although substantial biases of the NDVI were observed in specific pixels after the transformations. Figure 9c,d show the results for Hokkaido on 14 July 2017. After the transformation, negative biases were observed at certain points where the MODIS NDVIs were less than approximately 0.5 (Figure 9d). These negative biases might be attributable to the fact that directional FVCs from two sensors showed biases over areas that included sparsely distributed vegetation. In such areas, the light interception by the vegetation canopy may have sharply increased with increasing viewing zenith angle, thereby increasing the NIR transmittances and reflectances of vegetation as well as the top-of-canopy reflectances. The increase in the AHI NIR reflectances subsequently resulted in FVC biases between the sensors [23]. The FVC biases would have propagated into the transformed AHI NDVI in such cases. Finally, Table 3 summarizes the MAD (µ orig,m and µ trans,m ) and the SD (σ orig,m and σ trans,m ). The MAD decreased after the transformations, whereas the SD slightly increased. The relatively large increase in the SD in Hokkaido is attributable to the directional FVC biases between the sensors. Previous studies regarding vegetation index (VI) inter-comparisons between GEO and LEO sensors have typically applied BRDF correction/adjustment of reflectances to mitigate the effects of surface anisotropy [15,16,18,48]. Two sensors' VIs derived using the corrected/adjusted reflectances were then used to compare phenologically significant transition dates [11]. In the present study, our intention was to reduce the absolute NDVI differences between sensors. A few studies quantitatively investigated the absolute NDVI differences using BRDF-adjusted reflectances (e.g., [16,48]); for example, Tran et al. [16] reported that mean absolute differences between atmospherically corrected, near-solar-noon AHI NDVI and atmospherically corrected and BRDF-adjusted MODIS NDVI (adjusted to solar noon and an oblique AHI viewing angle) was only 0.018-0.030 over Australian grassland/pasture sites. The remaining errors between the VIs are likely attributable to discrepancies in pixel footprints, geolocation accuracy, cloud mask methods, atmospheric correction methods, SRFs, and calibration uncertainties. The method developed in the present work is distinct from the common approaches using the BRDF method. In particular, our method is capable of standardizing the GEO NDVI to the near-nadir viewing direction and mitigating the effects of differences in the SRFs.
In the present study, we used the directional FVC as a parameter to derive the NDVI transformation equation on the basis of the equivalence of the FVCs between sensors with the error term. Similar approaches assuming the equivalence of physical parameters for deriving inter-sensor VIs or reflectances equations have been reported [22,[49][50][51][52][53][54]. To quantitatively estimate the coefficients of the equations, which are functions of physically meaningful parameters, nonlinear regression was successfully applied for global applications [52]. On the contrary, coefficients of NDVI transformation equations could be obtained without using regression analysis. For example, atmospheric parameters and hyperspectral profile of surfaces [22] or pixel-by-pixel minimum and maximum NDVIs from smoothed NDVI products [55] were prepared and used to compute the coefficients for NDVI conversion equations. Compared to these, our method does not require any prior information, and coefficients of the transformation equation could be automatically computed using reflectance data.

Benefits and Limitations of the Developed Method
The method developed in the present study can not only be used to geometrically and spectrally transform the NDVI but may be less sensitive to several biases included in the data. This reduced sensitivity is attributable to the endmember spectra used to compute the coefficients of the transformation equation and the target spectra both being influenced by sensor-dependent biases (e.g., biases in gain coefficients for radiometric calibration) and, thus, directional FVCs tend to be sensor-independent, leading to correct transformations of the NDVI. The advantage of using such "image endmembers" for reducing the effects of biases included in the data has also been reported in studies for spectral mixture analysis [56] as well as relative radiometric normalization between different sensors [57].
The NDVI transformation reduced biases between NDVIs in our results; however, the variability of the NDVI differences remained or slightly increased, possibly because of the viewing angle dependency of the FVC and errors that cannot be corrected in the transformation. The errors likely arise from observation time differences between the sensors (within approximately 5 min [23]) causing cloud and atmospheric condition differences, pixel footprint differences, geolocation errors, and intrinsic sensor noise. Thus, consideration of the directional FVC biases, i.e., ω in the transformation, and improvements in the data quality (e.g., sensor noise and geolocation errors) would further reduce the variability of the NDVI differences after the NDVI transformation [51]. Additionally, to mitigate the effects of differences in FVC as a function of direction in numerical evaluations of the method, reference data can be chosen that utilize MODIS data obtained at similar viewing zenith and relative-azimuthal angles to AHI [17].
In real applications of the transformation algorithm, finding pure endmember spectra from coarse spatial resolution (1 km) data is challenging [23]. Even though reflectance spectra used for endmembers were not pure, we expect that the NDVI transformation using Equation (12) functions properly if the FVC values computed using those endmembers are known to be equivalent between sensors or the correct ω value is known. Meanwhile, the effects of non-purity for endmember spectra obtained from 1 km resolution data on the transformation should be quantitatively evaluated using high spatial resolution (e.g., ≤10 m) data. This evaluation would, however, require considerable efforts for data processing and should be conducted in future works.
Note that the NDVI transformation for full-disk scale is not available using a single implementation of the method because atmospheric contamination should be as spatially uniform as possible and because variations of the sun and viewing angles should be as small as possible in the data for processing. Therefore, regional-scale NDVI transformation with moving-window approaches would be required to cover broader areas.
Although our proposed method still faces several unresolved problems, it can be used to predict near-real-time LEO-compatible GEO NDVI toward time-series vegetation analysis for a specific region. In this application, regional-or full-disk-scale reference sensor (LEO) reflectances with a certain temporal resolution (e.g., 8-day) should be prepared and used to predict the endmember spectra of a destination sensor for the time when a GEO system makes an observation. In addition, GEO reflectances observed in various solar zenith and azimuth angles should be normalized prior to or during the NDVI transformations.

Conclusions
In this study, we derived a GEO-to-LEO NDVI transformation equation based on a linear mixture model of anisotropic vegetation and nonvegetation endmember spectra. Using the derived equation, we developed an NDVI transformation method in which endmember spectra from each set of sensor data are computed independently and subsequently combined to compute the coefficients of the transformation equation. The proposed method could reduce the effects of surface anisotropy caused by viewing angle differences and SRF differences at the scene level. Moreover, the NDVI transformation is insensitive to calibration biases because the endmember spectra are computed using data. Notably, our method does not require regression analysis for the NDVI transformation because the endmember spectra used in the coefficients of the transformation equations are independently extracted from each sensor dataset.
Numerical demonstrations of the method were conducted using Himawari 8 AHI at off-nadir-view and Aqua MODIS at near-nadir-view data in middle latitude test sites, where the directional FVCs of the sites were reported to be nearly identical except for a few cases. The results showed that, compared with the original AHI NDVI, the transformed AHI NDVI was closer to the MODIS NDVI. The magnitudes of the averaged NDVI biases between the AHI and MODIS were reduced from 0.016-0.026 to <0.01 after the NDVI transformation was applied to the AHI NDVI. However, biases between FVCs from GEO and LEO sensors, which were caused by the viewing angle dependency of the directional FVC, may have added errors during the NDVI transformation. We expect that inclusions of the directional FVC biases would reduce the transformation errors.
The developed method is universally applicable to NDVI transformation among distinct sensors. Thus, the method can be used for not only GEO and LEO sensors but also, for example, LEO sensors on different platforms.
In future studies, further evaluations of the NDVI transformation using data obtained using various viewing angles and land-cover types are necessary. In addition, although this study used TOA reflectances for evaluation and demonstration purposes, the proposed method is capable of transforming partially or fully atmosphere-corrected NDVIs, which should be evaluated in a separate study. Furthermore, NDVI transformation across different atmospheric conditions (e.g., TOA to bottom-of-atmosphere) should be investigated and evaluated to explore the possibility of avoiding atmospheric corrections of original sensor data. Finally, the proposed method should be expanded to be capable of near-real-time GEO-to-LEO NDVI transformation across broader spatial scales. This would make it possible to obtain the spatial information of vegetation index time series with high temporal resolution and NDVI compatibility between GEO and LEO sensors.
Parameter o band(i,j) is a function of the second-and higher-order interaction terms between the surface and atmosphere, and o em,band(i,j) is the interaction terms between the surface of an endmember and the atmosphere. For the fully atmosphere-corrected reflectances, ρ a,band(i,j) and S band are zero and T ↓ band(i,j) and T ↑ band(i,j) are unity; thus, o band(i,j) and o em,band(i,j) become zero.