Downscaling of MODIS leaf area index using landsat vegetation index

Abstract Several organizations provide satellite Leaf Area Index (LAI) data regularly, at various scales, at high frequency, but at low spatial resolution. This study attempted to enhance the spatial resolution of the MODIS LAI product to the Landsat resolution level. Four climatically diverse sites in Europe and Africa were selected as study areas. Regression analysis was applied between MODIS Enhanced Vegetation Index (EVI) and LAI data. The regression equations were used as input in a downscaling model, along with Landsat EVI images and land-cover maps. The estimated LAI values showed high correlation with field-measured LAI during the dry period. The model validation gave statistically significant results, with correlation coefficient values ranging from relatively low (0.25–0.32), to moderate (0.48–0.64) and high (0.72–0.94). Limited samples per vegetation type, the diversity of species within the same vegetation type, land-use/land-cover changes and saturated EVI values affected the accuracy of the downscaling model.


Introduction
The Leaf Area Index (LAI) is an important parameter characterizing vegetation and is widely used to describe activities within an ecosystem, such as the monitoring of plant growth stages (Clevers J 1989). LAI data are also used in hydrological models operating at river basin or sub-basin level (Schumann 1993) and especially in the estimation of water related parameters (Pedro et al. 2011), in order to define the water loss by evapotranspiration, as well as the roughness for surface water flow. Ideally, input LAI data should be acquired at as high spatial resolution and frequency (depending on the examined vegetation type) as possible, in order to describe accurately the development of vegetation.
LAI data can be acquired either through fieldwork, or remotely, using satellite or airborne sensors. Measurement of LAI in the field is a laborious and costly process and can be achieved with direct (ground-based), or indirect methods, such as taking vertical photographs of the canopy with a hemispherical lens and processing the acquired photos with specialized software (Alexandridis TK, Stavridou, et al. 2013). Satellite-based LAI is provided on weekly to monthly basis by a number of organizations, such as the Food and Agriculture Organization (FAO, http://api.data.fao.org/1.0/query-api-builder.html), the European Space Agency (ESA, https://land.copernicus.eu/global/products/lai) and the National Aeronautics and Space Administration (NASA, https://earthexplorer.usgs.gov/). Although the above-mentioned satellite LAI products are characterized by high temporal resolution (8-10 days), their low spatial resolution (300-1000 m) makes them useable mainly in regional and global studies, while they are not suitable for local studies. At higher resolution, the recently developed Biophysical Processor tool available through the Sentinel-2 toolbox of the SNAP software (http://step.esa.int/main/toolboxes/snap/), offers the ability to automatically create Sentinel-2 imagery-based LAI data at a spatial resolution of 20 m. However, it is not possible to create historical LAI datasets prior the launch period of the Sentinel-2 satellite in the orbit (2015), thus, prior the date Sentinel-2 entered into operational acquisition mode. Therefore, for the acquisition of historical LAI data (prior 2015) of high resolution, alternative methods and data sources should be investigated.
The LAI, as well as the Vegetation Indices (VIs), are some of the most valuable tools with a global use to monitor and evaluate the progress during plants' growth. Various studies have demonstrated the use of VI's for the estimation of LAI. The most common index used in relevant literature for the estimation of LAI is the Normalized Difference Vegetation Index (NDVI) . Major disadvantage is that the NDVI is not sensitive at medium to large LAI values. The Enhanced Vegetation Index (EVI) on the other hand, is considered to be a more stable and accurate VI compared to the NDVI, while in general these two indices have similar performance no matter what crops are analysed (Colombo et al. 2003). Moreover, the EVI has improved sensitivity in high biomass regions and according to  the values of EVI are more useful than those of the NDVI in estimating LAI. Even though NDVI is considered to be more appropriate for estimating chlorophyll responsiveness of the vegetation, the EVI is more sensitive to the variations of canopy structure (canopy architecture, crown type, vegetation physiognomy and LAI) (Chen and Sun 2010). An additional advantage is that EVI normalizes and adjusts the reflectance in the red band as a function of the reflectance in the blue band, minimizing this way the residual atmospheric influence (Huete et al. 1997).
Several studies have attempted to estimate LAI at high resolution using as input satellite imagery. These studies are generally narrowed to only certain types of crops (Walthall et al. 2004;Clevers et al. 2017), certain types of forest vegetation (Eklundh et al. 2001), or are limited to narrow-sized study sites. The presented methods are usually dependent on in-situ observations for the training of the algorithms and limited to a few sampled locations. In relevant literature Onojeghuo et al. (2018) developed a methodology for MODIS LAI downscaling after analysing a 2-year timeseries of LAI and NDVI data. However, this method can only find application in rice paddies and may be influenced by the low sensitivity of NDVI at medium to large LAI values. Houborg and McCabe (2018) exploited a Cubesat Enabled Spatio-Temporal Enhancement Method (CESTEM) to optimize the utility and quality of very high-resolution CubeSat imagery, resulting to the daily retrieval of NDVI and LAI at 3 m resolution via the fusion of CubeSat, Landsat, and MODIS data. Although this methodology provides LAI data of very high resolution (3 m), the acquisition of the input data used (CubeSat) is considered as costly. Houborg et al. (2015) managed to retrieve LAI directly from at-sensor radiances of Landsat data, using a fully integrated system of radiative transfer models (applicable over most regions without the need for site-specific calibration) and information extracted entirely from image-based and readily available datasets. Houborg et al. (2015) proposed a LAI downscaling methodology integrating a rule-based regression tree approach for estimating Landsat-scale LAI using the MODIS LAI product at 1 km spatial resolution and the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM), to intelligently interpolate the downscaled LAI between Landsat acquisitions. In their method, they also use MODIS NDVI data in order to first generate a finer resolution LAI map (250 m instead the initial LAI map of 1 km spatial resolution), which could assist to improve predictability in heterogeneous agricultural landscapes and improve the accuracy of the downscaling model. Similarly, Yu et al. (2018) obtained high resolution LAI by using a regression tree approach and the spatial and temporal adaptive reflectance fusion model (STARFM). Sun et al. (2017) developed a regression tree method for retrieving Landsat LAI for grape yield prediction, using homogeneous and high quality MODIS LAI data as reference. Ma et al. (2019) implemented a revised version of STARFM to produce high-resolution LAI data for use in SWAT hydrological models. Zhang and He (2016) (Silleos et al. 2014).
This study attempted to develop a methodology for downscaling the coarse spatial resolution of MODIS LAI images to higher resolution, using VI data of MODIS and Landsat satellites. The focus was on developing a simple and easily applied downscaling methodology, taking into account all vegetation types present in a river catchment, the diversity of climatic conditions found in study sites around the globe, as well as the use of satellite imagery that can be accessible free-of-cost.

