Correction of Directional Effects in VEGETATION NDVI Time-Series

: Land surface reﬂectance measurements from the VEGETATION program (SPOT-4, SPOT-5 and PROBA-V satellites) have led to the acquisition of consistent time-series of Normalized Difference Vegetation Index (NDVI) at a global scale. The wide imaging swath (>2000 km) of the family of VEGETATION space-borne sensors ensures a daily coverage of the Earth at the expense of a varying observation and illumination geometries between successive orbit overpasses for a given target located on the ground. Such angular variations infer saw-like patterns on time-series of reﬂectance and NDVI. The presence of directional effects is not a real issue provided that they can be properly removed, which supposes an appropriate BRDF (Bidirectional Reﬂectance Distribution Function) sampling as offered by the VEGETATION program. An anisotropy correction supports a better analysis of the temporal shapes and spatial patterns of land surface reﬂectance values and vegetation indices such as NDVI. Herein we describe a BRDF correction methodology, for the purpose of the Copernicus Global Land Service framework, which includes notably an adaptive data accumulation window and provides uncertainties associated with the NDVI computed with normalized reﬂectance. Assessing the general performance of the methodology in comparing time-series between normalized and directional NDVI reveals a signiﬁcant removal of the high-frequency noise due to directional effects. The proposed methodology is computationally efﬁcient to operate at a global scale to deliver a BRDF-corrected NDVI product based on long-term Time-Series of VEGETATION sensor and its follow-on with the Copernicus Sentinel-3 satellite constellation.


