# Temporal patterns in illumination conditions and its effect on vegetation indices using Landsat on google earth engine

Pablo Martín-Ortega, Luis Gonzaga García-Montero, Nicole Sibelet

► To cite this version:

Pablo Martín-Ortega, Luis Gonzaga García-Montero, Nicole Sibelet. Temporal patterns in illumination conditions and its effect on vegetation indices using Landsat on google earth engine. Remote Sensing, 2020, 12 (2), pp.211. 10.3390/rs12020211. hal-05174630

HAL Id: hal-05174630

https://hal.science/hal-05174630v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Temporal Patterns in Illumination Conditions and Its Effect on Vegetation Indices Using Landsat on Google Earth Engine

Pablo Martín-Ortega<fnref n="1" />,<fnref n="2" />,<fnref n="3" />*, Luis Gonzaga García-Montero<fnref n="1" /> ID and Nicole Sibelet<fnref n="2" />,<fnref n="3" /> ID

Received: 29 November 2019; Accepted: 6 January 2020; Published: 8 January 2020

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** Vegetation indices (VI) describe vegetation structure and functioning but they are affected by illumination conditions (IC). Moreover, the fact that the effect of the IC on VI can be stronger than other biophysical or seasonal processes is under debate. Using Google Earth Engine and the latest Landsat Surface Reflectance level 1 data, we evaluated the temporal patterns of IC and two VI, the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) in a mountainous tropical forest during the years 1984–2017. We evaluated IC and VI at different times, their relationship with the topography and the correlations between them. We show that IC is useful for understanding the patterns of variation between VI and IC at the pixel level using Landsat sensors. Our findings confirmed a strong correlation between EVI and IC and less between NDVI and IC. We found a significant increase in IC, EVI, and NDVI throughout time due to an improvement in the position of all Landsat sensors. Our results reinforce the need to consider IC to interpret VI over long periods using Landsat data in order to increase the precision of monitoring VI in irregular topography.

**Keywords:** illumination condition; Landsat; Google Earth Engine; EVI; NDVI; topographic effects

## 1. Introduction

Vegetation indices (VI), defined as the arithmetic combination of two or more bands related to the spectral characteristics of vegetation [1,2], have been used in a variety of fields including phenology, classification of vegetation, photosynthetic activity, aboveground net primary productivity and land surface temperature [3,4]. VI, particularly the Normalized Difference Vegetation Index (NDVI), are essential components of any study aiming to investigate ecosystem services especially those where vegetation, water, and biodiversity are involved [5,6].

However, VI sensitivity is affected by the changing radiance that accompanies changes in orientation of the vegetation surface being sensed [7]. The radiance changes at different times in the year and between years, due to different solar incidences over the surface, the so-called sun-sensor geometry [8]. Radiance is further changed in rough terrain, where a combination of the orientation of the terrain and the position of the satellite will determine high or low illumination conditions (IC) [7–10]. The ratio properties of NDVI enable a large proportion of noise caused by changing sun angles, topography, clouds or shadow to be cancelled, making this index less susceptible to IC. However, this index is still affected by atmospheric, soil and background vegetation cover feedback conditions. Furthermore, in tropical forests with high levels of biomass or high Leaf Area Index (LAI) values, NDVI also tends to saturate [7,11]. The Enhanced Vegetation Index (EVI), in turn, is more

<footnote n="1">¹ Department of Forest and Environmental Engineering and Management, Universidad Politécnica de Madrid UPM, Ciudad Universitaria s/n, 28040 Madrid, Spain; luisgonzaga.garcia@upm.es</footnote>

<footnote n="2">² Innovation, University Montpellier, CIRAD, INRAE, Montpellier SupAgro, 34000 Montpellier, France; nicole.sibelet@cirad.fr</footnote>

<footnote n="3">³ CIRAD, UMR Innovation, F-34398 Montpellier, France</footnote>

<footnote>* Correspondence: pablo.martin@geodata.consulting</footnote>
sensitive than NDVI to biophysical attributes such as LAI [11,12] and much more affected by IC than
NDVI, because is not a ratio-based VI and cannot compensate for variations in IC. In addition, EVI was
proven to be five times more sensitive than NDVI to changes in near-infrared reflectance (NIR) [11–14].

Due to the effect of IC, the use of EVI as an indicator of vegetation functioning or forest productivity has been put under debate in recent years. Some authors claimed that the unusual greening effect observed in the dry season in the Amazon forest using the EVI Moderate-resolution Imaging Spectroradiometer (MODIS) imagery [15] was induced by changes in sun-sensor geometry and not as a result of canopy structure, phenological patterns, or vegetation functioning [16]. This effect has been confirmed in similar ecosystems at different times of the same season, with different IC [14,17]. However, in similar tropical forests, after the removal of the IC effects, seasonal patterns were still present and seemed to be correlated with gross primary production (GPP), although authors recommended being cautious with this correlation [13]. Effects related to sun-sensor geometry and topography were also found at different times of the year using EVI and NDVI in subtropical deciduous forests under different IC. Sunlit and shadowed surfaces showed, respectively, different intensities of decrease and increase in reflectance even after topographic correction [12].

