Temporal Patterns in Illumination Conditions and Its E ﬀ ect on Vegetation Indices Using Landsat on Google Earth Engine

: Vegetation indices (VI) describe vegetation structure and functioning but they are a ﬀ ected by illumination conditions (IC). Moreover, the fact that the e ﬀ ect 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 Reﬂectance level 1 data, we evaluated the temporal patterns of IC and two VI, the Normalized Di ﬀ erence 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 di ﬀ erent 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 ﬁndings conﬁrmed a strong correlation between EVI and IC and less between NDVI and IC. We found a signiﬁcant 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.


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][8][9][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 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. 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.
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.
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. Two common VI were calculated for every image NDVI [2] (Equation (1)) and EVI [1] (Equation 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 C1 (6.0) and C2 (7.5) are the coefficients of the aerosol resistance term [33]. Two common VI were calculated for every image NDVI [2] (Equation (1)) and EVI [1] (Equation (2)) and added as new bands.
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 1 (6.0) and C 2 (7.5) are the coefficients of the aerosol resistance term [33].

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 = Illumination Condition (IC) = cos z × cos s + sin z × sin s × cos (a − o) 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.

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

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 = Illumination Condition (IC) = cos z × cos s + sin z × sin s × cos (a − o) 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.

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

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

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

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 Remote Sens. 2020, 12, 211 7 of 17 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.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18 (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.
(a) (b)  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. 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. Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 18  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 2 = 0.5715, p < 0.01) (Figure 11), whereas NDVI and IC showed a very small to almost no dependence (R 2 = 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.  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 2 = 0.5715, p < 0.01) (Figure 11), whereas NDVI and IC showed a very small to almost no dependence (R 2 = 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. 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 2 = 0.5715, p < 0.01) (Figure 11), whereas NDVI and IC showed a very small to almost no dependence (R 2 = 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.   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.

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

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).  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]. 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. 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 2 = 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. 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 2 = 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

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.

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.

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.