Introduction
The normalized difference vegetation index (NDVI) is the most widespread vegetation index [1,2]. It is defined as the difference-normalized by the sum-between the absorption of radiation in the RED spectral region, mainly caused by chlorophyll pigments, and the reflectance in the NIR (near-infrared) spectral region as a result of canopy structure. The computation of NDVI is as follows.
where ρ denotes the bidirectional reflectance factors (hereafter referred to as reflectances) in their respective spectral bands. Since soil and snow reflectance spectral values do not show a clear spectral difference between these bands [3], therefore NDVI was deemed relevant to highlight and monitor vegetated surfaces. Moreover, NDVI shows a significant relationship with fractional vegetation cover and leaf area index (LAI) [4]. Therefore, consistent long-term time-series of NDVI play a crucial role for the survey of the biosphere, in pointing out drought episodes [3] or anomalies in vegetation health [5]. The Copernicus Global Land Service ((CGLS) https://land.copernicus.eu/global/ index.html, accessed on 15 March 2021) is a component of the Land Monitoring Core Service that provides a series of bio-geophysical products about the status and the evolution of land surface for the whole Earth. The production and delivery made available in near real time conditions finds complement in long-term archive. The portfolio of the Copernicus Global Land Service which contains an NDVI product at 1 km spatial resolution. The so-called CGLS Collection 1 km NDVI Version 2 product 2 is a 10-daily maximum value composite (MVC [6] To ensure a temporal continuity in surface reflectance and NDVI time-series, the three VEGETATION instruments' designs were consistent in terms of radiometric definition and spectral bands characteristics and performance [7,8]. Since a prime requirement was to offer daily global coverage of land surfaces, the three platforms carrying the VEGETATION instruments flew on Sun-synchronous circular orbit at an altitude of ∼820 km to image the continental surfaces with a swath of ∼2200 km. The acquisition geometry for targets with such a wide swath can significantly differ between satellite overpasses. Moreover, the imaging Sun-sensor geometry for a given target changes according to the day of the year and step during the mission lifetime [9,10]. The high level of anisotropy of most land surface reflectance triggers spurious deviations in the NDVI time-series [11], which may significantly hamper further analysis. Minimizing the impact of directional effects in NDVI time-series is possible by considering satellite observations within a reduced range of viewing zenith angles (e.g., close to nadir). However, such a method discards clear sky observations that are mandatory to tune the time-series. The constraint implemented in the MVC composites of the CGLS Collection 1 km NDVI Version 2 product is aimed to reduce the impact of directional effects due to changing viewing conditions whereas the varying illumination is disregarded. Since the normalization of TOC time-series to a standard Sun-sensor geometry before NDVI computation improves the interpretation of seasonal variations of land surface vegetation derived from VEGETATION [12,13] and other sensors [14][15][16][17], there is an obvious need to remove directional effects from the Copernicus Global Land Service NDVI derived from VEGETATION instruments.
The primary goal of this paper is to describe a methodology to perform a Bidirectional Reflectance Distribution Function (BRDF) correction (removal of directional effects inferred by land surface anisotropy) to the NDVI time-series derived from the three VEGETATION sensors. The work is structured as follows: Section 2 briefly describes the datasets used in the methodology's processing and for performance evaluation. Section 3 details the different processing steps implemented in the methodology. The results of performance inherent to the BRDF correction are presented in Section 4 and further discussed in Section 5 while our conclusions are listed in Section 6. Satellite scenes images at 333 m and 100 m spatial resolution are acquired by the VGTPV instrument, albeit with less coverage for the second one (for a full description of the VGTPV instrument see [7]). The VEGETATION instrument is equipped of four wide spectral bands (BLUE, RED, NIR and SWIR) with only two considered here (see Table 1). VGT1 and VGT2 TOA reflectances used in this processor are the standard VGT-P Collection 3 [18] products provided by the SPOT/VEGETATION program through the VITO Product Distribution Facility (http://www.vito-eodata.be/ accessed on 15 March 2021). VGTPV TOA reflectances are the Collection 1 Level2A 1 km products issued from the PROBA-V nominal processing chain. The PROBA-V product user manual [19] (http: //proba-v.vgt.vito.be/sites/proba-v.vgt.vito.be/files/products_user_manual.pdf accessed on 15 March 2021) contains a full description of TOA reflectances measured by VGTPV.

VEGETATION Data Sets
A subdataset of the NDVI 1 km V2 product serves to illustrate the performance of the BRDF correction (see Section 4). As aforementioned, NDVI 1 km V2 product inputs are VGT1, VGT2 and VGTPV RED and NIR 1 km 10-days synthesis (using a Maximum Value Composite approach) of atmospherically corrected reflectance data sets. Please note that a spectral correction is applied only to the NDVI derived from VGT1. A full description of the Copernicus Global Land Service NDVI 1 km V2 product, hereafter referred to as NDVI MVC, can be found in the Algorithm Theoretical Basis Document (ATBD) [20] (https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/ GIOGL1_PUM_NDVI1km-V2.2_I2.31.pdf accessed on 15 March 2021).  Figure 1 shows the methodology's flowchart to obtain long-term series of BRDFcorrected NDVI from VEGETATION sensors radiometry. An atmospheric correction is first performed on RED and NIR TOA reflectances to obtain Top of Canopy (TOC) reflectances. Then, a spectral BRDF retrieval is performed on a pixel basis accumulating TOC values over a short time period. The retrieved BRDF models serve to set the RED and NIR VEGETATION TOC reflectances and NDVI values to a standard Sun-sensor geometry. A detailed description of the processing stages is shown in Figure 1 and detailed in the next subsections.

Atmospheric Correction
The VEGETATION TOA L2A orbital segments are cut out in tiles with a spatial extent of 10°. The atmospheric correction is performed in parallel for all tiles. The status map (SM) of the Level2A orbital segments is assessed before performing the atmospheric correction. Only land pixels flagged as clear and with good radiometric quality for the RED and NIR spectral channels are atmospherically corrected using the Simplified Model for Atmospheric Correction, SMAC [21]. The coefficients used in SMAC to perform the atmospheric correction in data from the VEGETATION sensors can be found in the CESBIO's repository (http://tully.ups-tlse.fr/olivier/smac-python/tree/master/COEFS, accessed on 15 March 2021). This algorithm is well-suited for operational application and is thus widely used in various projects, such as the Land Surface Analysis Satellite Application Facility (LSA-SAF) [22] and the VEGETATION operational ground segments [18,23]. Additional to the observed TOA reflectance, SMAC requires the angular information and aerosol optical thickness (AOT), ozone and water vapor columnar contents as input parameters.
The aerosol input consists of AOT monthly climatology data (applied to day 15 of each month) at a mid-visible (λ = 0.55 µm) wavelength. This climatology was generated from re-analysis data covering the period 2003 to 2012 at 0.125 • × 0.125 • resolution provided by the Copernicus Atmosphere Monitoring Service (CAMS), with bilinear interpolation applied to get monthly AOT data at 1 km resolution. The CAMS-based climatology is hereafter referred to as AOT CAMS .
A fall-back static AOT is used to replace the AOT CAMS when the SMAC calculations using AOT CAMS result in negative TOC reflectances. The static AOT is solely based on latitude and is hereafter referred to as AOT lat , which is defined as: Worth emphasizing that the AOT lat calculated from Equation (2) does not capture any longitudinal and temporal variation and has a limited dynamic range of 0.05 to 0.20. Although these limitations in AOT variation and magnitude are physically unrealistic, an asset of using AOT lat is before all to avoid negative TOC reflectances. For example, the static latitudinal AOT has been used to derive other surface reflectance datasets, such as in LSA-SAF, and the AOT fall-back mechanism in the PROBA-V Collection 1 atmospheric correction.
Ozone information is retrieved from monthly mean climatology data (applied to day 15 of each month) calculated from an 11-year Earth Probe Total Ozone Monitoring Satellite (EP-TOMS) dataset, available at 32 km resolution. A bilinear interpolation is applied to obtain the ozone information at one km resolution. Total column water vapor data are obtained from the European Centre for Medium-Range Weather Forecast (ECMWF) Numerical Weather Prediction (NWP) data. For each day, the column-integrated water vapor amounts are available at 6-h time steps (at 0, 6, 12, and 18 UTC) at a nominal resolution of 0.25 • × 0.25 • (∼27 km). Per segment, the nearest NWP time step is taken and subsequently a bilinear interpolation is performed spatially to obtain the water vapor columnar data at 1 km resolution. The SMAC TOC reflectance calculations include two steps, with the second one only to be performed if the first SMAC calculation using AOT CAMS delivers negative TOC reflectances. In such situations, the SMAC calculation is repeated, using this time AOT lat as input. This fall-back solution is necessary, as SMAC is somewhat limited in calculating accurate TOC reflectances for high AOT CAMS and large solar/viewing zenith angles.

BRDF Correction
The land surface reflectance varies with the observation and illumination geometry (defined by the sensor and Sun positions) and wavelength. The BRDF quantifies the reflectance anisotropy, which can be approximated by a semi-empirical BRDF linear model with three free parameters, commonly known as BRDF parameters or descriptors [24][25][26][27]. BRDF parameters are optimally retrieved based on an inversion method (see [28] for a review) over a period of time during which the surface properties are assumed not to undergo significant changes. The spectral BRDF models are then used to estimate normalized (to a standard Sun-sensor geometry) RED and NIR reflectances. These normalized reflectances are used to compute the BRDF-adjusted NDVI. The processing mentioned above, which is pixel-based, is summarized in a flowchart (see Figure 2). The processing of BRDF-adjusted NDVI correction module runs in three main stages (blue blocks in Figure 2) which can be resumed as follows: (i) First step consists of collecting a set of clear observations (output of the atmospheric correction module) which is achieved within the assembling dataset stage; (ii) Second step is achieved during the BRDF correction stage when the sets of observed RED and NIR TOC reflectances are fitted independently using a semi-empirical linear BRDF model [24] to obtain normalized reflectances; (iii) In a last step, the NDVI is computed. Each module of the BRDF correction is exemplified as blue blocks in Figure 2, which comprises individual processes (listed in Figure 2) that are described in the following subsections.

Accumulation
Normalizing surface reflectance time-series to a standard Sun-sensor geometry requires a knowledge of the surface BRDF. This is an inverse problem, and its solution involves the finding of wavelength-dependent kernels. The quality of the inversion procedure and further reliability of BRDF kernels strongly depends on the number and diversity of angular observations. Persistent cloud coverage can seriously jeopardize the collection of surface reflectance observations, thereby shrinking the ensemble of observations to invert. Gathering observations over a long period of time (e.g., >30 days) would provide enough data to constrain the inversion, but vegetation properties cannot be longer assumed stationary and temporal trend would contaminate the BRDF retrieval. On the other hand, if a too short composite period is retained, the risk exists to lack of sufficient observations to ensure a successful inversion.
Based on numerical experiments, an accumulation period of 16 days appeared as a fair trade-off between the objective to gather enough clear surface reflectance observations and to get the stationary properties of the surface. Please note that such strategy is consistent with MODIS MCD43A4 product [27,29] for which a period of 16 days was deemed appropriate to prevent smoothing out short timescale variations. Since the cloud/shadow mask product of VEGETATION could still be improved, accumulated clean observations may incorporate outliers resulting from some undetected clouds and shadows in the scene. Therefore, we performed a simple screening to remove these outliers using a modified z-score following the guidelines proposed by [30].

Compositing
To avoid introducing any delay on the timing of key phenological stages, we considered new observations as if they were acquired in the last ten days of the accumulation period. If there is sufficient new information, i.e., three new clear observations at least, the BRDF inversion uses only those. On the other hand, if there are not enough new data (less than 3) acquired during the last ten days, the entire dataset of observations accumulated over 16 days is used. This approach is what we called an adaptive window, implemented to ensure that NDVI key phenological stages are not too much garbled by data accumulation. If no clear observation is available in 16-days, the inversion is not achieved. After assembling the set of clear observations, the median date of observations used in the BRDF modelling is computed. The latter is provided as an independent layer (time-grid) and indicates the most representative date of the observations ensemble used to compute the BRDF-adjusted NDVI. The analyses of BRDF-adjusted NDVI time-series presented in Section 4 use the aforementioned time-grid.

Weights
By assigning weights, we favor certain observations through the inversion process. The weights assigned to the ensemble of observations have both temporal and angular components. The analytical expression is given below: It was proposed by [31] to estimate uncertainties associated with TOC reflectances as a function of solar zenith (θ s ) and viewing zenith (θ v ) angles. The coefficients c 1 and c 2 in Equation (3) are wavelength-dependent and are taken from Table 1 in [31]. Thus, the angular weights are the inverse of the TOC uncertainty assessment with Equation (3). These weights allow setting less importance to observations that were acquired at large viewing and solar zenith angles.
An effective temporal weight function is introduced by multiplying the prior covariance matrix (see Equation (10)) with a factor which is a function of the elapsed time since the last successful inversion. This temporal weight function has a decreasing exponential shape. It gives an enhanced weight to new observations. The timescale of such an exponential function is arbitrary set to 10 days, which means that a temporal weight obtained ten days before the last day of the accumulation period is 0.5. In contrast, the most recent observation has an effective temporal weight of 1 (maximum). This approach is consistent with the BRDF inversion methodology presented in [31,32]

BRDF Model
We use the Roujean semi-empirical BRDF model [24] which approximates the surface reflectance of a target by a linear combination of three kernels as in expression below: where θ s and θ v are the Sun and viewing zenith angles, respectively. The relative azimuth is the difference between viewing azimuth angle (φ v ) and Sun azimuth angle (φ s ). Functions f 0 , f 1 and f 2 are commonly known as the BRDF kernels. The f 0 is an isotropic function accounting for bidirectional reflectance with nadir viewing and the overhead Sun ( f 0 = 1); f 1 is a geometric function accounting for the effects of shadows and the geometrical structure of protuberances on the surface; f 2 is a function that simulates the volume scattering properties of the surface (for a comprehensive review of linear BRDF model see [28]). Full expressions for kernels f 1 and f 2 are presented in [24]. The geometric and scattering kernels are computed using the illumination and observing geometry associated with the assembled dataset. The wavelength-dependent kernel weights (k 0 , k 1 and k 2 ) from Equation (4) are then found via the BRDF inversion procedure, which runs independently for the RED and NIR spectral bands. The selection of the Roujean BRDF semi-empirical model is in line with the BRDF retrieval methodology of other Copernicus Global Land Service products [32].

BRDF Model Inversion and Adjustment
Given a set of n surface reflectance observations obtained under varying illumination conditions and viewing angles, the BRDF model can be inverted to further proceed to the directional correction. Its solution involves finding the wavelength-dependent kernel weights: k 0 , k 1 and k 2 (see Equation (4)). Moreover, the solution of the inverse BRDF problem allows us to obtain a model of the reflectance anisotropy of the surface; that can then be used to estimate reflectances under any observing and illumination conditions. Thus, based on an ensemble of n observations obtained from the compositing phase, the model serving for the inverse problem is expressed as: The expression above is an equation system which can be rewritten as follows: where F is an (n × 3) matrix containing the Roujean kernels for the ensemble of n observations. However, k is a (3 × 1) vector containing the three wavelength-dependent kernel weights and r is an (n × 1) vector with the observed reflectances. As described before, uncertainties associated to individual reflectance observations (angular component as defined in [31]) and their relevance within the compositing time window) are quantified by means of weighting factors (w j for j = 1, ...n). Thus, Equation (6) can be rewritten as: with a weighted reflectance vector defined as b with elements b j = r j w j and a design matrix A with elements A ji = F ji w j . The BRDF inversion is known for being an ill-posed problem e.g., [33][34][35], meaning that there is no unique set of k parameters that can solve Equation (7). However, significant variations of the surface BRDF shape parameters over a short timescale due to changes in surface cover are usually not expected. An approach to ensure a smooth evolution of the inverted k parameters is to make use of prior information obtained from previous BRDF inversions (see Figure 2). This approach has also been implemented in previous works e.g., [31][32][33]. Following Geiger et al. [31], the equation system to solve becomes: where k prior = [k 0 prior , k 1 prior , k 2 prior , ] is a vector containing the values of the k parameters found in the previous BRDF inversion and the diagonal matrix contains the uncertainties of the uncorrelated prior information on the k parameters. The previous estimates of the kernel weights (k prev ) and their uncertainties (C prev ) are used as follows as prior information: where t − t prev is the time difference (number of days) between the current production date and the date associated with k prev and C prev . As mentioned in Section 3.2.3, in consistency with Roujean et al. [32], the timescale τ is chosen as 10 days. So that the factor ∆ (defined as ∆ = 2 2/τ − 1) when applied to the prior covariance matrix could be equivalent to run the inversion using a temporal weight function with a decreasing exponential shape and timescale of 10 days [31,32]. This multiplicative factor applied to the prior covariance matrix reduces the confidence in the prior estimate. Moreover, it allows the BRDF inversion procedure to find the kernel weights that could best represent the current surface conditions captured by the observations, albeit being constrained by prior information on BRDF surface. Then, Equation (8) is solved as: where C = (A T A + C −1 prior ) −1 is the covariance matrix. Assuming uncorrelated errors, we consider the covariance matrix's diagonal elements as the variances associated with the wavelength-dependent k parameters obtained via Equation (11). In other words, It is important to notice that the quality of the BRDF inversion attested by the uncertainties associated with the retrieved BRDF parameters is strongly correlated with the number and diversity of observations used in the inversion.
Once we have approximated the surface anisotropy with a semi-empirical BRDF model [24], we use the wavelength-dependent models to estimate RED and NIR reflectances at a standard Sun-sensor geometry. Using the k parameters obtained from solving Equation (11), the BRDF-adjusted reflectances are calculated using Equation (4) for the solar zenith angle corresponding to 10:00 local time and viewing zenith angle at nadir.