Most of the studies referenced above used coarse spatial resolutions (MODIS) and/or limited dates to account for seasonal or intra-annual variability using VI [7,12,14,16]. When tried to correct IC in Landsat sensors, for either global coverages or long time-series only some of them covering limited time spans were included in the analysis [18–20]. In general, IC effects were calculated using the information on solar zenith and azimuth angles for MODIS available imagery or using this sensor but with information taken from Landsat [18,21,22]. Transforming image angle information into a model of IC using data on digital elevation models (DEM) provides higher detail at the pixel level, especially in irregular terrain [10,23] and is important to compare different spatial resolutions [24,25].

The monitoring of IC and VI spatiotemporal patterns can be improved using higher spatial resolutions, longer time series, and the best-calibrated data to minimize errors in the analysis. However, the computation of large datasets of imagery, preprocessing, and downloading can be time-consuming and limited for traditional processing methods. Platforms like Google Earth Engine (GEE) provide free access to high-performance computing resources allowing to produce a vast amount of computation analysis to unprecedented time and extent scales [26]. Regarding quality in terms of spatial and temporal resolution, Landsat collection since the early 1970s [27] became freely available which widely expanded its use among researchers [28]. In May 2018 in a new effort, the United States Geological Survey (USGS) reorganized the archive into a tiered collection launching the Landsat level 1 product. This product provides a consistent quality of Landsat images that are comparable across Landsat sensors 4, 5, 7 and 8, suitable for time-series analysis and therefore ideal for monitoring [29].

In this study, we used Landsat level 1 surface reflectance data to evaluate the spatiotemporal
patterns of IC, EVI, and NDVI over a tropical forest in Costa Rica in irregular terrain. We tested GEE
capacity to calculate IC over a large number of images between the years 1984–2017. The specific
objectives of this study were: (1) evaluate the patterns of IC change over time in Landsat and the factors
behind these patterns, (2) explore the relationships between IC, EVI, and NDVI in changing conditions
of topography and time with improved resolution, (3) assess the importance of IC when monitoring
forested areas in irregular terrain. Our study supports previous evidence about the importance of
using IC when monitoring forests and gives additional insights about the effect of the topography in
Landsat sensors.

**2. Materials and Methods**

## 2.1. *Study Area*

The study area is a private property located in the Central Valley in Costa Rica within the province of Cartago (9°46′31.17″N, 83°45′27.31″W, WGS 84). The vegetation is classified as an evergreen premontane rain forest bordering with Tapanti-Cerro de la Muerte, an important national conservation
area (Figure 1). The altitude ranges between 986–1302 m.a.s.l. The property is devoted to conservation and ecotourism with 88% of the area covered by little disturbed native forest. We visited the property in August 2017. After screening the area for forest changes using aerial photo single frames available for free at https://earthexplorer.usgs.gov/ and high-resolution imagery at different years available at Google Earth pro (Figure 2) the area was considered as an invariant forest spot to carry out our long-term analysis.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Geographical location of the area of interest (AOI), the red polygon corresponds to the Landsat tile for path and row 15/53.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Area of Interest (AOI) (a) Aerial single frame of the area 11 August 1992, (b) Very high-resolution satellite image 16 January 2017 (DigitalGlobe), (c) Aspect (°), (d) Slope (%).</figcaption></figure>
## 2.2. Landsat Datasets and Image Processing

Data were available and processed with the GEE platform. Landsat Surface Reflectance Tier 1 collections for Landsat 4 Enhanced Thematic Mapper (ETM), Landsat 5 ETM, Landsat 7 ETM+ and Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) were used in the analysis. Calibration between sensors has been done, so they are suitable for time-series analysis. Correspondent ImageCollection IDs in GEE are: 'LANDSAT/LT04/C01/T1_SR'; 'LANDSAT/LT05/C01/T1_SR'; 'LANDSAT/LE07/C01/T1_SR' and 'LANDSAT/LC08/C01/T1_SR'. The collections correct illumination/viewing geometry and atmospheric effects using LEDAPS in ETM sensors and LaSRC in OLI/TIRS sensors. [30,31]. We filtered the collection for all images available during the period 1 January 1984 to 31 December 2017. Path and row for all images were 15/53.

For the removal of clouds, we used the C Function of Mask CFMask algorithm [32]. After the removal of cloudy pixels in the area we obtain clear and masked cloudy pixels for every image. We used GEE to count the highest number of contiguous clear pixels inside our AOI, for the highest number of available images. This reduced the initial area to a limited number of pixels (n = 1228) (Figure 3). Finally, these pixels were cloud-free for 39 images (Figure 4). We inspected every individual image visually to search for cloud or haze remnants.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Landsat pixels available (30 × 30 m), n = 1228 for the 39 selected images in the AOI after cloud masking, and screening.</figcaption></figure>

Two common VI were calculated for every image NDVI [2] (Equation (1)) and EVI [1] (Equation (2)) and added as new bands.

$$NDVI = (NIR - RED)/(NIR + RED) \qquad (1)$$

$$EVI = G * ((NIR - RED)/(NIR + C_1 * RED - C_2 * BLUE + L)) \qquad (2)$$

where NIR, RED and BLUE are Near-infrared; red and blue bands for Landsat images. G (2.5) is a scaling factor; L (1.0) is the canopy background adjustment factor; and C₁ (6.0) and C₂ (7.5) are the coefficients of the aerosol resistance term [33].
<figure><img src="image_7.png" /><figcaption>Figure 4. Dates expressed as year and day of year when the 39 selected images were available for the analysis.</figcaption></figure>