Study areas
Four catchment areas in Europe and Africa ( Figure 1) characterized by climatic diversity were selected as study sites, covering a wide selection of natural vegetation types and agricultural systems. This choice was made in order to investigate the performance of the proposed methodology under various climatic conditions, since intra-species variability of vegetation is highly climate-dependent. The main characteristics of each study area are summarized in Table 1.

Data collection and pre-processing
Modis lai-evi MODIS data (products MCD15A2-LAI and MOD13A2-EVI at 1 km spatial resolution) were downloaded for each study area for the years 2012 and 2013 (according to in-situ LAI data availability, see Table 4), from the United States Geological Survey platform (USGS, https://earthexplorer.usgs.gov).
All pixels of good quality, according to the MODIS Quality Assurance (QA) data, from the 8-day (MCD15A2) and the 16-day (MOD13A2) image composites were used in subsequent analyses, excluding those falling on -or close to -the borders of different vegetation types (mixed pixels). A total of 20 MODIS images (10 for the EVI product and 10 for the LAI product) have been used for the analysis, as presented in Table 2.
Landsat EVI. Several atmospherically corrected and geo-referenced (UTM/WGS84) Landsat 7 and 8 EVI images were downloaded for each study area and for dates coinciding with the field surveys (Table 3), through the Land Satellites Data System (LSDS) Science Research and Development (LSRD) platform provided by USGS (https://espa.cr. usgs.gov/). The Landsat 7 EVI images included lines of no data, due to the sensor's scanline corrector malfunction. Gap-filling was not attempted to avoid inserting errors (Alexandridis TK, Cherif, et al. 2013) during the comparison of the proposed methodology results with the field measurements. Mainly pasture for grazing and other single-cycle crops. Greenhouses also present. Land use/land cover data. A high resolution and a moderate resolution Land Use/Land Cover (LULC) map was acquired for each study area, through EU FP7 project "MyWater" (see acknowledgments) for the study purposes. The Landsat derived LULC maps followed the Land Cover Classification System (LCCS) developed by FAO. The initial classification of 12 LCCS classes was changed for the study purposes, reducing the vegetation classes to 6 types of biomes, similar to those used by the MODIS land-cover product (which is used as input in the MODIS LAI algorithm): irrigated crops, rainfed crops, broadleaved forests, needleleaved forests, shrubs and grassland, as presented in Figure 2. The vegetation type maps were produced using spectral classification techniques (Nunes et al. 2013), such as the Maximum Likelihood Classifier (MLC), the majority voting fuser classifier and the multi-stage classifier, on Landsat multispectral images and field surveyed data. Two Landsat images acquired in different seasons were used per study site in order to account for different phenological stages of vegetation and improve classification accuracy. The overall accuracy of the high resolution vegetation type maps was 89.43% (estimated kappa 0.8808) for Nestos, 88.16% (kappa 0.8331) for Rijnland, 86% (kappa 0.792) for Tamega and 85% (kappa 0.849) for Umbeluzi.
The moderate resolution (300 m) GlobCover land-cover product for year 2009 was also acquired for each study area, in order to be used together with the MODIS datasets. The choice of the GlobCover dataset instead of the relevant MODIS LULC product, was made in order to avoid the use of interdependent data in the proposed methodology, since the MODIS LAI algorithm uses as input the MODIS LULC product. Additionally, the GlobCover dataset uses a spatial resolution approximately three times finer than the MODIS pixel (300 m instead of 1 km), so it was also feasible to determine the effect of mixed vegetation pixels and further exclude any sample points representing such a state during the analysis of data (see section Data collection and pre-processing -MODIS LAI-EVI). The Globcover data were reclassified according to the 6 biome types also used in the Landsat LULC maps presented above (irrigated crops, rainfed crops, broadleaved forests, needleleaved forests, shrubs and grassland). In-situ LAI. Field measured LAI data of 2012-2013 were acquired for all study areas through EU FP7 project "MyWater" (see acknowledgments), recorded at a seasonal timestep (Table 4). Representative sites in each study area (22-37 locations, see Annex I for details) were selected, covering all types of vegetation present (according to the 6 vegetation classes presented in section Data collection -Land Use/Land Cover). For the selection of a suitable sample design strategy, four key characteristics were used: the number of sample sites required within each study area, the spatial distribution of the sample sites, the required size of each sample site selected, as well as the number of subplots required (Brogaard and Olafsd ottir 1997). A stratified random sampling approach was preferred rather than a random sampling approach, in order to maintain a high level of accuracy and additionally reduce the number of samples. The number of locations selected as samples is highly dependent on the size of the study area, as well as on its spatial variability. Following the guidelines of Townshend and Justice (1988) for the estimation of the required size of each sample site, the pixel size and the geometric accuracy of the satellite images used were taken into consideration. Three subplots were measured at each location, forming the corners (vertices) of a triangle with dimensions equal to 2/3 of the satellite image's pixel size and then the three results were averaged. In cases of sparse vegetation or plantation, where great variability was expected, an approach of taking measurements on the corners and the centre of each sample site was followed. As a result of the sampling design, error inherent to sampling ranges from 7.5 (for Nestos) to 12.5% (for Umbeluzi). The error could not be further reduced as this would lead to a higher number of sample sites to be visited during field survey thus increasing the cost of field work in terms of time, money, and effort.
A protocol for measuring LAI using hemispheric camera photographs was used in order to ensure that the quality of the measurements would be as high as possible (Alexandridis TK, Stavridou, et al. 2013). The technique used was to estimate the canopy thickness by measuring the light transmission through the canopy. Hemispherical (use of fish-eye lens) photographs in upward and downward direction were processed with specialized software to estimate the field LAI values. In order to ensure that LAI estimates would be comparable in time and space, it was considered essential to compensate for the variable sun irradiance throughout the day and year, as well as under variable tree canopies. Therefore, the exposure of the photographic camera used was calibrated relatively to the sky, and this exposure was increased by two stops during image acquisition under the canopy (Welles and Norman 1991). Underexposure in near horizontal directions was avoided by using a zenith angle range between 0 to 75 (van Wijk and Williams 2005). In cases of low stature vegetation, where distortions might be present to the LAI estimated values (van Wijk and Williams 2005), the hemispherical photographs were taken from above in a downwards direction (Br eda 2003). The hemispherical photographs, were processed with CAN-EYE (freely available for downloading at https://www4.paca.inra.fr/can-eye/) software for the downward direction and with HEMISPHER software (http://www.wsl.ch/dienstleistungen/produkte/software/ hemisfer/index_EN) for the upward direction, following a different approach for each of the two cases. The upward direction photographs were processed based on gap-fraction identification, separating the canopy (dark) from the sky (light) pixels. The photos of downward direction were processed by using reflectance ratios, in order to identify the features on the photographs based on the spectral behaviour of the examined wavelengths. The analysis of the downward direction photographs was performed by grouping them according to the date each photograph was taken and to the light conditions present. An automatic class pre-selection was applied, with the user identifying or correcting after this step the final sample of vegetation and non-vegetation pixels on a hemispherical photograph. As a next step, this sample was used for all the photographs of each group, with an image classification process following, in order to acquire photographs indicating two or three possible states: vegetative state, non-vegetative state and mixed state. When this step was finished, CAN-EYE software automatically calculated LAI values in multiple photographs, with the final results provided in a spreadsheet. During the analysis of the upward direction photographs, the pixels were first classified as sky (white) or canopy (black) by the application of a brightness threshold to each picture analysed. The optimal threshold was estimated using the LAI-2000 method, as it is reported to be closer to a standard reference (Stenberg 1996). After this step, the light transmission (T) was calculated as a proportion of the white pixels present in the hemispherical photograph. Then the average number of times that a light ray would touch the canopy when travelling a distance equal to the thickness of the canopy was estimated, providing this way the contact number (K), which is described by the equation: K ¼ -cos H ln T. Finally, the K values were integrated over the rings in order to estimate the LAI values.
Data processing MODISlandsat EVI data comparison. The MODIS and Landsat EVI products are computed according to the following equation: where G is the Gain factor (¼ 2.5), near infra-red (NIR), red and blue are atmospherically corrected values of surface reflectance, L (¼ 1) is the factor for the canopy background adjustment and C1 (¼ 6) -C2 (¼ 7.5) are the coefficients of the aerosol resistance term used for the blue band, in order to correct aerosol influences in the red band.
Since there are minor differences in the bandwidth values between the Landsat and MODIS bands participating in the formulation of the EVI equation (see Table 5), the EVI products of the two satellites intended to be used in the same downscaling model were compared, to examine their compatibility. Therefore, for two of the study areas (Umbeluzi and Tamega) and for dates with minimum possible cloud presence, Landsat and MODIS EVI images of the same period were compared. The rescaling methodology involved grid polygons (fishnets, see Figure 3) which were created per study area, MODIS LAI-EVI data quality assurance. Only the pixels indicated from the MODIS Quality Assurance (QA) data as of "good" quality were taken into consideration during the analysis. Additionally, extreme values for LAI (<0, >20) and for EVI (<0.1) were excluded, as they clearly did not represent any possible vegetation state.
Landsat data quality assurance and LULC classification. For each Landsat image, its corresponding cloud mask (C implementation of Function of Mask, CFMask) layer was used to exclude pixels of low quality and those not representing vegetation (e.g., covered by snow/ice). Extreme EVI values (<0.1) were excluded in this case too, as they clearly did not represent any possible vegetation state. Using the Landsat LULC maps (presented in section Data collection -Land Use/Land Cover), each pixel of the above-mentioned processed Landsat EVI images was classified according to one of the six vegetation types described in section Data collection -Land Use/Land Cover.
MODIS LAI-EVI regression equations for the downscaling model. The MODIS LAI-EVI relationship was examined through regression analysis, taking into account both the characteristics of seasonality (Day of Year, DOY) and the vegetation type (Alexandridis et al. 2019). The regression analysis was applied using a linear fit (y ¼ a þ b Ã x) for each available pair of a MODIS LAI and EVI image of the same period and for each study area, while the data were grouped per vegetation type.

Downscaling model
The estimated MODIS LAI-EVI regression equations per study area and vegetation type, along with the processed Landsat EVI images, were used as input data in the downscaling model (Figure 4), in order to estimate the LAI values at the Landsat spatial resolution. In this model, the equations calculated per study area, vegetation class and date during the MODIS LAI-EVI regression analysis, were used as formulas to estimate the LAI values for each Landsat pixel, using the corresponding Landsat EVI values (classified per vegetation type) as input in the regression equations.

Downscaling model accuracy assessment
The created LAI maps of all study areas and for all available time-periods were validated using the available in-situ measured LAI data (presented in section Data collection -In-situ LAI).

MODIS-Landsat EVI data compatibility
The comparison between EVI values derived from MODIS data and the equivalent values derived from Landsat data showed high R 2 (0.89 for Tamega and 0.87 for Umbeluzi, Figure 4. Flowchart of the methodology followed for the downscaling model. p < 0.001, see Figure 5). The mean difference was 0.031 (m 2 /m 2 ). In both cases examined, the differences between the MODIS and Landsat EVI products were minor and thus, the data were evaluated as directly comparable and compatible for use in the downscaling model.

MODIS LAI-EVI regression analysis
During the quality assurance check of the MODIS datasets, only pixels flagged as of "good quality" were retained. After the application of the regression analysis, the equations that connected MODIS LAI with MODIS EVI per date and vegetation type were generated for each test site. An example of the regression analysis results for the study area of Umbeluzi on DOY 193-2013 can be found in Table 6, while a detailed presentation of all results is available at Alexandridis et al. (2019). In all vegetation types of Umbeluzi, the probability value (p) was <0.0001. The R 2 ranged between 0.39 (lowest value for the broadleaved forests due to the mixed tree and grass/shrubs vegetation) and 0.79 for the grassland vegetation type. The moderate to high R 2 values found were estimated using a large number of samples (N) and the offset and slope factors had similar values in all cases examined.
In the related scatter plots (Figure 6), most of the vegetation types of Umbeluzi show low LAI and EVI values due to the dry climatic conditions present. Considering the rainfed crops, which represent the main agricultural activity of the study area, the moderate R 2 found is due to the existence of two different species of rainfed crops dominating the study area (maize and pulses). Shrubs and grasslands show the lowest pairs of LAI and EVI values. They nevertheless have high R 2 , while broadleaved forests have the lowest R 2 . The scatter plots of rainfed crops and broadleaved forests show that two different  clouds of sample points were found (in similar location of the plots, see Figure 6), implying the presence of two different vegetation species in each vegetation class. After examining all regression results between MODIS LAI and EVI, the highest R 2 was estimated for the irrigated crop type, while the grassland showed moderate to high R 2 in all study areas except for Rijnland, probably due to water inundation phenomena in this flat landscape. Grasslands are characterized by vertical and horizontal homogeneity as well as a high ground cover, and because of this, the view angle influence is considered to be minor for this vegetation type (Myneni RB et al. 2002). Grassland areas are regularly mown for hay production; however, this could not be verified with the available field information.

Downscaled LAI maps at the landsat resolution
The regression equations estimated using the MODIS LAI-EVI datasets were applied to the Landsat EVI data for the matching dates and vegetation types. LAI maps were created at the Landsat spatial resolution level for all Landsat scenes available (see Table 3). An example per study area is presented in Figure 7, in which minimum cloud contamination of the images allowed the model to estimate LAI values for almost the entire study area in Tamega and Umbeluzi, and over 2/3 of the study area in Nestos and Rijnland.
A visual comparison between the MODIS LAI map of DOY 177-2013 for the study area of Umbeluzi and the equivalent downscaled Landsat LAI map (DOY 177-2013) is presented in Figure 8. The two LAI maps are characterized by similar patterns of LAI values; however, the Landsat LAI is more detailed.

Accuracy of the landsat resolution LAI maps
Validation of the downscaling model was performed by comparing the LAI values calculated at the Landsat resolution with the available in-situ LAI data. The results of the correlation analysis are presented below per study area and a full list of the associated results is provided in Table 7.

Nestos
During April 2013, the 2 LAI images of DOY 096-2013 and 104-2013 included 4 and 5 in-situ sampled locations respectively (Table 7). The correlation coefficient (r) was high for both cases (0.72 for DOY 096-2013 and 0.93 for DOY 104-2013) and statistically significant at 0.05 level, even though the number of sampled pixels was low compared to other cases examined. The RMSE was low in both cases (approximately 0.3 m 2 /m 2 ), while the mean difference was also low (0.1-0.6 m 2 /m 2 ), revealing a systematic overestimation of LAI.
For the period of July (DOY 184-2013), r dropped to approximately 0.25, which was the lowest between all datasets examined. The RMSE was also high (1.484 m 2 /m 2 ), while the mean difference was much lower (0.054 m 2 /m 2 ) than the results of April (for both DOY 096 and 104).

Rijnland
The validation between in-situ and Landsat calculated LAI data was not possible for the period of June (DOY 2013-177), since none of the field measured locations coincided with a pixel of a Landsat calculated LAI value, due to heavy cloud contamination of the Landsat EVI image used as input.
During September, the map of DOY 265-2013 included 7 samples (Table 7), for which r was 0.8065, statistically significant at 0.05 level. The RMSE remained moderate (0.6511 m 2 / m 2 ), while the mean difference was high (1.776 m 2 /m 2 ) and statistically significant at 0.001 level. The number of available samples increased to 15 for DOY 273-2013 resulting into a moderate r (0.4801), statistically significant at 0.05 level. The RMSE increased to 1.15 (m 2 / m 2 ), but the mean difference decreased to approximately half the value of DOY 265-2013 (0.87 m 2 /m 2 ), statistically significant at 0.05 level. This is due to the fact that in some cases grasslands were falsely identified as agricultural areas; however, the associated in-situ LAI values were similar to those estimated at the high spatial resolution level. Tamega During May (DOY 147-2013), 21 samples were used for the correlation analysis (Table 7) giving moderate-low r (0.3246) and RMSE of 0.88 m 2 /m 2 . The mean difference of LAI values (0.455 m 2 /m 2 ) was statistically significant at the 0.05 level.
The measurements of September (DOY 243 and 251) showed the highest correlation within all areas located in Europe (Nestos, Rijnland and Tamega). Both the created LAI maps of that period correlated with moderate-high r (0.636 for DOY 243-12013 and 0.643 for DOY 251-2013) statistically significant at 0.001 level and over 30 samples were included in both cases examined. The mean difference was characterized by similarly low values (0.188 and 0.248 m 2 / m 2 , respectively), whereas the RMSE was just above 1 m 2 /m 2 in both cases.

Umbeluzi
During November (DOY 309-2012) the accuracy of the downscaled LAI was very low (Table 7), although the sample size was larger compared to other cases examined. From the 10 samples examined, 3 showed a LAI value difference greater than 2 m 2 /m 2 and fell into the mixed vegetation category according to the field measurements (shrubs, trees and grass), while Landsat classified these areas as rainfed crops in two cases and shrubs in the third.
Correlation analysis between the in-situ LAI values of March and those estimated at the Landsat level (DOY 071-2013) was not possible, since only 2 of the in-situ LAI locations coincided with the areas for which a LAI value was calculated, due to heavy cloud contamination of the Landsat EVI image used as input.
For DOY 175-2013, r was the highest found for all cases examined (0.935, statistically significant at 0.001 level), with the mean difference (À0.2754 m 2 /m 2 ) showing underestimation of the field-measured LAI. The RMSE was low (0.2013 m 2 /m 2 ).
For DOY 191-2013, 16 samples were available showing an r value of 0.874, statistically significant at 0.001 level, and a mean difference of À0.192 m 2 /m 2 , statistically significant at 0.05 level. The RMSE was low (0.312 m 2 /m 2 ), although it was slightly higher as compared to DOY 175-2013.

Discussion
This study developed an easily applicable methodology for enhancing the coarse spatial resolution LAI product of MODIS (1 km) to the Landsat resolution (30 m), using free-ofcost and accessible satellite data of VIs. The accuracy of the proposed downscaling methodology was evaluated across various seasons at four hydrologically and climatically diverse river basins in Europe and Africa. In-situ LAI measurements of 2012-2013 available for all study areas were used for the validation of the downscaling model.

Overall accuracy of the downscaling model
The patterns of LAI distribution in the downscaled LAI maps were evaluated as reasonable, after visually examining the image pairs. Correlation analysis between the downscaled and in-situ LAI for the same locations and dates, showed variability in the deviations from field measurements across study areas and seasons. In general, the correlation results were found to be statistically significant, even in cases with low availability of samples (less than 10). In all cases where a Landsat 7 image was used in the model, several sampled locations were unusable when they coincided with a data gap. In most cases examined, the mean difference showed overestimation of the field measured LAI compared to the calculated LAI, except for the Umbeluzi study area. The study area of Umbeluzi showed consistently better results compared to the others. This mainly relates to the climate of Umbeluzi, where during dry season the vegetation has uniformly low LAI and EVI values. Comparison of the correlation coefficients with a similar study (Silleos et al. 2014) showed improved results (for the same study areas and field-measured data), with high r values (>0.7) in 5 out of 11 cases examined. Additionally, there was a slight improvement in the moderate r values (0.4-0.7), as most of the cases examined in present study showed moderate values over 0.5.

Study limitations
Several global studies have validated the MODIS LAI product using related fieldwork, for a wide range of vegetation types and time-periods, with the accuracy of the product estimated at 0.66 RMSE (Root-Mean-Square-Error) LAI units in case all types of biomes were taken into consideration, falling to a value of 0.5 RMSE units, when broadleaved forests were excluded (Myneni R 2012). Researchers warn that there has to be increased awareness when using MODIS data, since the relation between LAI and VIs is not applicable everywhere and all the time, due to measurement geometry and spatial resolution of the examined plant canopies (Tian et al. 2000). Seasonality during the growth cycle of each vegetation type seems to be another crucial factor when regression analysis is applied between MODIS LAI and VI data. As  support, during the regression analysis between LAI and VI data of MODIS, R 2 values are constantly increasing during different growing stages, but they also warn that assumptions and conclusions should be made only per growth-cycle (over one year period), as year to other years comparison between LAI and VIs may be erroneous . Meteorological conditions during different seasons is also a characteristic affecting the overall quality of such data, however, by the 8-day image compositing method that MODIS uses, it is possible to eliminate the contamination from cloud cover (Knyazikhin et al. 1999).
A possible limitation of the proposed methodology could be the different days (periods) of data acquisition within the composite MODIS images used for the regression analysis. For instance, MODIS 8-day composite LAI values might describe vegetation conditions of a different DOY compared to the values of the 16-day composite product of MODIS EVI. The MODIS LAI algorithm chooses the day with the highest fPAR value to assign the LAI pixel values in the associated maps. On the other hand, the highest MODIS EVI value, also closer to nadir, is selected to represent the pixel values in EVI maps. Consequently, in vegetation types such as irrigated or rainfed crops, where agricultural activity is present, a high LAI value of a crop during the last days before harvesting might be paired with a low value of EVI for the same period, if the closer to nadir EVI value coincides with the days just after harvesting. The same phenomenon could also occur in dates of rapid leaf green-up occurring in a period of 10-12 days (Ahl et al. 2006).
Low accuracies might have also occurred due to saturation of the EVI signal at high biomass densities, which has been documented for other VIs too (Turner et al. 1999). Pixels in low resolution data are likely to contain an amount of radiative contribution from the background (Myneni RB et al. 2002), and while saturation could explain low R 2 in the regression analysis results of needleleaved forests vegetation type, this was not the case for shrubs and grassland vegetation types, since saturation is not likely to affect sparse biomes (Fensholt et al. 2004). Considering needleleaved forests, the results could also be affected by the presence of lower layers of vegetation in areas of dense canopies, causing unsuitability of reflectance to model vegetation (Fensholt et al. 2004).
An identified source of error was the misclassification of vegetation types across the two sources of LULC data used in the downscaling model (Globcover and Landsat), which was expected to a certain degree due to the different spatial resolutions used between the two products. This was present in all study areas and contributed to overestimation of LAI in cases where dense canopies were present (such as the broadleaved forests) and underestimation of LAI in cases of low and dry vegetation (e.g., dry grassland during the summer period). Relevant studies also support that the effect of misclassification is larger when data of coarse spatial resolution are used (Tian et al. 2000). Moreover, for various in-situ LAI locations it was proven that the vegetation type was different from the one indicated by the Landsat vegetation type maps, which might have affected the correlation analysis results associated with the validation of the downscaling model. LAI measurements from fisheye hemispherical images might have also affected the validation results of the downscaling model, since this technique has been reported to show an underestimation of 25-50% as compared to direct measurements (Br eda 2003). This is mainly due to the deviation from the assumption of random distribution of the leaves, that the processing software is based on. Another issue that could significantly lower the accuracy of LAI measurements is the number of hemispherical photos taken in each location (Weiss et al. 2004). However, this has been accounted for by taking from 8 to 15 hemispherical photos in each location. Finally, the issue of calibration for variable illumination conditions has been accounted for with the inclusion of a procedure for calibration, using the sky as reference.

Future work suggestions
Variability of species within the same vegetation type, but of a different phenology, could have affected the results of the regression analysis. Heterogeneity causes greater underestimation levels of LAI and the magnitude of underestimation increases as vegetation heterogeneity increases (Fensholt et al. 2004). Therefore, the ideal situation would be to classify the data that fall into that category into more detailed vegetation classes (depending on the diversity of the dominating species per study area), in order to improve the linearity of the LAI-EVI relationship and succeed in estimating more accurate LAI values through the downscaling model. For example, Houborg et al. (2015) in their downscaling methodology which uses similar datasets with the present study (MODIS LAI-NDVI and Landsat products) achieved in their validation results over homogenous agricultural areas in Nebraska (corn and soybean crops) high levels of agreement between their calculated high resolution LAI and in-situ LAI (R 2 > 0.80 for all cases and timesteps examined). Applying a regression model across scales (MODIS to Landsat) is involving all the noises existing in high-resolution input at the final high-resolution output. This may be a reason explaining some low performances. Hence, the application of a smoothing technique before applying the downscaling model on high-resolution data might have been beneficial for the accuracy of the relevant results. An alternative method for downscaling LAI in cases of weak regression between LAI and EVI could be to use Landsat EVI regression as weights to redistribute MODIS LAI values at the Landsat pixels. This has been used effectively for downscaling evapotranspiration maps at similar scales (Chemin and Alexandridis 2006;Alexandridis et al. 2014), but could also be effective for LAI.
In order to improve the accuracy of the downscaling model results, Cohen et al. (2003) propose coupling more than one VI in the regression analysis. Although this could improve the overall accuracy of the downscaling model, it is believed that the relevant data and the effort demanded to develop the methodology would dramatically increase in volume, since the present study was applied in four geographic regions and for various seasons during a whole year.

Conclusions
A methodology to downscale MODIS LAI data to the Landsat spatial resolution was demonstrated, proving that MODIS and Landsat data are compatible for integration in the same downscaling model. Limited number of pixels available for some vegetation types in specific DOYs, diversity of vegetation species within the same vegetation type and misclassification of vegetation types, were the most probable factors affecting the downscaling model's accuracy. The created LAI maps at the Landsat spatial resolution showed high visual similarity in the patterns of LAI values distribution, compared to the MODIS LAI maps. Validation using insitu LAI data showed moderate to high correlation for most cases, while lower accuracies were observed during the rainy season for all the study areas. The presented methodology is useful in the development of a tool for LAI estimation in space and time, that is adaptable to different land cover types and geographic areas and could be used for monitoring vegetation development of both natural and agricultural areas.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by the Research Committee of the Aristotle University of Thessaloniki (Greece), grant number 90773 "Improvement of the estimation of Leaf Area Index (LAI) at basin scale using satellite images".