NDVI Computation
In addition to the surface properties, the reflectances are also largely dependent on the specifications of the sensor spectral bands defined by the spectral response functions (SRFs). This means that there is not a true NDVI for a specific surface, but only an NDVI according to sensor-specific SRFs e.g., [36][37][38]. Consequently, the VGT1 NDVI is spectrally corrected to VGT2 NDVI (It corrects for the fact that the SRF of VGT1 in the RED band overlaps considerably with the transition zone from leaf absorption to foliage reflectance, which was adapted in the design of VGT2 and PROBA-V), by correcting the BRDF-normalized RED reflectance according to the relationship below: With a = 0.277052, b = −13.1103, c = 33.03465, d = −41.0406, and NDVI calculated from the uncorrected bands. The methodology used to derive these numerical coefficients is detailed in [20]. The VGTPV NDVI is not spectrally corrected to VGT2 NDVI, because both NDVI data sets have already a good matching that could not be improved by the estimated correction functions. Thus, the final step is the computation of NDVI by Equation (1).

Results
The primary purpose of this work is to describe the methodology developed to remove directional effects from the long-term VEGETATION NDVI time-series rather than present a validation assessment of the BRDF-adjusted NDVI product as this can be found in the Copernicus Global Land Service Quality Assessment report [39] (https://land.copernicus. eu/global/sites/cgls.vito.be/files/products/CGLOPS1_QAR_NDVI1km-V3_I1. 10.pdf, accessed on 15 March 2021). This section describes the results of the evaluation of the performance of the algorithm implementation. The methodology described in Section 3 was applied to a selection of land surface reflectances tiles (10 • × 10 • ) acquired with the VEGETATION series of sensors. A subdataset from the Copernicus NDVI 1 km V2 products (MVC NDVI), was prepared for the same study regions to assess the impact of the BRDF adjustment on the NDVI time series.
A notable example of directional effects on MVC NDVI time-series is a site located near Bukedea (latitude = 1.25 • , longitude = 34.08 • ), Uganda. Data for a pixel near the location above is extracted from the MVC NDVI product and the 2019 NDVI profile evolution along with its long-term average (computed from 1998 to 2020) are shown on the top panel of Figure 3. Despite the viewing angle constraints implemented in the synthesis of MVC NDVI product, directional effects on the 2019 NDVI time-series can be identified as saw-like patterns. These latter can seriously hamper the interpretation of observed NDVI temporal changes and derived information. For example, the computation of NDVI anomalies (comparison of current NDVI values to long-term statistics via a z-score) identifies whether vegetation growth conditions are ordinary. In this respect, the bottom panel of Figure 3 reveals the NDVI anomalies for a pixel located near Bukedea where bands denoting favorable (light green), very favorable (dark green), unfavorable (light red) and very unfavorable (dark red) conditions are also displayed.
From the bottom panel of Figure 3, it can be seen that the temporal evolution of NDVI anomaly follows abrupt changes with timescales shorter than the characteristic vegetation activity timescales. A specific example is the NDVI anomaly fluctuations between dekad 12 and 19. Within this short period which lasts about two months, the evolution of the NDVI anomaly shows a significant variability. The erratic switching of NDVI anomaly values, as seen in bottom panel of Figure 3, is related to the saw-like patterns in the NDVI time-series acquired with a changing viewing and illumination conditions. Thus, the Bukedea site illustrates the occurrence of directional effects in NDVI timeseries (see Figure 3) and allows us to verify whether our methodology can successfully minimize directional effects (short-term variability introduced by changing viewing and illumination conditions) in NDVI time-series. The NDVI profile and anomaly temporal evolution for the Bukedea pixel is shown in the top and bottom panels of Figure 4, respectively. At first glance, the saw-like pattern is not noticeable anymore in the BRDF-adjusted NDVI profile evolution ( top panel of Figure 4). The uncertainties provided in the BRDF-adjusted NDVI dataset allow us to identify less reliable normalized NDVI values. The later can occur during periods of persistent clouds or in the presence of outliers associated with misclassified clouds. The uncertainties on the NDVI are used as a weight in the computation of the long-term NDVI average, allowing down-weighting of the impact of less reliable NDVI values. The uncertainty associated with the current NDVI values is further propagated through the anomaly computations. Unlike, anomalies calculated from MVC NDVI (see the bottom panel in Figure 3), the NDVI anomaly evolution derived from the BRDF-adjusted NDVI dataset does not go through the anomaly regions rapidly. As a result, the bottom panel of Figure 4 shows an anomaly signal with variability timescales, and an overall trend that seems to be consistent with seasonal variations expected for a vegetated site. More details on the spatial and temporal consistency of the NDVI anomalies derived from the presented methodology can be found in Section 6 of its Copernicus Global Land quality assessment report [39]. To go further in the appraisal of the methodology, two additional sites (Belmanip 151, Belmanip 332) were considered to be they show marked NDVI seasonal patterns. MVC and BRDF-adjusted NDVI time-series for a pixel near Belmanip 151 site (Tanzania) are shown in Figure 5. The uncertainties associated with the BRDF-adjusted NDVI are not shown to avoid a cluttered plot. From Figure 5, NDVI time-series from both datasets show an overall agreement in terms of phenological variations with no noticeable lag among time-series. However, it can be noticed that during periods where high NDVI levels are observed, very rapid variations are noticeable on the Belmanip 151 MVC NDVI profile. The impact of directional effects on NDVI time-series is revealed by noise on the datasets [40], which leads to spurious deviations for a timescale less than the characteristic timescale of the vegetation activity. The rapid variability on NDVI time-series inferred by the BRDF effect is even more conspicuous at the end of the growing canopy periods [16] when NDVI reaches a peak. The lack of high-frequency variations in BRDF-adjusted NDVI Time-Series, especially in the summertime, demonstrates that directional effects have been properly removed.
MVC NDVI and BRDF-adjusted NDVI time-series for Belmanip 332 site (India) are shown in Figure 6 with the same line and color code as Figure 5. The NDVI time-series of Belmanip 332 site show a strong seasonal component which is characteristic of an area that is cultivated twice a year. From a visual inspection of Figure 6, it can be gleaned that the BRDF-adjusted NDVI dataset does not lag the directional NDVI time-series. The BRDFadjusted NDVI time-series, from Belmanip 151 and Belmanip 332 sites, follow the overall long-term phenological variations seen in MVC NDVI. From the examples above, it is noticeable that the onset of the growing season in BRDF-adjusted NDVI time-series is affected by the process of data accumulation. Since no significant delay on the NDVI time-series is being introduced by the BRDF adjustment, the agreement between both NDVI datasets, in regard to the temporal evolution, adds confidence on the adaptive window method that is proposed here (see Section 3.2.2). This latter acts for avoiding abrupt temporal NDVI variations and to overcome persistent gaps in time-series. However, an apparent difference between MVC NDVI and BRDF-adjusted NDVI in terms of NDVI peak levels can be noticed for specific growing periods of Belmanip 151 and Belmanip 332 sites. This apparent disagreement is discussed in Section 5.  Following previous works [13,16,[40][41][42][43], we quantified the impact of the BRDF correction on the VEGETATION NDVI time-series by computing an estimate of the statistical noise. According to Bréon and Vermote [41], the statistical noise in time-series can be defined as follows: where the NDVI values in a Time-Series are represented by ρ i . The first and last element are ρ i=0 and ρ i=n , respectively.
We use Equation (13) to estimate the statistical noise in MVC NDVI and BRDFadjusted NDVI Time-Series for pixels in three geographical subregions belonging to Eastern Africa, Western Europe and Northern India for a period of 5 years (January 2014-December 2018). These estimates are presented as noise maps in Figures 7-9. The noise in MVC NDVI and BRDF-adjusted NDVI Time-Series for a subregion located in East Africa is shown at the top and bottom panel of Figure 7, respectively. As shown in Figure 7, there is a noise reduction in the VEGETATION NDVI time-series after applying the BRDF correction. It is more demonstrated for vegetated areas than for pixels associated with barren soil areas (scaled noise < 3), for which the noise in MVC NDVI is already low. Figure 8 illustrates the artefacts for a subregion of Western Europe where the top and bottom panels display the statistical noise estimates for MVC NDVI and BRDF-adjusted NDVI time-series, respectively. The noise levels for MVC NDVI time-series display a dependence on the latitude. In particular, the MVC NDVI compositing strategy favors satellite observations within a reduced range of viewing zenith angles. Therefore, this compositing strategy-aimed to minimize directional effects in NDVI time-series-is significantly affected by the lack of clear observations. Since the frequency of cloud coverage is higher for northern latitudes, it is expected to observe noisier NDVI time-series as moving towards North. Top panel of Figure 8, revels this smooth trend with a residual effect in the BRDF-adjusted NDVI. Comparison between noise images displayed on top and bottom panels of Figure 9 (Northern India) shows at first glance an apparent reduction of noise for a region dominated by mountainous landscapes.

Discussion
Directional effects are inherent properties of time-series of land surface reflectance and vegetation indices derived from satellite observations. A comparison between directional and normalized NDVI time-series exemplifies the necessity for removing directional effects prior any analysis devoted to point out possible natural anomalies during vegetation growth, peak and decay. Two sites (Belmanip 151 and Belmanip 332) were selected to appraise the BRDF adjustment performance on the NDVI temporal patterns. Figures 5 and 6 show the NDVI time-series for Belmanip 151 and Belmanip 332, respectively. Despite an adequate agreement in terms of temporal evolution, a significant difference between MVC NDVI and BRDF-adjusted NDVI in terms of NDVI peak levels can be noticed. Since directional effects can significantly impact the levels and timing of maximum greenness [15], we explore whether the difference mentioned above is associated-or not-to a removal of directional effects.
For such, we performed a qualitative comparison between VEGETATION NDVI timeseries and BRDF-adjusted NDVI derived from the MODIS MCD43A41 product (https: //lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd43a4_v006, accessed on 15 March 2021) [27,29]. Nadir BRDF-adjusted reflectances (NBAR) for MODIS band1 (RED) and band2 (NIR) were retrieved for Belmanip 151 and Belmanip 332 sites via Google Earth Engine [44]. NDVI time-series derived from MCD43A4 NBAR data are displayed in Figures    Nevertheless, bottom panel of Figures 10 and 11 show a general agreement between BRDF-adjusted NDVI time-series. This confirms the efficiency of the directional correction and that the high peak levels of MVC NDVI (see top panels in Figures 10 and 11)-occurring mostly in periods of canopy development-are very likely linked to the anisotropy of the land surface reflectance. The consistency between both BRDF-adjusted NDVI Time-Series (MODIS and VEGETATION) relies on both the NDVI intensity and chronicle. Since the MODIS MCD43A product is provided with a daily cadence, the consistency of the NDVI temporal evolution seen in the bottom panels of Figures 10 and 11 adds confidence in the timing of BRDF-adjusted NDVI dataset key phenological stages achieved by the adaptive window approach proposed in our methodology.
To quantitatively evaluate the impact of the BRDF correction in NDVI time-series, it is of primary interest to perform a comparison in terms of statistical noise between MVC NDVI and BRDF-adjusted NDVI time-series. Herein, we define the noise reduction percentage as follows: where µ noise BRDF and µ noise MVC are the median of the noise distributions from MVC NDVI and BRDF-adjusted NDVI Time-Series, respectively. Figure 12 shows the normalized distribution of the noise estimates for all pixels (in the subregions selected) with a scaled noise estimate from MVC NDVI time-series larger than 3 i.e., noise×100 > 3). The latter allows us to remove pixels associated with barren soil areas for which the BRDF effect should not be considered significant, and a reduction of noise due to removal of directional effects is not relevant. The left panel in Figure 12 displays the normalized distribution of noise from MVC NDVI and BRDF-adjusted NDVI time-series color coded as indicated in its legend. From left-panel in Figure 12, the corresponding values of µ noise BRDF and µ noise MVC are 3.25 and 5.24, respectively. From Equation (14), the impact of the BRDF on the NDVI time-series for the subregion in East Africa can be approximated as a median noise reduction in VEGETATION NDVI time-series ∼38%. Although BRDF correction mitigates the effect of variations in solar and viewing zenith angles on surface reflectance, NDVI values are still affected by numerous factors, including sensor-specific and atmospheric conditions. The quality of the BRDF inversion is strongly correlated with the number and diversity of observations used in the inversion. There might be locations where, at periods, the gathering of clear observations is limited (or not possible), thereby yielding to an unsuccessful BRDF retrieval and thereby a corrected product with limited quality. Western Europe is a geographical region affected by periods of persistent clouds or natural darkness. Thus, noisy and scarce data lead to a noticeable latitudinal noise gradient in MVC NDVI (top panel of Figure 8). The lack of reliable observations translates into periods of less reliable BRDF-adjusted NDVI time-series. However, the latitudinal noise gradient seen in NDVI MVC is strongly quenched in the BRDF-adjusted NDVI dataset (see bottom panel of Figure 8). The noise distributions in MVC and BRDF-adjusted NDVI time-series for the selected subregion in West-Europe are shown in the middle panel of Figure 12. It is clear from this figure that the dispersion of the noise distribution in BRDF-adjusted NDVI Time-Series is significantly reduced when compared with MVC NDVI. Moreover, from the median values of the noise distribution, the noise reduction achieved with the directional correction is 33%, which is meaningful.
The right panel of Figure 12 shows the normalized distributions of noise in MVC NDVI and BRDF-adjusted NDVI time-series for the subregion in northern India; its median noise reduction is 28%. Although this noise reduction percentage is lower than the other two subregions, the difference in the shape of noise distributions (right panel in Figure 12) allows us to consider the noise reduction as significant. The latter is supported by a Kolmogorov-Smirnov test, resulting in a probability of both noise distributions being drawn from the same parent population <0.05. It should be noticed that the noise reduction on NDVI time-series achieved with the BRDF adjustment presented here does not depend on time (period of the time-series considered). As shown in Figure 13, the noise reduction (27%) also holds for the early period of the VEGETATION NDVI time-series (1999 to 2005) comprising VGT1 and VGT2 observations. Therefore, the noise reduction in BRDFadjusted NDVI datasets, as described above, is to verify that the NDVI BRDF adjustment methodology can significantly reduce the short-term variability in VEGETATION NDVI Time-Series acquired with changing viewing and illumination conditions.
The computation time to generate a tile (10 • × 10 • ) of the BRDF-adjusted NDVI product has been estimated by running the code on the MEP's cluster [45] in parallel mode (Spark). An average time to process a dekad is about 5 min per tile. In the VEGETATION tile grid, there are 336 tiles. Thus, the processing time required to generate all dekads of a time-series with 20-years is less than three days-using a single core per tile. The somewhat limited computation time required to implement our directional effects correction methodology to the full VEGETATION long-term archive is related to the linear kernel-driven models, which offer a trade-off in computational complexity and physically motivated parametrization of the surface reflectance anisotropy. The presented methodology is flexible enough to be upscaled to process data from another space-borne sensor with a wide imaging swath (allowing for adequate sampling of the BRDF) such as Sentinel-3.

Conclusions
An implementation of a methodology aimed to remove directional effects from the long-term VEGETATION NDVI Time-Series has been presented here with a detailed description of the core processing stages of BRDF correction and computation of a BRDFadjusted NDVI. In addition, VEGETATION data for subregions in Eastern Africa, Western

Conclusions
An implementation of a methodology aimed to remove directional effects from the long-term VEGETATION NDVI Time-Series has been presented here with a detailed description of the core processing stages of BRDF correction and computation of a BRDFadjusted NDVI. In addition, VEGETATION data for subregions in Eastern Africa, Western Europe and Northern India have been used to evaluate the expected BRDF correction performance on NDVI Time-Series. For each case, the impact of the BRDF correction on NDVI Time-Series is evaluated in taking the reduction of statistical noise as criterion. These examples show that an overall noise reduction of ∼ 30% in NDVI Time-Series is achieved after applying the BRDF correction. This estimate is consistent with previous works [40,43] that used the noise reduction in NDVI (or reflectance) time-series to assess the performance of the BRDF correction. Thus, the noise reduction in the processed BRDF-adjusted NDVI time-series suggests that the high-frequency noise component identified in MVC NDVI time-series is likely associated with directional effects. Moreover, it can be significantly diminished by the BRDF correction described here.
Since the methodology presented can successfully mitigate directional effects introduced by the very large swaths of the VEGETATION series of sensors and is computationally efficient, it was adapted to process the whole VEGETATION NDVI data archive (March 1998-June 2020) at a global scale. A full validation of this BRDF-adjusted NDVI dataset will be presented in a forthcoming article based on the Copernicus Global Land Service validation report of the methodology here presented [39]. The archival value of the long-term NDVI Time-Series dataset after BRDF correction allows a seamless extension of the ongoing vegetation activity monitoring by the Copernicus Sentinel 3 satellite constellation.