## 2.3. Illumination Condition (IC)

The incidence angle (i), is the angle between the normal to the pixel surface and the solar Zenith angle (z) [34,35] (Figure 5). The radiance detected by the sensor after interacting with surfaces of different slopes and aspects is dependent on this incidence angle (i). If (i) is low, the directions of the sun and the sensor are closer and radiance will be higher, the scene will be well illuminated and shadows due to topography minimized. Conversely, if (i) is large, the light will reach the surface in an oblique direction, illumination will decrease in the scene and more shadows will be cast due to the topography effect. Illumination condition (IC) is the cosine of the incidence angle (i) and ranges from 0 to 1, respectively, with 0 being values for poorly illuminated and 1 for well-illuminated areas. IC is calculated in Equation (3).

$$ \cos i = \text{Illumination Condition (IC)} = \cos z \times \cos s + \sin z \times \sin s \times \cos (a - o) \quad (3) $$

where z is the solar zenith angle; s is the terrain slope; a is the solar azimuth angle and o is the terrain aspect [8]. To calculate slope and aspect we used a DEM, the Shuttle Radar Topography Mission (SRTM) which is also available as a dataset in GEE (ID: 'USGS/SRTMGL1_003') [36]. Zenith and azimuth angles are available in the metadata of each Landsat image. A script was created in GEE using SRTM and metadata solar angles in the image to calculate IC at the pixel level and incorporate it in each image as a new band.

<figure><img src="image_8.png" /><figcaption>Figure 5. Terrain and solar angles involved in the calculation of the incidence angle (i) and illumination condition (IC).</figcaption></figure>
## 2.4. Statistics Across the Collection of Images

The 39 selected images, each one with the three bands IC, EVI, and NDVI (Figure 6b) were organized as a GEE Image Collection object. We used “ reducers ” which are a powerful tool in GEE that allow us to make computations at the pixel level across a stack of images or limited to a region in space (Figure 6a). We used the ee.Reducer algorithm to calculate descriptive statistics such as mean and standard deviation and correlations such as Pearson correlation. Mean IC and standard deviation IC were calculated at the pixel level for the 39 images. Mean IC shows which pixels are more frequently shadowed or illuminated with respect to the terrain and the sensor. The standard deviation image indicates the magnitude of variation of the IC for each pixel. Finally, we obtained two images, a mean IC image, and a standard deviation IC image.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Example of computation of descriptive statistics using the ee.Reducer algorithm across the stack of 39 images (mean in the example) (a). Each image has three bands IC, Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) (b).</figcaption></figure>

Because the amount of reflected radiation from each pixel is dependent on IC and this influences VI [7,37], we analyzed the correlation between EVI~IC and NDVI~IC using Pearson correlation to check if changes in IC are related to changes in VI [38]. Pearson correlation has been shown to be a good descriptor of the linear relationship between IC and VI in previous studies [17,37,39]. The correlation is also computed using the ee.Reducer algorithm and the 39 pairs of values present at each pixel using IC as X and VI as Y (EVI or NDVI). The result for each VI~IC correlation is two images, one containing the values of the Pearson correlation coefficients and the other containing their significance values (p) at the pixel level.

We also calculated the mean IC, EVI and NDVI values for all pixels available at each individual image of the selected images (n = 39). This allowed us to compare mean values between images at different times. Finally, because IC can be modeled for all the images available (n = 397), we calculated the mean IC for all of them to establish the temporal trend. Results were exported to CSV files and analyzed using either R software or EXCEL. Raster images and maps were generated using ArcGIS 10.6.1.

# 3. Results

## 3.1. Illumination Condition and Vegetation Indices

Mean and standard deviation of IC values for the selected images were calculated for each pixel and transformed into an image (Figures 6 and 7). Mean IC values ranged from 0.47 to 0.95, whereas the standard deviation ranged from 0.027 to 0.208. Values for both images were unevenly distributed across the image as a result of the effect of topography and sun-sensor geometry. Pixels with low mean IC values correspond to frequently shadowed areas across the 39 images (dark pixels in Figure 7a). Due to the sun-sensor geometry, these areas are placed obliquely or even opposite to the sun direction and the radiance received by the sensor is lower or null. Pixels with high mean IC values (bright
pixels in Figure 7a) are better placed in the sun direction and the radiance received by the sensor is higher. Using the standard deviation as a measure of the magnitude of change in IC for each pixel, we observed that pixels with high mean IC values experience little variation (dark pixels in Figure 7b), compared to those with low mean IC values which experienced greater variation (bright pixels in Figure 7b). The visual effect looking at the images is that, overall, both mean and standard deviation IC images show opposite patterns, except for an area with low-medium mean IC values and low-medium IC variation (green ellipse in Figure 7). This area would correspond to the terrain where the surface is parallel to sunlight beams (Figure 8), therefore, keeping low mean IC values and low variation. Figure 8 shows how mean IC is related to aspect and slope in the terrain. Aspect shows the highest mean IC values in the interval of 80–160° with maximum values around 120°. The lowest mean IC values are in the interval 320–360° and 0–40°, with the minimum around 360° (Figure 8a). Areas with aspect values around 120° would correspond to areas well oriented to direct sunlight, whereas those around 360° would be placed opposite to the sun (shadowed). We hypothesize that the areas described as parallel to sunlight beams might be those around 200–280°, where low mean and standard deviation IC values are found. These areas can be spatially located in Figure 2. The slope shows high mean IC values (~0.8) for flat areas (slope = 0°), and then increasing slope values show a variety of IC depending on how the surface is oriented to the sun (Figure 8b). This effect can be seen easily between the highest mean IC value (0.95) found at 34° of slope and 119° of aspect, whereas the lowest mean IC value (0.47) is found at 28° of slope and 321° of aspect, both slope values are close, but they have opposite orientations (202° of difference). Although variations in the IC have proven to be normal as an effect of seasonality and sun-sensor geometry variability across the year, we would expect to find IC variations with a similar degree of change across the whole image. These results suggest that variations depend on terrain conditions (Figure 8). However, if we consider that terrain, topography and cover are constant, the uneven variation in different IC values must be affected by changes in sun-sensor geometry. To investigate the magnitude of these changes and their influence on VI, we computed the correlations between VI~IC.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Spatial distribution of the mean illumination condition (a) and its standard deviation in the area (b). Each pixel shows the result of the 39 selected images. Red and blue circles show areas with opposite patterns. The green ellipses show areas with low mean IC and variation.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Scatter plot showing mean IC and aspect (°) in (a) and IC and slope (°) in (b) at the pixel level.</figcaption></figure>

The analysis of the Pearson correlation between IC and VI on a pixel-basis showed both negative and positive trends in the area clearly associated with IC. The correlation EVI~IC showed a higher range of values (−0.55–0.97) than the correlation NDVI~IC (−0.25–0.73), confirming that EVI is more sensitive than NDVI to changes in IC. Positive correlation values were more abundant than negative ones, describing overall positive trends for both VI, meaning that for an observed increase in IC, VI also increases. However, not all correlations were significant. The EVI~IC correlation was significant (*p* < 0.05) with positive correlations for most of the pixels and negative correlations for a few pixels (Figures 9 and 10). Areas with positive EVI~IC correlations correspond to pixels with low mean IC in Figure 7a, suggesting that EVI is very sensitive to increase when IC increases in poorly illuminated areas, showing the strongest correlation in pixels with the lowest IC (Figure 8a,b). A small number of pixels showed significant values for negative correlations between EVI~IC in well-illuminated areas. This can be related to the appearance of tree shadows between canopy levels when IC conditions change, which has been found in other studies [40]. The response of EVI to trees shadowed by other trees would be lower. However, higher resolution imagery would be needed to describe the canopy structure and the distribution of canopy shadows with different IC to account for this effect. The correlation NDVI~IC was significant only for some areas that showed positive correlations and not in areas with negative or no correlations (Figures 9 and 10). The strength of the positive correlations NDVI~IC was also lower compared to those in EVI~IC. The increase of NDVI with increased IC was also found in pixels with low mean IC and low IC variation, previously described as having a parallel orientation to that of the sunlight and limited to a well-defined aspect range. It is known that NDVI saturates rapidly with high biomass levels and dense canopies, and is less sensitive to changes in IC, which would explain the low correlations NDVI~IC found [7,11]. However, the fact that NDVI seems to be sensitive to changes in areas with low mean IC and low IC variation is a new finding that needs to be further investigated. Because the correlation is calculated with all IC and VI values available for each pixel and these pixels correspond to different dates in different seasons and years, the correlations for both VI are reflecting intrinsic temporal patterns associated with IC.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Spatial distribution of the Pearson correlation values between EVI~IC (a); Pearson correlation values between NDVI~IC (b); significance values (*p* < 0.05) for Pearson correlation between EVI~IC (c); and significance values (*p* < 0.05) for Pearson correlation between NDVI~IC (d). Each pixel shows the result of the 39 selected images.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Scatter plot showing p values and Pearson correlation coefficient values for the correlation EVI~IC (a) and NDVI~IC (b). The black vertical line represents *p* = 0.05, being all points to the left significant at *p* < 0.05. Colors represent low mean IC values (purple) to high mean IC values (yellow).</figcaption></figure>

When the relationship between EVI~IC and NDVI~IC is described comparing the mean values for each selected image (n = 39), EVI and IC showed a significant and positive dependence (R² = 0.5715, *p* < 0.01) (Figure 11), whereas NDVI and IC showed a very small to almost no dependence (R² = 0.075, *p* < 0.09) (Figure 12). While these results corroborate the higher sensitivity of EVI and the lower sensitivity of NDVI to the effect of IC, they are not able to explain alone the VI~IC variations found when we examine the image at the pixel level.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Scatter plot of the mean illumination condition (IC) versus the mean EVI value for each of the 39 selected images.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Scatter plot of the mean illumination condition (IC) versus the mean NDVI value for each of the 39 selected images.</figcaption></figure>

Furthermore, the images used in this analysis correspond to different dates in different years. Because it is known from other studies that IC change at different dates [13,14,16,17,39], the images were chronologically ordered to investigate temporal trends on IC in Section 3.2.

## 3.2. Time Series for IC, Enhanced Vegetation Index (EVI), and Normalized Difference Vegetation Index (NDVI) from 1984–2017

Temporal trends for IC, EVI, and NDVI show positive and similar slopes and are highly significant ($p < 0.01$) (Figure 13). NDVI shows the best adjustment ($R^2 = 0.57, p < 0.01$) followed by EVI ($R^2 = 0.54, p < 0.01$) and then by IC ($R^2 = 0.28, p < 0.01$), the latter which showed the poorest adjustment due to the high variation as a result of known seasonal changes. For the 39 selected images, the minimum mean IC value (0.71) corresponds to the date 31 December 1989 with mean EVI and NDVI values of 0.48 and 0.84 respectively. The maximum mean IC value (0.88), corresponds to the date 4 September 2016 with mean EVI and NDVI values of 0.59 and 0.88 respectively. A comparison between two close images, both from the dry season but separated 28 years also depicts the effect of the increase of mean IC in both mean EVI and NDVI (Figure 14).

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Time-series plot of the mean IC, mean EVI and mean NDVI values for all the 39 selected images. The solid black line represents the mean IC calculated for all the images available ($n = 397$).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Barplot comparing mean values of IC, EVI, and NDVI for two close images in date but 28 years apart.</figcaption></figure>

Observing an increase in IC, EVI, and NDVI over time despite having a limited number of images but still representative of different seasons and years to make this analysis would also indicate that there is a general improvement in IC and hence in EVI and NDVI from old to new Landsat images.

The modeling of IC for all images in the period studied is shown in Figure 13. The mean IC value in the area shows an overall increase over time. This can be clearly observed comparing images close in date but separated in years, where the increase in IC occurs to a similar degree for every date and month compared (Figure 15). In more detail, the mean IC trend experiences an overall increase with a slight drop between the years 1989–1993 and then an increase again. It can also be observed that due to missing images in some dates between the years 1984–2012, the IC seasonality pattern cannot be as clearly seen as after the year 2012 to present when the coverage of images is more complete. Overall, EVI and NDVI show higher values in the recent period 2014–2017 than in the previous period 1988–2002 regardless of the seasonal effects. The lack of clear images between the period 2003–2013 is due partly to cloudiness and to the Scan Line Corrector (SLC) failure in Landsat 7 [29].

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Comparison of mean IC for images close in date but separated in years. The criteria of selection were first close dates, then the highest number of years in difference.</figcaption></figure>

Monthly trends are shown in Figure 16. Clear seasonal trends can be observed. The lowest mean IC occurs in the dry season (December–February) and higher IC can take place in two different times of the year, at the beginning of the rainy season (April–May) and in the short dry period within the rainy season (August–September). Every month shows a similar and steady increase in IC over time from old to recently acquired images, although the strength of the increase varies between months. March shows the best adjustment in the increase of IC with time (R² = 0.71), whereas November showed the
lowest ($R^2$ = 0.18). Besides the drops in IC between 1989–1993 that we commented on previously and that can also be seen in the figure, anomalous changes are occurring in the month of December for the latest imagery available. However, this behavior would require further study. The magnitude of the increase reflects constant trends across years for all months. Although beyond the scope of this study, further research would be needed to explore these trends with phenology patterns.

Figure 17 shows the mean IC for all Landsat images in the area in chronological order indicating the sensor capturing the scene. Most of Landsat 4 and 5 images covered the period 1987–2001 and then Landsat 5 alone covered some dates between 2007–2010. The period for Landsat 7 ranges from 1999–2017 with a period of synchronicity with Landsat 8 between 2013–2017 that has produced more frequently available images during this time. The lowest mean IC value (0.69) was found on 31 December 1992 using Landsat 4 and the highest (0.89) on 15 September 2017 using Landsat 7. More detailed analysis also shows improved IC within the same sensor. For a variety of close dates in different years, Landsat 5 and Landsat 7 show a general increase in IC.

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Temporal trend of mean IC in the study area by month and year for all Landsat images selected (n = 397). Blue lines represent regression lines and grey areas confidence interval (CI) at 95%. The number of images available for each month(n) and $R^2$ are displayed.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Temporal trend of IC in the study area for all Landsat images selected (n = 397) and all Landsat sensors.</figcaption></figure>
# 4. Discussion

## 4.1. Illumination Conditions and Vegetation Indices

The IC model provided a detailed view of the effect of sun-sensor geometry and topography at a pixel level. Our research reveals that frequently shadowed areas experience more variation in IC than sunlit areas. The increase in IC values with the improvement of sun-sensor geometry has been reported in other studies, but usually comparing two images at different times or several images across a season in the same year [7,12,41]. The effect of terrain in the correlations between EVI~IC and NDVI~IC has important implications when interpreting VI, especially EVI, which showed higher sensitivity. Because EVI~IC correlations were positive in shadowed areas and neutral to negative in sunlit areas, it can be interpreted that the areas with different IC experienced "greening" or "browning" when in reality there was no change in the vegetation conditions. The explanation for this would be that sunlit areas keep high IC values with low variation due to more constant sun-sensor geometry conditions. This would have a low effect on the variation of EVI and NDVI which correspond to neutral or even negative correlation values. However, shadowed areas experience a higher range of IC values which would have a higher effect on EVI and NDVI values, showing positive correlations. The strong effect observed in EVI in this irregular terrain corroborates the need for removing the topographic effect in the reflectance data before calculating this index [7]. Therefore, the results using EVI for modeling phenology patterns, vegetation functioning or GPP would be biased if the IC effect is omitted, in accordance with research supporting this effect [3,17,40,42]. In relation to the research that claims the effect of IC on VI, most of these studies used coarser resolution than Landsat and did not employ a DEM to simulate IC. In addition, most of them were carried out in the Amazon region in relatively flat terrain and yet IC effects were patent [12,14,16,17]. We agree with previous studies that found drastic forest changes comparing various images in irregular terrain [10] and we suggest that using IC is critical to understanding these patterns. Although using a DEM with the same scale as the satellite data improved the results of our interpretation, caution is required because of known issues in SRTM data such as shadows that can bias the results, further investigation is needed [23]. Despite this, the spatial resolution of the images used, allowed more detailed and better interpretations of changes in IC than sensors with coarser spatial resolution such as MODIS or methodologies relying upon angle information of the image [12,14,16,17].

VI are also sensitive to vegetation structure factors such as LAI or canopy shadow which very likely have an effect on the variations observed, but lower resolutions than those used in this study would be required, especially for EVI [40]. Factors such as the background influence of soil and leaf litter or LAI are probably responsible for quickly saturating NDVI [43], hiding the effect of IC variation on this index.

## 4.2. Temporal Analysis of Illumination Conditions

Temporal analysis for the whole Landsat reflectance collection in the area showed significant and increasing trends for IC, EVI, and NDVI. We employed information on zenith and azimuth solar angles already present in each available image, eliminating the need for merging information between Landsat and MODIS sensors to model or simulate data. Intra- and inter-annual variations of VI as a response to changes in sun-sensor geometry and variations in sun-sensor geometry itself have been conducted using MODIS and Landsat collections, but covering shorter periods of time [18,19,21]. Some of the observed changes in IC agree with observations reported after analyzing the orbit change in Landsat 5 between the years 1995–2000, but here we describe the IC variations better at the terrain level, something that was not possible to reflect in that study [18]. Although the terrain is believed to remain constant, the use of the SRTM dataset to model IC throughout the whole period could be a source of bias in some areas.
## 5. Conclusions

In this study, we have developed a novel approach to using Landsat temporal series to evaluate changes in the illumination conditions and vegetation indices in forested areas in irregular terrain. We used Google Earth Engine to analyze 397 images for 28 years (1984–2017) of the latest Landsat Surface Reflectance collections. The main conclusions of our work can be described as follows:

(1) In summary, the illumination condition of Landsat scenes can be easily calculated using Landsat image information and an elevation model of the same resolution. Changes in illumination condition and its relation to vegetation indices can be evaluated effectively at the pixel level without needing to use ancillary data from other sensors or simulating data.

(2) The analysis at the pixel level showed correlations between illumination conditions and vegetation indices with a strong effect in EVI and less in NDVI. Moreover, positive correlations were found in shadowed areas, whereas neutral to negative correlations were found in sunlit areas. This, is a crucial finding to interpret vegetation indices derived from Landsat data in irregular terrain. In this regard, our analysis at the pixel level revealed patterns that could not be detected using only solar angle information.

(3) Our long-term analysis revealed that there is an increasing trend in illumination conditions, EVI, and NDVI associated with the improvement in the position in Landsat sensors over time. These trends calculated at selected images at varying seasons were significant when placed in chronological order. This effect cannot be overlooked when employing vegetation indices in the study of phenology patterns or aboveground net primary productivity. The incorporation of illumination condition into time-series analysis according to this method can provide additional data to understand the behavior of vegetation indices in irregular terrain.

**Author Contributions:** L.G.-M. and N.S. conceptualized the research ideas and supervised the research activity. P.M.-O., designed the methodology, analyzed the data and wrote the paper. L.G.-M. and N.S. revised the manuscript and provided useful suggestions on the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been conducted as part of a Ph.D. thesis project supported by the Agricultural Transformation by Innovation (AGTRAIN) Erasmus Mundus Joint Doctorate Programme, EU Framework Partnership Agreement 2011-1530/001-001-EMJD (2011.0019)—Contract AGTRAIN NO. 2015-005 funded by the EACEA (Education, Audiovisual and Culture Executive Agency) of the European Commission

**Acknowledgments:** The authors would like to thank the anonymous reviewers in their effort of reading and giving constructive feedback on the manuscript. We also acknowledge the institutional support of the Department of Forest and Environmental Engineering and Management of Universidad Politécnica de Madrid, Montpellier SupAgro, CIRAD, UMR Innovation and Centro Agronómico Tropical de Investigación y Enseñanza (CATIE). We would also like to thank feedback from the Google Earth Engine team on the code and from anonymous users in the GEE forum.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Liu, H.Q.; Huete, A. Feedback based modification of the NDVI to minimize canopy background and atmospheric noise. *IEEE Trans. Geosci. Remote Sens.* **1995**, *33*, 457–465. [CrossRef]

2. Rouse, J.W.; Hass, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation systems in the great plains with ERTS. *Third Earth Resour. Technol. Satell. Symp.* **1973**, *1*, 309–317.

3. Cao, S.; Sánchez-Azofeifa, G.; Duran, S.; Calvo-Rodriguez, S. Estimation of aboveground net primary productivity in secondary tropical dry forests using the Carnegie–Ames–Stanford approach (CASA) model. *Environ. Res. Lett.* **2016**, *11*. [CrossRef]

4. Liu, K.; Wang, S.; Li, X.; Li, Y.; Zhang, B.; Zhai, R. The assessment of different vegetation indices for spatial disaggregating of thermal imagery over the humid agricultural region. *Int. J. Remote Sens.* **2020**, *41*, 1907–1926. [CrossRef]

5. Cord, A.F.; Brauman, K.A.; Chaplin-Kramer, R.; Huth, A.; Ziv, G.; Seppelt, R. Priorities to Advance Monitoring of Ecosystem Services Using Earth Observation. *Trends Ecol. Evol.* **2017**, *32*, 416–428. [CrossRef] [PubMed]
6. De Araujo Barbosa, C.C.; Atkinson, P.M.; Dearing, J.A. Remote sensing of ecosystem services: A systematic review. *Ecol. Indic.* **2015**, *52*, 430–443. [CrossRef]

7. Matsushita, B.; Yang, W.; Chen, J.; Onda, Y.; Qiu, G. Sensitivity of the Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) to topographic effects: A case study in high-density cypress forest. *Sensors* **2007**, *7*, 2636–2651. [CrossRef]

8. Teillet, P.M.; Guindon, B.; Goodenough, D.G. On the slope-aspect correction of multispectral scanner data. *Can. J. Remote Sens.* **1982**, *8*, 84–106. [CrossRef]

9. Hantson, S.; Chuvieco, E. Evaluation of different topographic correction methods for landsat imagery. *Int. J. Appl. Earth Obs. Geoinf.* **2011**, *13*, 691–700. [CrossRef]

10. Tan, B.; Masek, J.G.; Wolfe, R.; Gao, F.; Huang, C.; Vermote, E.F.; Sexton, J.O.; Ederer, G. Improved forest change detection with terrain illumination corrected Landsat images. *Remote Sens. Environ.* **2013**, *136*, 469–483. [CrossRef]

11. Peng, D.; Zhang, H.; Yu, L.; Wu, M.; Wang, F.; Huang, W.; Liu, L.; Sun, R.; Li, C.; Wang, D.; et al. Assessing spectral indices to estimate the fraction of photosynthetically active radiation absorbed by the vegetation canopy. *Int. J. Remote Sens.* **2018**, *39*, 8022–8040. [CrossRef]

12. Galvão, L.S.; Breunig, F.M.; Teles, T.S.; Gaida, W.; Balbinot, R. Investigation of terrain illumination effects on vegetation indices and VI-derived phenological metrics in subtropical deciduous forests. *GIScience Remote Sens.* **2016**, *53*, 360–381. [CrossRef]

13. Maeda, E.E.; Heiskanen, J.; Aragão, L.E.O.C.; Rinne, J. Can MODIS EVI monitor ecosystem productivity in the Amazon rainforest? *Geophys. Res. Lett.* **2014**, *41*, 7176–7183. [CrossRef]

14. Maeda, E.E.; Galvão, L.S. Sun-sensor geometry effects on vegetation index anomalies in the Amazon rainforest. *GISci. Remote Sens.* **2015**, *52*, 332–343. [CrossRef]

15. Huete, A.R.; Didan, K.; Shimabukuro, Y.E.; Ratana, P.; Saleska, S.R.; Hutyra, L.R.; Yang, W.; Nemani, R.R.; Myneni, R. Amazon rainforests green-up with sunlight in dry season. *Geophys. Res. Lett.* **2006**, *33*, 2–5. [CrossRef]

16. Morton, D.C.; Nagol, J.; Carabajal, C.C.; Rosette, J.; Palace, M.; Cook, B.D.; Vermote, E.F.; Harding, D.J.; North, P.R.J. Amazon forests maintain consistent canopy structure and greenness during the dry season. *Nature* **2014**, *506*, 221–224. [CrossRef]

17. Galvão, L.S.; dos Santos, J.R.; Roberts, D.A.; Breunig, F.M.; Toomey, M.; de Moura, Y.M. On intra-annual EVI variability in the dry season of tropical forest: A case study with MODIS and hyperspectral data. *Remote Sens. Environ.* **2011**, *115*, 2350–2359. [CrossRef]

18. Zhang, H.K.; Roy, D.P. Landsat 5 Thematic Mapper reflectance and NDVI 27-year time series inconsistencies due to satellite orbit change. *Remote Sens. Environ.* **2016**, *186*, 217–233. [CrossRef]

19. Zhang, H.K.; Roy, D.P.; Kovalskyy, V. Optimal solar geometry definition for global long-term landsat time-series bidirectional reflectance normalization. *IEEE Trans. Geosci. Remote Sens.* **2016**, *56*, 3624. [CrossRef]

20. Roy, D.P.; Zhang, H.K.; Ju, J.; Gomez-dans, J.L.; Lewis, P.E.; Schaaf, C.B.; Sun, Q.; Li, J.; Huang, H.; Kovalskyy, V. A general method to normalize Landsat reflectance data to nadir BRDF adjusted reflectance. *Remote Sens. Environ.* **2016**, *176*, 255–271. [CrossRef]

21. Nagol, J.R.; Sexton, J.O.; Kim, D.H.; Anand, A.; Morton, D.; Vermote, E.; Townshend, J.R. Bidirectional effects in Landsat reflectance estimates: Is there a problem to solve? *ISPRS J. Photogramm. Remote Sens.* **2014**, *103*, 129–135. [CrossRef]

22. Sims, D.A.; Rahman, A.F.; Vermote, E.F.; Jiang, Z. Seasonal and inter-annual variation in view angle effects on MODIS vegetation indices at three forest sites. *Remote Sens. Environ.* **2011**, *115*, 3112–3120. [CrossRef]

23. Li, F.; Jupp, D.L.B.; Thankappan, M. Issues in the application of Digital Surface Model data to correct the terrain illumination effects in Landsat images. *Int. J. Digit. Earth* **2015**, *8*, 235–257. [CrossRef]

24. Adhikari, H.; Heiskanen, J.; Maeda, E.E.; Pellikka, P.K.E. Does topographic normalization of Landsat images improve fractional tree cover mapping in tropical mountains? *Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.* **2015**, *XL*, 11–15. [CrossRef]

25. Jin, H.; Li, A.; Xu, W.; Xiao, Z.; Jiang, J.; Xue, H. Evaluation of topographic effects on multiscale leaf area index estimation using remotely sensed observations from multiple sensors. *ISPRS J. Photogramm. Remote Sens.* **2019**, *154*, 176–188. [CrossRef]

26. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. *Remote Sens. Environ.* **2017**. [CrossRef]
27. Pasquarella, V.J.; Holden, C.E.; Kaufman, L.; Woodcock, C.E. From imagery to ecology: leveraging time series of all available Landsat observations to map and monitor ecosystem state and dynamics. *Remote Sens. Ecol. Conserv.* **2016**, *2*, 152–170. [CrossRef]

28. Young, N.E.; Anderson, R.S.; Chignell, S.M.; Vorster, A.G.; Lawrence, R.; Evangelista, P.H. A survival guide to Landsat preprocessing. *Ecology* **2017**, *98*, 920–932. [CrossRef]

29. USGS. Landsat Collection 1 Level 1 Product Definition; EROS: Sioux Falls, SD, USA, 2019.

30. USGS. Landsat 4-7 Surface Reflectance (Ledaps) Product Guide; EROS: Sioux Falls, SD, USA, 2019.

31. USGS. Landsat 8 Surface Reflectance Code (LASRC) Product Guide; EROS: Sioux Falls, SD, USA, 2019.

32. Foga, S.; Scaramuzza, P.L.; Guo, S.; Zhu, Z.; Dilley, R.D.; Beckmann, T.; Schmidt, G.L.; Dwyer, J.L.; Joseph Hughes, M.; Laue, B. Cloud detection algorithm comparison and validation for operational Landsat data products. *Remote Sens. Environ.* **2017**, *194*, 379–390. [CrossRef]

33. Didan, K.; Munoz, A.B.; Solano, R.; Huete, A. MODIS Vegetation Index User's Guide (MOD13 Series); The University of Arizona: Tucson, AZ, USA, 2015.

34. Smith, J.A.; Lie Lin, T.; Ranson, K.J. The Lambertian assumption and Landsat data. *Photogramm. Eng. Remote Sens.* **1980**, *46*, 1183–1189.

35. Holben, B.N.; Justice, C. The topographic effect on spectral response from nadir-pointing sensors. *Photogramm. Eng. Remote Sens.* **1980**, *46*, 1191–1200.

36. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The Shuttle Radar Topography Mission. *Rev. Geophys* **2007**, *45*, 1–33. [CrossRef]

37. Ponzoni, F.J.; da Silva, C.B.; dos Santos, S.B.; Montanher, O.C.; dos Santos, T.B. Local illumination influence on vegetation indices and plant area index (PAI) relationships. *Remote Sens.* **2014**, *6*, 6266–6282. [CrossRef]

38. Zar, J.H. *Biostatistical Analysis*, 5th ed.; Prentice Hall: Upper Saddle River, NJ, USA, 2014; ISBN 9780131008465.

39. Ponzoni, F.J.; Galvão, L.S.; Liesenberg, V.; Santos, J.R. Impact of multi-angular CHRIS/PROBA data on their empirical relationships with tropical forest biomass. *Int. J. Remote Sens.* **2010**, *31*, 5257–5273. [CrossRef]

40. Brede, B.; Suomalainen, J.; Bartholomeus, H.; Herold, M.; Brede, B.; Suomalainen, J.; Bartholomeus, H.; Herold, M. Influence of solar zenith angle on the enhanced vegetation index of a Guyanese rainforest. *Remote Sens. Lett.* **2015**, *7058*. [CrossRef]

41. Wu, J.; Bauer, M.E.; Wang, D.; Manson, S.M. A comparison of illumination geometry-based methods for topographic correction of QuickBird images of an undulant area. *ISPRS J. Photogramm. Remote Sens.* **2008**, *63*, 223–236. [CrossRef]

42. Saleska, S.R.; Wu, J.; Guan, K.; Araujo, A.C.; Huete, A.; Nobre, A.D.; Restrepo-Coupe, N. Dry-season greening of Amazon forests. *Nature* **2016**, *531*, E4–E5. [CrossRef]

43. Gao, X.; Huete, A.R.; Ni, W.; Miura, T. Optical-biophysical relationships of vegetation spectra without background contamination. *Remote Sens. Environ.* **2000**, *74*, 609–620. [CrossRef]

<figure><img src="image_21.png" /></figure>