Leaf area index retrieval using gap fractions obtained from high resolution satellite data: Comparisons of approaches, scales and atmospheric effects

International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Contents lists available at ScienceDirect

International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag

Leaf area index retrieval using gap fractions obtained from high resolution satellite data: Comparisons of approaches, scales and atmospheric effects Alemu Gonsamo ∗ Department of Geography, University of Helsinki, Gustaf Häallströmin katu 2, P.O. Box 64, FIN-00014 Helsinki, Finland

a r t i c l e

i n f o

Article history: Received 13 August 2009 Accepted 16 March 2010 Keywords: Atmospheric correction High resolution satellite data Large scale LAI inversion Scale effect Vegetation index

a b s t r a c t This study is aimed at demonstrating the feasibility of the large scale LAI inversion algorithms using red and near infrared reflectance obtained from high resolution satellite imagery. Radiances in digital counts were obtained in 10 m resolution acquired on cloud free day of August 23, 2007, by the SPOT 5 high resolution geometric (HRG) instrument on mostly temperate hardwood forest located in the Great Lakes – St. Lawrence forest in Southern Quebec. Normalized difference vegetation index (NDVI), scaled difference vegetation index (SDVI) and modified soil-adjusted vegetation index (MSAVI) were applied to calculate gap fractions. LAI was inverted from the gap fraction using the common Beer–Lambert’s law of light extinction under forest canopy. The robustness of the algorithm was evaluated using the ground-based LAI measurements and by applying the methods for the independently simulated reflectance data using PROSPECT + SAIL coupled radiative transfer models. Furthermore, the high resolution LAI was compared with MODIS LAI product. The effects of atmospheric corrections and scales were investigated for all of the LAI retrieval methods. NDVI was found to be not suitable index for large scale LAI inversion due to the sensitivity to scale and atmospheric effects. SDVI was virtually scale and atmospheric correction invariant. MSAVI was also scale invariant. Considering all sensitivity analysis, MSAVI performed best followed by SDVI for robust LAI inversion from high resolution imagery. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Land surface parameters retrieved from remote sensing are required for several applications such as: forestry, agriculture, landscape studies, land management, hydrological, meteorological and meso-scale weather circulation forecasts (e.g., Stern and Donald, 1961; Bounoua et al., 2006; GCOS, 2008; Gobron and Verstraete, 2008; Gonsamo Gosa, 2009). Among several land surface parameters, leaf area index (LAI) is one of the principal biophysical parameters in climate, weather, and ecological studies, and has been routinely estimated from remote sensing measurements (Asner et al., 2003). LAI is a part of the essential climate variables identified by the Global Climate Observing System (GCOS, 2006). LAI is a dimensionless value of the amount of foliage area of a vegetation canopy and is defined as one half the total radiation intercepting leaf area per unit ground horizontal surface area (Chen and Black, 1992; Gonsamo and Pellikka, 2008). The major physiological processes of vegetation including photosynthesis and evapotranspiration are determined by the vegetation biophysical parameters that describe the canopy structure as an exchange medium for energy and matter between the terrestrial ecosys-

∗ Tel.: +358 9 191 50770; fax: +358 9 191 50760. E-mail address: alemu.gonsamo@helsinki.fi. 0303-2434/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jag.2010.03.002

tem and atmosphere. Remotely sensed data offers a wide range of spectral, spatial, temporal, and angular resolutions to estimate LAI. Several numerical models require a continuous field of high spatial and temporal resolution LAI measurements due to heterogeneity and size of natural vegetation or agricultural patches, and the large seasonal dynamics of vegetation. To fulfil these needs, the LAI retrieval methods including the processing of the remotely sensed data are expected to be efficient and convenient for the modellers and end users. There are several global initiatives for ground-based measurement validation networks, and global satellite products and intercomparisons for mapping LAI on continuous bases (Weiss et al., 2004; Myneni et al., 2002; Morisette et al., 2006; Garrigues et al., 2008). The retrieval of LAI from remote sensing data is usually done based on statistical/empirical relationships between spectral vegetation indices (SVI) and ground-based measurements (Huete, 1988; Baret et al., 1989; Myneni et al., 1995). While these relationships work well under a particular viewing and illumination geometry and for specific vegetation phenology, they produce inaccurate results when applied over the broad range of land cover types and optical and geometrical conditions normally encountered in satellite images. Another approach for canopy LAI retrieval is to invert a physically based canopy reflectance models as they describe the transfer and interaction of radiation inside the canopy based on physical laws to provide explicit connection between the biophys-

234

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

ical parameters and the canopy reflectance (e.g., Vohland et al., 2010). However, there are no ideal canopy reflectance models that can be directly used for LAI inversion over large areas as the inversion process is ill-posed by nature due to measurements and model uncertainties since different combinations of model parameters may correspond to almost identical spectra (Combal et al., 2002). Additionally, the selection of the parameters for the inversion of canopy reflectance models is complicated, and parameterization technique typically relies on the existence of experimental data collected at the site of interest. In general, the success of LAI estimation from remotely sensed data remains cumbersome and there is always a need to calibrate remotely retrieved parameters with ground-based measurements. This led to the cascades of methodological studies for LAI retrievals from remote sensing observation (Richardson and Wiegand, 1977; Baret et al., 1989; Myneni et al., 1995; Gutman and Ignatov, 1998; Combal et al., 2002; Asner et al., 2003; Bréda, 2003; Xiao and Moody, 2005; Jiang et al., 2006; Sprintsin et al., 2007). Most of the LAI retrieval limitations associated with the empirical approaches can be addressed by using SVIs which normalize out the atmospheric, topographic, and scale effects. The normalized and scaled difference vegetation indices such as normalized difference vegetation index (NDVI, Rouse et al., 1974) and scaled difference vegetation index (SDVI, Jiang et al., 2006) minimize the atmospheric and topographic effects, and the scale effect, respectively. In view of soil background reflectance effect, Qi et al. (1994) put forward modified soil-adjusted vegetation index (MSAVI) with correction parameters which can be retrieved from the image scene. Thus, general outcome of interpretation is that normalized difference SVIs remain the only practical approach to the global analysis of vegetation using remotely sensed data owing to partial cancelation of the bi-directional, atmospheric and other interfering effects on observed radiances (Gutman and Ignatov, 1998). Based on these premises, a number of studies analysed SVI for estimation of fractional vegetation cover (fc ) aiming particularly to relate with LAI from satellite remote sensing with arguable convergence of results (Richardson and Wiegand, 1977; Baret et al., 1995; Gutman and Ignatov, 1998; Xiao and Moody, 2005; Jiang et al., 2006; Sprintsin et al., 2007). These studies ignored the sensitivity analysis of atmospheric, soil background and scale effects opposed to the major aim of finding LAI methodology with a lesser amount of parameterization. Moreover, the use of SVIs for gap extraction is limited often by the coarse resolution of satellite observations, which results in mixed pixels made up of more than one land cover types. Therefore, it was hypothesised that the robustness of canopy gap fraction estimations can be improved using high spatial resolution remote sensing imagery. This paper first presents the comparative analysis of the algorithms used to retrieve LAI and its sensitivity for soil background, atmospheric and scale effects. The second part of

this paper describes the validation of the retrieved LAI and interactive correction to optimize the algorithms. It should be noted that the LAI retrieved solely from SVIs is the ‘green leaf area index’ compared to the LAI reported in literature which is derived from optical field instruments as a structural parameter (Bréda, 2003). The green leaf area index derived from SVI is a photosynthetically active LAI which is strongly related with the canopy optical thickness, whereas LAI as a structural parameter from ground measurement represents ‘plant area index’ including both photosynthetic and non-photosynthetic canopy materials. In spite of this difference, ground-based LAI estimated from optical field instruments is being used as a ground-truth validation dataset even in recent studies (e.g., Ganguly et al., 2008; Houborg and Boegh, 2008). In this regard, ground-based measurements are used as comparison dataset rather than as the ground-truth. Besides, the effectiveness assessment was extended by comparing the retrieved LAI with MODIS LAI product, and 1D Radiative Transfer Model inversion. 2. Site description and datasets 2.1. Site description The study site is located in the Great Lakes – St. Lawrence forest in Southern Quebec, Canada. It is part of the Gatineau Park (Fig. 1), which is managed by the National Capital Commission (NCC) of Canada and centred at 45◦ 30 N, 75◦ 52 W. Gatineau Park is a vital ecological and conservation reserve, consisting of many endangered and rare wildlife and plant species, and is one of the largest and oldest Parks managed by NCC (NCC, 2005). More than 80% of the park’s approximately 36,000 ha is forested. The park is about 10 km by 50 km and is mostly temperate hardwood forest with a dominant overstorey of sugar maple (Acer saccharum Marsh.) and small patches dominated by American beech (Fagus grandifolia Ehrh.), trembling aspen (Populus tremuloides Michx.), and red oak (Quercus rubra L.). Small numbers of red maple (Acer rubrum L.), American basswood (Tilia americana L.), ironwood (Ostrya virginiana (Mill.) K. Koch), white ash (Fraxinus americana L.), black ash (Fraxinus nigra Marsh.), white birch (Betula papyrifera Marsh.), and black cherry (Prunus serotina Ehrh.) are also present. 2.2. Ground-based LAI measurements 20 m by 20 m plots were established in the Gatineau Park (n = 61) along two north-south oriented transects (Fig. 1), which their size was set by a semivariance analysis conducted on field and remote sensing data of two Quebec Ministry of Natural Resources (QMNR) permanent sample plots (Butson and King, 1999). Both transects have been part of ongoing research on monitoring forest damage, structure, health and succession following the ice storm of

Fig. 1. Location map of study area, a SPOT image (colour infrared composite) taken July 23, 2007, showing the Gatineau Park, Canada.

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

1998, and were sites of LAI determination from airborne imagery (Gonsamo et al., in press). A subset (n = 54) of previously inventoried forest plots that could be easily found in 2007 were selected for LAI measurements. Plot corners were surveyed using differential GPS to provide positional accuracy on the order of <1 m. The ground LAI measurements were collected using digital hemispherical photography between August 10th and 20th, 2007. The photographs were acquired using a high resolution (8 mega pixels) Nikon Coolpix 8800 VR digital camera equipped with a fish-eye Nikon FC-E9 lens adapter (Nikon Inc., Japan). In total, five photographs were acquired in each study plot, one at each corner and one at the centre. All photographic procedures are described in Gonsamo et al. (in press). LAI was computed for the range of 0–60◦ view zenith angle to reduce the growing effects of mixed pixels near the horizon, which result from light scattering and coarse resolution. All photographs were analysed using CAN EYE software (http://www.avignon.inra.fr/can eye/). CAN EYE allows interactively classifying RGB colours into vegetation elements and background. At the end of this process a binary image, background vs. vegetation elements (including both green and non-green elements) is obtained to retrieve gap fraction. LAI is inverted from gap fraction assuming an ellipsoidal distribution of leaf inclination (Campbell, 1986), and using the lookup table (LUT) from zenithal variation of gap fraction (Weiss et al., 2004). These LAI estimation procedures from hemispherical photography using CAN EYE software were proven to be a robust approach (Demarez et al., 2008; Garrigues et al., 2008). Further details on LAI derivation from gap fraction using hemispherical photographs can be found in Weiss et al. (2004). For simplicity, the effective LAI, which assumes random foliage distribution, was used in this study and is hereafter referred to as LAI. Effective LAI is a part which directly contributes to the canopy interception of incident radiation and therefore well related with SVIs. Plotwise LAI ranged from 2.6 to 5.7 with an average value of 4 and was close to normal distribution along the two sampling transects of Gatineau Park. 2.3. Remote sensing data pre-processing Radiances in digital counts in the green (500–590 nm), red (610–680 nm), near infrared (NIR) (780–890 nm) and shortwave infrared (SWIR) (1580–1750 nm) wavelength regions were obtained in 10 m resolution acquired on cloud free day of July 23, 2007, 10:49 local time by the SPOT 5 High Resolution Geometric (HRG) instrument for 60 km × 60 km image swath. The sun zenith ( s ) and azimuth (ϕs ) angles were 30◦ and 140◦ , respectively. The view zenith angle ( v ) was 20◦ . The SWIR band was originally acquired at 20 m pixel size and resampled into 10 m in order to match the other bands. The SWIR band was used only for classification of land cover types (Section 3.2).

235

2.3.2. Radiometric calibration Satellite radiances were obtained for SPOT HRG quantised calibrated pixel values in 8-bit digital numbers (DNs) in four spectral bands. The DNs were converted to at-satellite radiance based on: Lsat =

DN G

(1)

where Lsat is the satellite radiance in W m−2 sr−1 ␮m−1 for band , and G is band specific absolute calibration gain. Satellite radiances were converted to top of the atmosphere (TOA) reflectance following: TOA =

Lsat (ESUN /d2 )cos(z )

(2)

where TOA is TOA reflectance for band , d is the earth–sun distance in astronomical units, and ESUN is the mean solar exoatmospheric irradiance in W m−2 sr−1 ␮m−1 for band . 2.4. MODIS LAI product The MODIS Collection 5 LAI product was acquired in a form of HDF subset from the Warehouse Inventory Search Tool (WIST) client for searching and ordering earth science data from various NASA and affiliated centres (https://wist.echo.nasa.gov/). The subset LAI product was in 1 km resolution with 8 days composite based on the maximum FPAR measurement. The HDF subset has been reprojected and they matched the study area layout. The available LAI was downloaded from date 23/07/2007 (20–27/07/2007) to cover the same period of ground and SPOT data acquisitions. The final MODIS LAI and associated quality layers are clipped to preferred study area covering 40 km by 40 km (1600 MODIS pixels). MODIS LAI product quality flags were extracted for the subset as the time period of the data acquisition coincided with the peak of growing summer period. Under optimal circumstances, a look-up-table (LUT) method is used to achieve inversion of a three-dimensional radiative transfer model for MODIS LAI retrieval (Myneni et al., 2002). When this method fails to optimize a solution, a back-up algorithm based on a relationship between the normalized difference vegetation index (NDVI) and LAI is utilised (Myneni et al., 1995). The Quality Flags serve to determine the origin of the calculated value or mark pixels where no retrievals were made. In the subset of 1600 MODIS pixels used in this study, 87% were retrieved using the main LUT algorithm with or without saturation, 1% were retrieved with main LUT algorithm containing significant cloud, 6% were retrieved using the back-up algorithm and 6% were fill values containing non-vegetated land covers. Only pixels retrieved with main algorithm without cloud contamination were used in further analysis. 2.5. Simulated data using PROSPECT and SAIL model

2.3.1. Geometric correction The image was orthorectified using a 3 in. digital elevation model (DEM) obtained from Canadian Digital Elevation Data (CDED) (www.geobase.ca). The digital source data for CDED at scales of 1:50,000 were extracted from the hypsographic and hydrographic elements of the National Topographic Data Base and various scaled positional data acquired from the provinces and territories. The DEM was resampled based on the nearestneighbour interpolation into 10 m to match the SPOT HRG image resolution. The geometric correction model was based on the SPOT XS sensor colinearity equation. The SPOT image was registered using 139 ground control points collected from orthorectified digital aerial photography and national road network of Canada (www.geobase.ca) into UTM zone 18 with a NAD83 datum. The X and Y root mean square deviation was below 6 m (0.6 pixels).

The combined PROSPECT leaf model (Jacquemoud and Baret, 1990) and the SAIL canopy model (Verhoef, 1984) were used to simulate canopy reflectance data in order to test the robustness of the methodologies developed in this study and to derive correction factors as described in the latter sections. The PROSPECT model simulates the leaf optical properties from visible to mid infrared as a function of leaf characteristics. The SAIL model is a turbid medium model with the assumption of homogeneous semi-infinite medium crop canopy with Lambertian reflecting leaves. The canopy reflectance has been simulated for two wavelength bands centred around 645 and 834 nm, corresponding to red (R) and NIR bands of SPOT HRG, respectively, for a range of plausible input parameters. PROSPECT uses the following input parameters: chlorophyll a and b content (Cab , ␮g cm−2 ), dry biomass content (Cm , g cm−2 ),

236

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

equivalent leaf water content (Cw , g cm−2 ), and mesophyll structure parameter (N); and SAIL uses: LAI, leaf angle distribution (LAD), soil reflectance (Rs ) and external parameters like View, Sun zenith ( v ,  s ) and azimuth (ϕv , ϕs ) angles, hot-spot size parameter (HS), and ratio of diffuse to total incident radiation (DF). For this analysis, the range of input LAI was between 0.25 and 9 by interval of 0.25, and 12 as the largest LAI value. Three different soil backgrounds ranging from dark to bright were specified from Bowker et al. (1985). The other input parameters were kept constant: N = 1.55, Cab = 34.24, Cw = 0.0137, Cm = 0.0045, spherical LAD,  v = 0◦ ,  s = 30◦ , ϕv = 0◦ , ϕs = 0◦ , HS = 0, and DF = 0. The input parameters which have been used in the analysis are within the realistic range and adopted in other simulation studies (Jacquemoud and Baret, 1990; Jacquemoud et al., 2000). Several studies have investigated the effect of changing one input variable at a time for vegetation characterization using PROSECT and SAIL models (e.g., Jacquemoud et al., 2009). Therefore, in this study, the evaluation of the methods was restricted only to the changes of LAI values and soil reflectances which are required for determination of maximum LAI and background reflectance determination for the LAI estimation approaches as explained in the latter sections. 3. Methodology 3.1. Atmospheric correction It was anticipated that atmospheric correction would be a crucial step for any quantitative analysis of satellite observations. Two procedures of atmospheric corrections and TOA compared are (1) apparent reflectance model, (2) absolute correction using four dark object subtraction (DOS) methods, and (3) absolute correction based on radiative transfer model. The main aim of the atmospheric correction is to remove atmospheric effects on remotely sensed DN in order to get ground surface reflectance. The first procedure is explained in Section 2.3.2, as retrieval method for TOA reflectance estimate. Apparent reflectance model corrects only the effects caused by the solar radiance and sun zenith angle (Eq. (2)), ignoring the effects caused by atmospheric scattering and absorption. Subsequently, procedures which further account for atmospheric scattering and also atmospheric absorptions are presented. 3.1.1. DOS method Several techniques have been developed for image-based atmospheric correction by assuming cloud free images and uniform path radiance (Chavez, 1996; Song et al., 2001). Although more detailed atmospheric corrections are required for generic studies using in situ measurements of atmospheric parameters, image-based methods usually provide sufficient accuracy. It is often difficult to quantify atmospheric data during the satellite overpass time to apply pixel based variable path radiance corrections. Dark object subtraction (DOS) is the simplest yet most widely used imagebased atmospheric correction approach. The first step of DOS method for atmospheric correction is to convert remotely sensed DN values to at-satellite radiance based on band specific gain coefficients which comes along the image header (Eq. (1)). Then the next step is to convert apparent at-satellite radiance to surface reflectance which involves the correction of effects caused by solar and viewing angles, and the atmosphere. DOS approach corrects the path radiance caused by atmospheric molecules and aerosols. Song et al. (2001) studied the accuracy of four different DOS based methods with respect to classification and change detection. They found that the best overall accuracies were achieved by using the simpler DOS method, rather than the more complex atmospheric corrections that combine both atmospheric models and the dark object concept. The same nomenclatures as Song et al. (2001) have

Table 1 Settings for DOS corrections. Method

Tv

Ts

EDOWN

DOS1 DOS2 DOS3 DOS4

1 cos(v ) e−r /cos(v ) e−/cos(v )

1 cos(s ) e−r /cos(s ) e−/cos(s )

0 0 6S ␲Lp

been adopted in this study following the improvements of DOS approaches as: DOS1, DOS2, DOS3, and DOS4. The basic assumption for all four DOS approaches is that within the image some pixels are in complete shadow and their radiances received at the satellite are only due to atmospheric scattering (path radiance). In this study, techniques originally developed by Chavez (1996) were adapted with minor modifications. All four DOS approaches (DOS1, DOS2, DOS3 and DOS4; Table 1; Song et al., 2001) calculate surface reflectance from Eq. (3) using different assumptions for atmospheric transmittance along the path from the ground surface to the sensor (Tv ), atmospheric transmittance along the path from the sun to the ground surface (Ts ), and downwelling spectral irradiance at the surface due to scattered solar flux in the atmosphere (EDOWN) (Table 1). If the sky irradiance is assumed isotropic and only atmospheric scattering and absorption are considered, the general equation for describing atmospheric interactions and retrieving a Lambertian surface reflectance can be expressed as:  =

(Lsat − Lp )

(3)

Tv ((ESUN /d2 )cos(s )Ts + EDOWN )

where  absolute surface reflectance and Lp is estimated path radiance. The path radiance is estimated using Eq. (4) assuming 1% surface reflectance for dark object. This assumption explains the fact that very few targets on the Earth’s surface are absolute black, so an assumed 1% minimum reflectance is better used than 0%. For this study, the dark object radiance for all bands was retrieved from the deepest and clearest water found in the image scene. The path radiance is estimated as:



Lp = Lsat −



0.01 (ESUN /d2 )cos(s )Ts + EDOWN Tv

(4)



DOS1 assumes no atmospheric transmittance loss, the transmittance along the sun-target and target-sensor paths is assumed to be 100% (Tv and Ts = 1), and the downwelling solar spectral diffuse irradiance is ignored (EDOWN = 0). In DOS2, cos( v ) and cos( s ) are used as a substitute for Tv and Ts , respectively, whereas EDOWN is the same as DOS1. Unlike, Chavez (1996) and other studies who mainly used DOS approaches for nadir viewing Landsat TM images, Tv is approximated in DOS2 by cos( s ) as a surrogate for targetsensor transmittance computation. For DOS3, Song et al. (2001) recommended a more realistic interpretation of path transmittance assuming Rayleigh scattering atmosphere for Tv and Ts (Table 1). Optical thickness for such an atmosphere is defined as: r = 0.008569−4 (1 + 0.0113−2 + 0.00013−4 )

(5)

where  is wavelength in ␮m. For DOS3, EDOWN for Rayleigh atmosphere was estimated by 6S atmospheric radiative transfer code as zero aerosol optical depth at 550 nm (Vermote et al., 1997). DOS4 is the same as Song et al. (2001) DOS4 model, i.e., adding atmospheric aerosols for path transmittance computation assuming isotropic sky radiance. For path transmittance Tv and Ts in Table 1,  is solved iteratively as defined by Song et al. (2001) as:

  = −cos(s ) ln

1−

4[Lsat − 0.01((ESUN /d2 )cos(s )Ts + EDOWN )Tv /]



[(ESUN /d2 )cos(s )] (6)

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

237

3.1.2. Radiative transfer model Radiative transfer code developed by Vermote et al. (1997) called 6S (second simulation of the satellite signal in the solar spectrum) was used. 6S first produces three atmospheric correction coefficients called a, b, and c that convert the measured at sensor radiance in each band to atmospherically corrected reflectance (acr). Atmospherically corrected reflectance is then obtained using the following formulas: acr =

y 1 + cy

(7)

where y = aLsat − b

(7.1)

The middle latitude summer standard atmospheric model available in 6S code was utilised. This was due to the fact that the atmospheric characteristics data were not available for the time of satellite overpass and studies have shown that the effect of the atmospheric correction from standard atmospheric model is comparable to that obtained from an in site model (Moran et al., 1992). The continental standard aerosol model was used to define aerosol type. Aerosol optical depth (AOD) at 550 nm was interpolated from AOD (440) and Angstrom exponent ˛ (440–675 nm) by: AOD (550 nm) =

AOD (440 nm) [550 nm/440 nm]

˛(440–675)

Fig. 2. Land use map of the study area showing the analysis extent and the whole SPOT image scene.

(8)

The AERONET AOD at 550 nm was first integrated with spatiotemporal constraints that the SPOT scene is located to two nearest available AERONET sites and the measured AOD is computed from the closest time t satellite overpass. The two AERONET sites were: CARTEL, Sherbrook, Canada (45◦ 22 44 ’ N, 71◦ 55 51 ’ W, and 300 m elevation); and Egbert, Canada (44◦ 13 33 ’ N, 79◦ 45 00 ’ W, and 264 m elevation). No bi-directional reflectance distribution function (BRDF) correction was applied. 3.2. Land use map The 6S and geometrically corrected SPOT image was used to generate land use map of the study area. The classification method was based on a supervised algorithm using a conventional maximum likelihood procedure. The training sites (polygons) were selected visually from the scene and the data was classified into four general classes initially aimed to extract vegetated areas. Apart from vegetation areas, the land use and land cover classes include: urban areas (towns, buildings and roads), water, and wet areas (shallow water and wet land). The selected pixels of training sites (an average per class of 45,000 pixels) were evaluated for class separability using Transformed Divergence (TD) (Jensen, 1996). The scale of the TD values can range from 0 to 2000. As a general rule, if TD result is greater than 1900, then the classes can be separated. The minimum TD distance of the four classes was above 1999. Fig. 2 shows the final land use map of the study area. For purposes of LAI retrieval, classes such as urban, water and wet areas were masked in the SPOT image prior to inversion as discussed in the following sections. All the analysis both from SPOT and MODIS datasets are restricted to the extent of the classification image shown in Fig. 2.

and ecological applications is the normalized difference vegetation index (NDVI, Eq. (9)), which has been shown to be widely usable for vegetation analysis (Rouse et al., 1974). Several studies have successfully estimated gap fraction or fc from satellite remote sensing using the NDVI (Richardson and Wiegand, 1977; Baret et al., 1995; Gutman and Ignatov, 1998; Xiao and Moody, 2005; Jiang et al., 2006; Sprintsin et al., 2007). NDVI =

(9)

Jiang et al. (2006) demonstrated scale invariant fc retrieval approach from difference vegetation index (DVI, Eq. (10)) scaled using bare soil and dense vegetation DVI values. The study showed that the derivative vegetation index inferred from the linear spectral mixing of red and NIR reflectances which in turn are the scaled difference of DVI computed from bare soil and dense vegetation is more suitable and robust approach to compute fc over heterogeneous surfaces (Jiang et al., 2006). DVI = NIR − R

(10)

Because of diversity both in vegetation cover and soil optical properties in natural environment, soil background have great effect on the SVIs. Huete et al. (1988) introduced soil-adjusted vegetation index (SAVI), which requires a parameter, L. The determination of L needs some known information, such as vegetation cover. For different types of surface materials, one fitted value of L is not suitable in natural environment. In view of this, Qi et al. (1994) put forward modified soil-adjusted vegetation index (MSAVI, Eq. (11)). MSAVI replaces a constant, L, with L being the value that can be adjusted with the variation of vegetation cover. Given a possible range of L in SAVI from zero to infinity, the expression can be derived with mathematical induction method as follows:

3.3. Retrieval of LAI 3.3.1. Methods to retrieve fractional vegetation cover from vegetation indices Fractional vegetation cover (fc ) is the percentage of vegetation cover projected vertically onto the ground or in the context of a remotely sensed image, the fc values represent the percentage of vegetative cover present in each pixel or in a unit area (Gonsamo et al., in press). The most frequently applied SVI for agricultural

NIR − R NIR + R

MSAV =

2NIR + 1 −



(2NIR + 1)2 − [8(NIR − R)] 2

(11)

Earlier modification of SAVI by Baret et al. (1989) proposed Transformed SAVI (TSAVI) as a measure of the angle between the soil line and the line which joins the vegetation point and the intercept of the soil line. The vegetation point is defined as the vegetation reflectances in the red and NIR bands. This study has shown that TSAVI is exponentially related with LAI based on experimental data.

238

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Table 2 Fractional vegetation cover (fc ) determination from spectral vegetation indices (SVIs). Name

Formulation of fc

Gutman and Ignatov

(NDVI − NDVIback )/(NDVI∞ − NDVIback ) 2

Carson and Ripley

((NDVI − NDVIback )/(NDVI∞ − NDVIback ))

Baret et al.

(1 − (NDVI − NDVIback )/(NDVI∞ − NDVIback ))

SDVI

(DVI − DVIback )/(DVI∞ − DVIback )

MSAVI

MSAVI/MSAVI∞

0.6175

Description

Source

Dense vegetation mosaic-pixel model Based on radiative transfer model Based on semi-empirical relationship between vertical gap and NDVI Scaled difference vegetation index from DVI Soil background invariant fc determination from MSAVI

Gutman and Ignatov (1998) Carlson and Ripley (1997) Baret et al. (1995)

Jiang et al. (2006) This study

NDVIback and DVIback are the SVI values for the bare soil (LAI = 0); and NDVI∞ , DVI∞ , and MSAVI∞ are the SVI values for infinite LAI (asymptotic value) of a pure vegetation pixel, respectively.

MSAVI unlike TSAVI or SAVI, can be calculated directly without a scene dependent correction factor yet minimizing soil background variations. The approaches for retrieval of fc from NDVI, DVI and MSAVI are presented in Table 2. The fc retrieval is based on the assumption that the SVI value of a given pixel is the linear combination of SVI values of green vegetation and bare soil, weighted by their relative proportions. NDVIback and DVIback are the SVI values for the bare soil (LAI = 0); and NDVI∞ , DVI∞ , and MSAVI∞ are the SVI values for infinite LAI (asymptotic value) of a pure vegetation pixel, respectively. As discussed above, MSAVI was deemed to minimize soil background influence on the retrieval of fc thus its bare soil value was not incorporated in the equation (Table 2). 3.3.2. LAI inversion from high resolution remote sensing data It is very plausible to assume that a pixel is more homogeneous at high spatial resolution as SPOT HRG pixel size (Price, 1992). In such cases, the radiances attenuated by a certain amount of vegetation layers can be sufficiently described by LAI using Beer–Lambert’s law (Nilson, 1971; Campbell, 1986; Choudhury et al., 1994; Sellers et al., 1996; Baret et al., 1995; Gutman and Ignatov, 1998). Therefore, the fc estimated using the equations listed in Table 2 is a complement in unity of a canopy gap fraction (P0 ), i.e., fc = 1 − P0 . Nilson (1971) demonstrated from theoretical considerations that the P0 in canopies could be expressed as a simple exponential function of LAI following the Beer–Lambert’s law of light transmission under forest canopy: P0 = e−kLAI

(12)

Therefore, for a known value of fc and k, LAI can be calculated as: LAI = −

ln P0 ln(1 − fc ) =− k k

(13)

where, k is the extinction coefficient which is related to leaf spectral properties and angles in the canopy. Canopies with vertical leaves have k values of 0.3 or less, whereas canopies with flat, horizontal leaves have k values of 1.0 (Campbell, 1986). The range of k among plants may span 0.1 to 1.0 (Campbell, 1986). For a broadleaved forest, k typically ranges between 0.42 and 0.58 (Bréda, 2003; see summary of lists of different k values). In this study, k was assumed to be 0.5, which is a good approximation for a broadleaved forest considering the near nadir view. Additionally, k was determined using the hemispherical photograph measurement of LAI and gap fraction from aerial imagery in order to invert Eq. (12) (Gonsamo et al., in press). The estimated k ranges from 0.4 to 0.6 with an average of 0.5 which is in good agreement with assumption and previous studies (Bréda, 2003). Subsequently, LAI is inverted using Eq. (13) from the gap fractions derived from fc values listed in Table 2 (Richardson and Wiegand, 1977; Baret et al., 1995; Gutman and

Ignatov, 1998; Xiao and Moody, 2005; Jiang et al., 2006; Sprintsin et al., 2007; Gonsamo et al., in press). Foliage element clumping index (CI) measures the degree of aggregation of individual leaves into shoots, tree crowns and vegetation patches (Gonsamo and Pellikka, 2009). The retrieval of aggregation of canopy elements at various scales from remotely sensed data is practically not feasible. Clumping in the forest stand is typically related to stand density, which is expressed as the number of trees per unit area and is sometimes taken as a coefficient (Daniel et al., 1979). Nevertheless, measuring density by remote sensing is not easy and requires more rigorous analysis as even high resolution sensors have difficulty in detecting single trees in the stand and gaps within tree crown (Bunting and Lucas, 2006). There have been some studies mostly concentrated on the estimation of clumping index with multi-angular remotely sensed data (Lacaze et al., 2002; Chen et al., 2003, 2005) and detailed analysis using optical field instruments (Leblanc et al., 2005; Walter et al., 2003; Gonsamo and Pellikka, 2009). However multi-angular remote sensing data is not always easily available. The relationships and the meaning of CI estimates from various studies have been always ambiguous. In this study, CI was computed from SPOT HRG derived NDVI-based gap fraction in the same manner as finite-length logarithmic averaging technique (Lang and Xiang, 1986) as: CI =

ln[1 − (NDVI − NDVIback )/(NDVI∞ − NDVIback )] ln[1 − (NDVI − NDVIback )/(NDVI∞ − NDVIback )]

(14)

where [1 − (NDVI − NDVIback )/(NDVI∞ − NDVIback )] is the canopy mean gap fraction averaged over predefined analysis grid size, and ln[1 − (NDVI − NDVIback )/(NDVI∞ − NDVIback )] is the mean of logarithmic gap fractions of all SPOT HRG 10 m pixels within predefined analysis grid. This is one of the three commonly applied CI computation methods from optical data, the remaining two being one based on gap size distribution from ground measurements (Leblanc et al., 2005) and the other based on multi-angular remote sensing data (Lacaze et al., 2002; Chen et al., 2003, 2005). The finitelength averaging approach based on large scale gap fraction is the first attempt to map CI at larger scale using single view angle remote sensing data. This approach assumes that vegetation elements (tree crowns in this case) are locally randomly distributed, i.e., Poisson law is valid at the size of selected segment (SPOT HRG pixel size). The gap size here in Eq. (14) is defined as the dynamic range of gap fraction obtained from NDVI. The analysis grid size is 1 km by 1 km, equivalent of 1 MODIS pixel and containing 10,000 SPOT HRG pixels. Therefore, the CI is computed using 10,000 SPOT pixels by overlaying analysis grid of MODIS pixel size. As the logarithm of zero is undefined, for ‘empty pixels’ where vegetation is full field of view, to allow the calculation of the logarithm of gap fraction, a new value of gap > 0 is computed, using a local Poisson model by

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

assigning the maximum LAI to arbitrary maximum value of 12 in Eq. (12) (Gonsamo et al., in press). The physical meaning of the CI derived from SPOT serves as the large-scale clumping related to tree distribution in the forest. NDVI has a capability for partial self-cancelation of bi-directional, atmospheric, soil, shadow and other interfering effects (Gutman and Ignatov, 1998). Therefore, the only difference is from the difference of tree distribution as other effects are either partially cancelled or assumed to be not optically varying within 1 MODIS pixel size. The resulted CI from all vegetated grids ranges from 0.4 to 0.999 excluding 0 from open fields. In real vegetation stands, CI typically ranges between 0.50 and 0.97 (Chen et al., 2003) and in good agreement with our result. The mean CI from this study was 0.7 which is in very close agreement with Leblanc et al. (2005) which found the mean of 0.72 from closed deciduous forest over Canada. The calculated CI is not expected to be intrinsically related with the magnitude of NDVI. However, the CI is computed from the pixels which contain only forest and agricultural fields, therefore the higher NDVI pixels are from forest sites which are expected to be uniform and the CI is close to unity. There is also compelling potential that a better understanding of the relationship between clumping and density would be a convenient approach for extinction coefficient determination since stand density determination by remote sensing is likely to be easier than direct determination of the leaf angle distribution and the angle with which light transverses the canopy (required for calculating the extinction coefficient). In this regard, CI was related to the extinction coefficient. CI was also further related with scale dependency of NDVI and other SVIs-based LAI retrieval from remotely sensed image. Furthermore, the CI is used to correct the computed leaf area index by dividing LAI with CI and the resulting LAI is compared with MODIS LAI product retrieved from 3D Radiative transfer model. LAI reported in this study is not corrected unless stated ‘clumping corrected’ otherwise, as the method is experimental and yet to be validated with more analysis,.

239

Fig. 3. The effect of background and saturated NDVI variations for LAI estimation using hypothetical NDVI data for (A) varying NDVI background and saturated values, (B) varying NDVI background and fixed saturated values, and (C) fixed NDVI background and varying saturated values. [ %, %] are saturated and background percentiles from NDVI histogram, respectively.

4. Results and discussion 4.1. Background and saturated SVI determination for LAI retrieval Background reflectance can be a significant contributor to the canopy reflectance signal and care should be taken to avoid confounding effects among fc -SVI (e.g., NDVI, DVI and MSAVI)-LAI especially at low vegetation coverage with varying soil optical properties. The effect of NDVI∞ and NDVIback determinations for the LAI algorithm used in this study were first analysed in a hypothetical dataset. The effect is the same for all NDVI, DVI and MSAVI algorithms except that MSAVI does not require background values (Table 2). 400 cases of random hypothetical NDVI values were generated spanning from 0 to 1, which are theoretically expected values of vegetated environments. Eq. (13) and fc cover retrieval approach shown in Table 2 were applied on a hypothetical dataset to retrieve LAI. The assessment on the NDVI is presented on Fig. 3 from hypothetical data and on Figs. 4 and 5 from the real dataset. As shown in equations listed in Table 2, the NDVI and DVI background values are divided both from the nominator and denominator therefore have minor effects. Fig. 3(B) shows that the background values have insignificant effect for the LAI estimations. However, all the NDVI values below the cut-off bounds used for the background value in the histogram of the NDVI will result in 0 LAI. For the LAI value calculated between the NDVI∞ and NDVIback cut-off bounds, LAI is systematically higher by 0.2354 for 0%, 0.2139 for 1% and 0.1445 for 5% compared with 10% cut-off bound NDVIback . All these are in acceptable range particularly for larger LAI values.

The variation of saturated NDVI value changes considerably the LAI estimates particularly for higher values whether background NDVI is the same or different (Fig. 3(A) and (C)). The lower the saturated NDVI, the higher the LAI estimates and the more measurements will be assigned to maximum LAI those which have NDVI values above the saturated value. The difference obtained for LAI values calculated using the 100% and 99% cut-off bounds as a saturated value is below 1%. However, the contrast between the cut-off bounds of 100%, and 95% or less is directly proportional to the LAI values (e.g., LAI estimated 3 and 4.5 using 100% cut-off bounds were estimated as 4 and 9.3, respectively using 95% cut-off bounds, Fig. 3(C)). In general, from a theoretical point of view, for the equations listed in Table 2, the higher cut-off bounds are used as a saturated NDVI, MSAVI and DVI values, the lower LAI estimates will be. Subsequently, the effect of the background and saturated NDVI determination was analysed on the real dataset. As the one of the objectives of this study is to develop image-based robust algorithms, the focus was set on analysing the determination of both background and saturated values from the single date image scene. Some studies have acquired both parameters from independent field measurements which are usually not readily available and practical; others used theoretically fixed value from radiative transfer model (Carlson and Ripley, 1997); and the common practice is using the time series NDVI analysis (Gutman and Ignatov, 1998). However, leaf optical properties can vary during the year independently of LAI. Using theoretically fixed, or other values which are derived from independent dataset are very arguable as the calcu-

240

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Fig. 4. The variations of NDVI background and saturated values as computed using various atmospheric correction methods from SPOT image (A) from NDVI histogram of the whole analysis extent of the image, and (B) from NDVI histogram of the whole analysis extent of the image after masking out built-up and wet areas (Fig. 2). 0–10%, and 90–100% are background and saturated NDVI percentiles from NDVI histogram, respectively.

lated values can vary due to atmospheric effects (Figs. 4 and 5), sun-view geometry and other factors which are not intrinsically related with the vegetation and soil conditions on the ground. There is also compelling reason that both parameters should be determined from the data itself, since they are specific to the region and

Fig. 5. The effect of background and saturated NDVI variations for LAI estimation using various atmospheric corrections from SPOT image over transect 1 (A) and 2 (B). 100 and 0 percentiles are calculated as maximum and minimum NDVI values from the image subscene after masking out built-up and wet areas and the rest are calculated from the whole image subscene. The solid horizontal lines represent the mean ground LAI values. The error bars are standard deviations. [ %, %] are saturated and background percentiles from NDVI histogram, respectively.

season due to their dependency on the soil and vegetation types and vegetation chlorophyll content (Price, 1992). For instance, in deserts, background SVI can be retrieved reliably, and in evergreen vegetation saturated SVI can be retrieved reliably throughout the seasons (Gutman and Ignatov, 1998). This indicates that from measurements which are made in the peak of the growing season as this study, saturated SVI can be retrieved accurately while background values are less important as shown in Fig. 3. Therefore, studies have used the upper and lower 1–5% cut-off bounds from the NDVI histogram as NDVI∞ and NDVIback (Sellers et al., 1996; Tang et al., 2007; Gonsamo et al., in press). The histogram of the NDVI is affected by the mixture of soil, vegetation and other land cover types which are present in most pixels resulting in ‘mixed pixels’. Fig. 4 shows the NDVI cut-off bounds computed from 5 atmospheric correction methods and with TOA reflectance applied on the histogram of the whole image with and without masking nonvegetated pixels (see Fig. 2. for whole image scene analysis extent and mask). The atmospheric correction has considerable effect on the overall NDVI calculations while all the DOS approach resulted in very close values. The variations of the saturated values are relatively flat along the 90–99% cut-off bounds due to the fact that the forested area is in its growing peak. Due to high weight of the vegetated area (Fig. 2), the lower NDVI cut-off bounds for 1–5% are much higher than what is reported in other studies (0.02–0.05, Sellers et al., 1996; Gutman and Ignatov, 1998; Tang et al., 2007). Here again, Fig. 4(A) and (B) show that the saturated NDVI can be retrieved more accurately when the vegetation is in peak of growing season, the conclusion which is in good agreement with Gutman and Ignatov (1998). Atmospheric corrections have considerable effects on the calculation of NDVI in general and on NDVI∞ and NDVIback in particular. MSAVI∞ has similar sensitivity for atmospheric correction (not presented for the sake of brevity) but SDVI normalizes out the effect and therefore has no influence on the final LAI calculation (see the next section). NDVI derived from atmospherically uncorrected reflectance values (TOA) is lower compared to corrected value (Fig. 4). This is due to the fact that scattering tends to increase the amount of red radiation received by the satellite as red is more readily scattered in the atmosphere than NIR. The LAI calculated from varying atmospheric correction methods and NDVI∞ and NDVIback values are presented for transect 1 and 2 in Fig. 5(A) and (B). The errors induced due to the determination of NDVI∞ and NDVIback variations are higher than the atmospheric effect on LAI calculations. The [100%, 0%] in Fig. 5 represents the minimum and maximum NDVI values after masking out the build-up and wet areas whereas the other upper and lower 1–10% cut-off bounds were calculated from the entire analysis extent of the image. The mean ground LAI measurements were

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

241

Fig. 6. Histogram distributions of LAI estimations using varying gap fraction estimation techniques (Table 2) and atmospheric correction methods. (A) Gutman and Ignatov NDVI, (B) Carson and Ripley NDVI, (C) Barte et al. NDVI, (D) SDVI, (E) MSAVI, and (F) all five LAI retrieval methods from 6S atmospherically corrected image.

used for the evaluation of estimated LAI. In both transects, consistently, the LAI estimated from 6S atmospheric correction method with [100%, 0%] outperformed all other options (Fig. 5) whereas all the other cut-off bounds calculations and atmospheric corrections overestimated the LAI values. The [100%, 0%] are more reasonable to use after carefully classifying the image into vegetated (including open agricultural fields) and non-vegetated classes. This highly reduces the cut-off bound dependency on the histogram weight caused by high prevalence of very low or very high density of vegetation, and the heterogeneity of land cover types, and the uncertainty caused by the atmospheric correction. Hereafter, the [100%, 0%] cut-off bounds as saturated and background values from the image masked out for non-vegetated areas were used unless stated otherwise. 4.2. The effect of atmospheric correction on LAI retrieval The atmospheric correction is probably the most critical element for any quantitative remote sensing assessments. Remotely

sensed reflectance datasets are not free of errors and the transformation of sensor DN counts into surface reflectance values will likely be associated with some extent of uncertainty. The conversion of radiance data into surface reflectance in the visible and NIR spectrum is typically most sensitive to the aerosol optical depth and type of aerosol model (Vermote and Vermeulen, 1999). The nonlinear effect of different atmospheric corrections on SVIs has strong effect on the individual pixel, and on background and saturated NDVI values. The five atmospheric correction algorithms and TOA were applied to obtain the reflectance of SPOT HRG imagery over the study area. The reflectance of the visible bands (green and red) by every algorithm is reduced considerably in contrast with that of TOA (not presented here for the sake of brevity). The nonsystematic differences among the four bands of the reflectance obtained from the five atmospheric correction algorithms and TOA do affect the ratio computation of SVIs such as NDVI, MSAVI and probably SDVI. All NDVI and MSAVI based methods are highly affected by atmospheric correction (Fig. 6(A)–(C) and (E)). SDVI is

242

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Table 3 Comparison of estimated LAI as computed over 1 km grids using various averaging methods from SPOT image atmospherically corrected with 6S. Method

Statistics

1

Gutman and Ignatov

Mean STDEV Fit model R2

3.13 1.33 x

Mean STDEV Fit model R2

1.27 0.54 x

Mean STDEV Fit model R2

1.94 0.87 x

SDVI

MSAVI

Average LAI

2

Average SVI

1.03 0.57 y = 0.06x2 + 0.11x 0.984 1

1.13 0.55 y = 0.91x 0.978

3

Average reflectance

1.83 0.27 y = 0.17x + 1.31 0.656

1

1.13 0.55 y = 0.91x 0.977

4

Average SVI

3.00 1.96 y = 0.29x2 − 0.11x 0.946

1

1.94 1.14 y = 0.66x2 + 0.52x 0.977

5

Average reflectance

5.15 1.49 y = 0.13x2 + 0.21x + 2.95 0.708

1

1.94 1.14 y = 0.66x2 + 0.52x 0.977

1

1

1

1

1.60 0.85 1 y = 0.85x 0.961

1.72 0.87 1 y = 0.91x 0.981

2.11 1.27 1 y = 0.26x2 + 0.49x 0.970

2.26 1.32 1 y = 0.25x2 + 0.58x 0.980

STDEV stands for standard deviation. 1 LAI inverted in an original 10 m pixel size using image-based minimum and maximum as background and saturated SVI values, respectively; 2 LAI inverted from SVI values averaged over 1 km grid using the same background and saturated SVI value as 1 ; 3 LAI inverted from averaged reflectance values over 1 km grid using the same background and saturated SVI value as 1 ; 4 LAI inverted from SVI values averaged over 1 km grid using 1 km pixel size image-based minimum and maximum as background and saturated SVI values, respectively; and 5 LAI inverted from averaged reflectance values over 1 km grid using the 1 km pixel size image-based minimum and maximum as background and saturated SVI values, respectively. 1 y is the average 1 LAI as a dependent variable in the fit model.

not affected by the atmospheric correction (Fig. 6(D)). For all NDVIbased LAI estimation methods, the histograms of 6S corrected were very close to TOA. The maximum LAI difference averaged all over the image was occurred among MSAVI based methods (21.45%), followed by approximately the same value for NDVI-based methods (21.23%) and lastly SDVI based method (3.2%). SDVI appeared to be invariant for varying atmospheric correction methods. Nevertheless, as shown in Fig. 6(F) and the following sections, the differences in magnitude obtained among the methods are higher than the effect of atmospheric corrections. However, there was very strong correlation of estimated LAI among the three categories of methods, i.e., Gutman and Ignatov-SDVI (R2 = 0.94), Gutman and Ignatov-MSAVI (R2 = 0.95), and SDVI-MSAVI (R2 = 0.99). 4.3. The effect of spatial resolution on LAI retrieval The effect of spatial resolution on LAI retrieval was further analysed by applying varying methods of averaging of reflectance, SVI and LAI over 1 km grid on the data derived from [100%, 0%] cut-off bounds as saturated and background SVIs values from the image masked out for non–vegetated areas and corrected using 6S. Changing the scale of surface reflectance, or the SVI and retrieving the LAI using the same algorithms as the high resolution reflectance would alter the result due to the exponential relationships assumed between LAI-P0 , LAI-fc , LAI-SVIs, and LAI-reflectance (e.g., Lang and Xiang, 1986; Table 2; Eq. (12); Fig. 3). Ratio based SVIs such as NDVI are known to be sensitive to the scale effect due to nonlinear formulation of the indexes. On other hand, orthogonal SVIs such as DVI are less affected and study has shown that SDVI is not affected by the scale (Jiang et al., 2006). The results of Gutman and Ignatov NDVI, SDVI and MSAVI based methods were presented, as the other NDVI-based methods have the same response as Gutman and Ignatov NDVI approach. Three commonly used averaging techniques and two methods for the determination of saturated and background SVI values are described in Table 3. The reference method for the comparison was based on the LAI retrieved using high resolution dataset and averaged over 1 km grid. The other two averaging methods are: the SVIs are averaged over 1 km grid and the LAI algorithm is applied, and the reflectances are averaged over 1 km grid and the SVIs derived and LAI algorithm is applied. NDVI was found to be very sensitive for any scaling effect applied (Fig. 7(A)). All mean values of LAI obtained from the five averaging methods of NDVI-based method are statistically different (t-Test, p < 0.0001). The normalized variations (coefficient of variation, CV) were found to be different for each averaging methods, i.e., CV of

1 LAI = 0.42, 2 LAI = 3 0.56, 3 LAI = 0.15, 4 LAI = 0.65, and 5 LAI = 0.29 (see

Table 3 for the description of superscript). The two reflectance averaging methods seem not to be traceable to high resolution LAI value. As LAI is dimensionless variable (m2 /m2 ), all averaging methods are not expected to change the value of LAI at various scale levels. Particularly, the NDVI∞ and NDVIback calculated from the averaged reflectance from both the high resolution and the aggregated histograms are not expected to change the resulting LAI (3 LAI and 5 LAI, Table 3), as the averaged reflectance is equivalent to what would be measured by a sensor at this resolution. Up to this moment, the fundamental characteristics of NDVI averaging techniques have never thoroughly been discussed for LAI and other variables which are related non-linearly with NDVI. Some studies, however, addressed the spatial aggregation of NDVI as a function of land cover types (Hu and Islam, 1997; Jiang et al., 2006; Sprintsin et al., 2007; Zhou et al., 2009). There is still a great need for theoretically understanding of different averaging methods by considering spectral mixture model. The LAI retrieved from various averaging technique of NDVI, and reflectances for NDVI computations is not generally comparable, nor fc derived from the NDVI. Which averaging technique would result the NDVI that would be obtained from coarser resolution acquisition should be theoretically and experimentally defined from simultaneous acquisition of dataset at different resolution to avoid the miss-registration of the pixels. In what so ever cases, if NDVI is very sensitive for scale effect for LAI retrieval, using solely NDVI without supplementary information for LAI retrieval should be avoided as used in this study. Both SDVI and MSAVI are virtually scale invariant for LAI retrievals provided that the saturated and background values are the same (Fig. 7(B) and (C)). However, if the saturated and background values are computed from the aggregated SVI or reflectance, the resulting LAI increase exponentially in both cases compared to the high resolution LAI (Table 3; Fig. 7(B) and (C)). SDVI is more stable than MSAVI for both SDVI aggregation and reflectance aggregation (Table 3). Detailed discussion of the fundamentals of SDVI as a scale invariant index is given in Jiang et al. (2006). The assumption used for LAI inversion algorithm from the gap fraction obtained from various SVIs assumes the uniform pixel model. This is valid for 10 m pixel SPOT HRG data as it is considered high resolution and the pixels are uniform. However, NDVI of a mixed pixel as demonstrated in Jiang et al. (2006) does not vary in consistent or traceable manner with fc or LAI. This is shown also in Fig. 7 by applying varying averaging technique for NDVI-based LAI retrieval. Fig. 7(D) shows that there is a trend between the computed clumping index and the differences obtained among high resolution LAI and LAI from

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

243

Fig. 7. Comparison of estimated LAI values using 10 m pixel size SVI values (1 Table 3) on X-axis averaged over 1 km grids with: LAI inverted from SVI average (2 Table 3); LAI inverted from reflectance average (3 Table 3); LAI inverted from SVI average (4 Table 3); and LAI inverted from reflectance average (5 Table 3). (A) Gutman and Ignatov NDVI, (B) SDVI, (C) MSAVI, and (D) LAI values using 10 m pixel size SVI values (1 Table 3) on X-axis averaged over 1 km grids plotted against the computed clumping index. All the legends for superscripts: 1, 2, 3, 4, and 5 are given in Table 3.

various averaging methods shown in Fig. 7(A), (B) and (C). For example, regression analysis reveals that: (a) from NDVI-based method, 1 LAI = 0.4343(2 LAI/CI) (R2 = 1, RMSD = 0); (b) from SDVI based method, 1 LAI = 2(4 or 5 LAI/CI) (R2 = 0.98, RMSD = 0.06); and (c) from MSAVI based method, 1 LAI = 1.4251(4 LAI/CI) (R2 = 0.99, RMSD = 0.06). Therefore, the clumping index developed in this study can be used for correction of SVIs or LAI calculated using various averaging methods. Since there is no information about the extinction coefficient (k), for Eq. (12) and (13), it is indispensable to approximate k from the image itself therefore all the variables needed for the LAI calculation based on the information obtained from image is acquired analytically. k from satellite observation depends on the leaf orientation and distribution, as well as stand and canopy architecture. The former two are difficult to obtain from satellite observations. The stand and canopy architecture can be approximated to stand and canopy distribution which in turn can be expressed with the clumping index. Light extinction under forest canopy decays exponentially with amount of green biomass as does fractional of absorbed photosynthetically active radiation. As shown in Fig. 7(D), the clumping index also decays exponentially with LAI (the higher value indicating the lower clumping). This information from clumping index can be related with k for gap fraction correction rather than LAI, as also Baret et al. (1995) and Carlson and Ripley (1997) tried to correct fc . The gap fraction in Eq. (12) is related to contact

frequency: P0 = e−kLAI = e−N

(15)

where K is the mean number of contacts between a light beam and a vegetation element. −ln(P0 ) = kLAI = K

(16)

Carlson and Ripley (1997) described that the variation of NDVI with respect to the global LAI in partially vegetated areas would be mostly controlled by the variation in the fraction of vegetated surface area illuminated by the sun and visible to the radiometer. Verstraete and Pinty (1991) discussed the nature and extent of NDVI variations in partially vegetated lands, and argued that NDVI is more strongly controlled by changes in vegetation cover than by changes in the optical thickness of canopies. Therefore, theoretically, it is fair to assume that the NDVI-based method (Gutman and Ignatov) gives K rather than LAI if there is no correction for fc or gap fraction applied as demonstrated by Baret et al. (1995) and Carlson and Ripley (1997). Fig. 7(A) confirms these assertions as the LAI from reflectance averaging methods (3 and 5 Table 3) overestimated compared to high resolution LAI for lower LAI values (1 Table 3). The overestimations and the deviation of the LAI from high resolution LAI are clearly explained by the CI (Fig. 7(D)). Therefore, if the clumping index shown in this study is related with correction exponential parameter for gap fraction, the LAI from corrected gap

244

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Fig. 8. Comparison of estimated LAI values from corrected and uncorrected gap fraction (P0 ). (A) LAI estimated using 10 m pixel size NDVI values (1 Table 3) on X-axis averaged over 1 km grids plotted against LAI estimated from NDVI computed from averaged reflectance (5 Table 3) and corrected P0 using Eq. (17), and (B) LAI estimated from NDVI computed from averaged reflectance (5 Table 3) with corrected P0 using Eq. (17) on Y-axis plotted against LAI from uncorrected P0 (5 Table 3).

fraction without the knowledge of k will be: LAI = −ln(P0CI )

(17)

If CI = 1, then the average, non-intercepted fraction (P0 ) of a beam of radiation travelling is related to fc from NDVI without correction, therefore the corrected and non-corrected P0 are the same. Baret et al. (1995) showed that the exponential correction parameter varies depending on the density of vegetation, and it becomes small in areas with a low vegetation density which is in good agreement with the CI computed in this study (Fig. 7(D)). The corrected gap fraction in Eq. (17) may explain the confounding effects of crown clumping, extinction coefficient and various averaging differences. Eq. (17) was applied for the gap fraction data obtained from NDVI by averaging reflectance over 1 km grid (5 Table 3). The result is compared with the high resolution LAI estimated from NDVI based on constant k (1 Table 3) (Fig. 8(A)). The result shows exponential relationships between two LAI estimations with the similar trend of reflectance averaging LAI with constant k (Fig. 7) with high resolution LAI. The exponential relationship is mainly due to the differences of the NDVI∞ and NDVIback . The LAI values were lowered systematically when the gap fraction is corrected compared with uncorrected gap fraction and constant k of the same reflectance averaging method (Fig. 8(B)). The estimated LAI from corrected gap fraction is in average 42% lower than the high resolution LAI. The corrected gap fraction method is believed to be a better estimate since the high resolution LAI without correction overestimates the true LAI as shown in the following section. This study has shown how the CI from the gap fraction derived based on NDVI is related with varying averaging methods for LAI and SVI calculations and approximation of gap fraction correction in attempt of replacing extinction coefficient. However, strictly speaking, the CI cannot be directly used as a real foliage clumping index or extinction coefficient.

as the average deviations of SVIs-based LAI from ground-based LAI, and the percentage of root mean square deviation (RMSD), standard deviation and bias of all five SVIs methods are presented in Table 4. The mean error of MSAVI was the lowest and the estimates were the closest to the ground-based measurement with smallest underestimation (-4% of bias) and with coefficient of variation (CV) close to the ground-based measurement (Table 4). MSAVI also resulted in least RMSD of 16.14%. SDVI resulted in the highest absolute bias (Table 4) indicating that the method underestimated the LAI. Jiang et al. (2006) also found out that SDVI has underestimated fc cover compared to NDVI-based methods. The underestimation of the SDVI was attributed to a non-linear increase of NIR reflectance for higher fc due to rapidly accumulating green biomass with only gradual lateral percent cover increase (Huete et al., 1985). However, SDVI based method explained the variation of LAI on the field better than all other methods (Table 4). All NDVI-based method (Gutman and Ignatov, Carson and Ripley and Baret et al., the latter two being derived from the first) resulted in very narrow range of LAI compared to the ground-based and other SVIs-based methods. NDVI-based method without correction (Gutman and Ignatov) was the only method which resulted in positive bias. Nevertheless, except that of the coefficient of variation, the other evaluation parameters (Table 4) may give bias for the comparison of SVIs-based LAI with ground-based measurements due to the difference of the definitions assumed and practically could be achieved. Hence, the evaluation assessment was further extended by comparing with MODIS LAI product and applying the methods developed here in the independently simulated data.

Table 4 Evaluation of different methods to derive LAI (n = 54).

4.4. Validity of high resolution LAI (SPOT) estimations as compared to ground measurements Since the entire model parameters for LAI estimations used in this study were derived from the image itself, the validity assessment of the product is crucial. The LAI estimated from SVIs corresponds to the green part of the vegetation which is slightly different from the classical definition and the LAI obtained from the ground-based hemispherical photography which is “plant area index”. Therefore, the validity assessment was focused more on the comparison rather than direct validation. The mean error calculated

Method

Mean error (%)

RMSD (%)

CVa (%)

Bias (%)

Gutman and Ignatov Carson and Ripley Baret et al SDVI MSAVI

28 21 5 36 2

29.93 17.98 28.01 40.24 16.14

3.21 3.21 4.14 12.45 12

25.01 −7.32 −22.81 −37.30 −4.23

a

Coefficient of variation (CV) of ground-based measurement = 16.39%. The percent mean error, root mean square deviation (RMSD), CV, and Bias are defined below. P is predicted and O is ground-based measured LAI value of plot i and n is the number of measurements and  is average value. Mean error (%) = (1/n)

n

[(Oi − Pi )/Oi ·

i=1  1/2 n n 2 100%], RMSD (%) = (1/n) (Oi − Pi ) /(1/n) Oi · 100%, CV (%) = i=1 i=1  1/2 n n n 2 (1/n) (O − ) /(1/n) Oi · 100%, and Bias (%) = (1/n) Pi − i=1 i=1 n i=1 i n

(1/n)

i=1

Oi /(1/n)

i=1

Oi · 100%.

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

245

Fig. 9. Comparison of MODIS LAI product with SPOT derived LAI based on (A) Gutman and Ignatov, (B) Carson and Ripley, (C) Baret et al., (D) SDVI, and (E) MSAVI methods. Upper row are not corrected and lower row are corrected with clumping index. The comparison with MODIS is made only for those non-fill pixels derived from the main LUT method.

4.5. Comparison with MODIS LAI product The comparison between MODIS LAI product and LAI estimated using the five formulas shown in Table 2 are given in Fig. 9. The comparison is made only for those pixels which are retrieved with main algorithm without cloud contamination. The LAI from SPOT data is derived using the high resolution and averaged over the MODIS pixels. Fig. 9 (lower row) shows the SPOT LAI which is corrected for clumping index. It is theoretically possible that the LAI of a single MODIS pixel is not retrieved accurately due to uncertainties in inputs and the algorithm. MODIS LAI product is overestimated compared to all five methods of LAI retrieval used in this study. However, the correlations are statistically significant and RMSD are reasonably low. The deviation from regression line in Fig. 9 (above) can partially be explained by the uncertainties in coregistration between MODIS and SPOT pixels. A statistical comparison of the overall averages shows that, MODIS LAI overestimated by 43% for Gutman and Ignatov, 100% for Carson and Ripley, 131% for Baret et al., 232% for SDVI, and 134% for MSAVI methods. For the LAI corrected for clumping index, MODIS LAI overestimated by 4% for Gutman and Ignatov, 46% for Carson and Ripley, 68% for Baret et al., 140% for SDVI, and 70% for MSAVI methods. The overestimation is very severe for higher LAI values (Fig. 9). During the late summer, MODIS LAI products show great and progressive overestimation (Pisek and Chen, 2007). MODIS LAI products are known to overestimate LAI, for example, up to 75% (Yang et al., 2006), 200% over Canada (Abuelgasim et al., 2006), by 2–3 LAI over BOREAS study area in Canada for broadleaved forest (Fang and Liang, 2005) and up to 323.6% across north American forest sites (Heinsch et al., 2006). MODIS LAI resulted in large standard deviation compared to any of the methods used in this study and in agreement with Pisek and Chen (2007). The spatial distribution of MODIS LAI measured as a coefficient of variation was 0.33, with the five methods without correction for clumping was 0.34–0.39 and with correction 0.23–0.28. This indicates that the range of LAI with both MODIS and the five methods are in good agreement to explain the distribution of LAI. There was better relation with MODIS LAI and the LAI estimated using SPOT data for lower LAI values where the MODIS pixels are retrieved without saturation. The overestimation in higher LAI values and in general can be explained by prevailing amount (72%)

of MODIS pixels were retrieved with saturation. The new clumping index might help better relate the LAI derived from SVIs with structural property of the forest canopy, and the geometry and the shade of the crown distribution. MODIS uses a stochastic radiative transfer model to account for 3D effects of vegetation heterogeneity (clumping and species composition) to assign a single clumping value to each cover type in LAI algorithm (Wang et al., 2004). This does not account for its spatial variation and relies on biome classification and the level of abstraction inherent in the classification procedure. Therefore, theoretically, however, there was no reason that the LAI corrected with clumping should be better related with MODIS product. However, the correction has improved the overall mean agreement with MODIS LAI but decreased slightly the correlation and increased the RMSD. The algorithm for the clumping index should be further tested particularly with high resolution multi directional remote sensing data and refined with more ground control parameters if there is any need for such kind of parameter for functional–structural plant models to represent physiological processes in 3D tree representations. There are few studies which estimated clumping index from remotely sensed multi-angular data (Lacaze et al., 2002; Chen et al., 2003, 2005; Leblanc et al., 2005) even though the meaning of the index is often arguable. The methodology presented in this study for clumping index estimation have a physical meaning to explain the architecture of the vegetated canopy, which remains to be discussed in further studies. 4.6. Correction factor for SVI-LAI relationship and robustness assessment using simulated data The NDVI, SDVI and MSAVI based methods applied on the simulated spectrum data in order to assess the performances and derive the correction parameters for fc -LAI relationships as implemented in Baret et al. (1995) and Carlson and Ripley (1997). NDVIback and DVIback were calculated as average SVIs from the three soil spectrum as described earlier. NDVI∞ , DVI∞ , and MSAVI∞ were calculated with PROSPECT + SAIL simulated spectra with an LAI value of 12. The interactive correction was determined by applying different exponential factors on the derived fc for the three methods. The performance of interactive determination for correction factor was controlled using the true and estimated LAI regression slope, R2 , and

246

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Fig. 10. The performance evaluation of the three LAI retrieval methods as the simulated true LAI is plotted against retrieved LAI based on NDVI, SDVI and MSAVI methods. The ‘corrected’ indicates that the correction factor is applied to fractional vegetation cover in order best fit the retrieved LAI with true LAI. The solid line is the one-to-one line. The three distinctive strips of scattering points show the three soil background used for simulation.

RMSD. The factors which gave the regression slope close to unity, highest R2 and lowest RMSD were deemed to be the correction factor. As shown in Fig. 10, SDVI and MSAVI based methods were found to be the best fit without any correction and outperformed NDVIbased method. Therefore, there was no need to apply any correction factor for SDVI and MSAVI based methods. LAI estimated using NDVI-based methods both with and without correction showed sigmoidal relationships with true LAI with smooth curves. The sigmoidal relationship between the estimated and true LAI from NDVI-based method can be explained by that as NDVI is ratio based, it is very susceptible for scaling using background and saturated values for the exponential fit of LAI estimation. Based on the interactive determinations, the correction factor for LAI estimation using NDVI-based method was 6.15, i.e., the fc derivation should be corrected as [(NDVI − NDVIback )/(NDVI∞ − NDVIback )]6.15 . The correction lowered the RMSD, increased the R2 and resulted in the regression slope close to unity (Fig. 10). However, the correction factor has underestimated the lower LAI and overestimated the higher LAI which were otherwise without correction (Fig. 10). The mean LAI obtained using the three methods computed over the three soil types for uncorrected models are statistically different (t-Test, p < 0.001). The soil has small effect for the LAI (Fig. 10), however except for that of NDVI-based method, the mean LAI for the three soil types are statistically different (t-Test, p < 0.001). The differences of mean LAI among the three soil types were 0.5 3% for NDVI, 6 22.5% for SDVI and 4 15% for MSAVI based methods. These indicate that SDVI and MSAVI based methods are vulnerable for the influence of the soil; although, the differences are in acceptable ranges.

The relationships between the gap fractions derived from the three SVIs with LAI are robust for the simulated dataset. However, the result is dependent on the choice of the input background and saturated SVI values as explained in Section 4.1. The exponential correction factors which would be applied on fc derived using the equations in Table 2 in order to best fit best LAI estimations were however found out to be less of importance. Using entirely different approaches and data sources, Choudhury et al. (1994) and Gillies and Carlson (1995) and Carlson and Ripley (1997) independently obtained identical square root relationship between fc and NDVI-based method (Table 2). Baret et al. (1995) obtained varying correction factors both for the same and among different SVIs for the simulated and real datasets. Baret et al. (1995) had studied solely gap fraction without inverting into LAI, and the others have studied fc (Choudhury et al., 1994; Gillies and Carlson, 1995; Carlson and Ripley, 1997). However, there is no uniform relationship between fc and LAI. Besides, the exponential correction factors are dependent on the choice of different combination of background and saturated SVI values. This indicates that, using any of the generic correction factors introduce biases and should be avoided. The focus rather should be on the choice of appropriate background and saturated SVI. The vegetation SVIs which give indication of the optical thickness of vegetated biomass can be exponentially related with LAI using Eq. (12) without any correction. As demonstrated in this study and earlier by Baret et al. (1995), applying the exponential correction term which is fitted in a certain season, vegetation type, and atmospheric condition and correction approach, may not be directly applicable to other areas.

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

5. Conclusions On the basis of simulated datasets, ground-based and satellite measurements, the validity of an approach for the computation of LAI based on the information solely contained on image scene itself was proven to be reasonably accurate. The approach uses the classical red and NIR radiometric channels and is applied to mostly temperate hardwood forest and open fields. The methodology developed here can be extended to other types of vegetation, provided that there is accurate estimation of SVI for the upperend value of saturated LAI. The role of accurately specifying the background and saturated SVIs, the spatial scale and atmospheric effects, the robustness test, and fitness with global LAI product was considered for effectiveness assessments. The simple approach was also demonstrated to characterize the spatial heterogeneity of canopy architecture from a single view angle radiometric data. The varying definitions and assumptions used for LAI obtained from ground-based measurements, SPOT image retrieval and MODIS product, and the validity of using simple 1D radiative transfer model for robustness assessment make any complete validation of the approaches almost impossible. The ground-based measurements give the “plant area index” which includes nonphotosynthetic parts of plants such as braches and stems. On the other hand, the LAI retrieved from SPOT imagery is related to the “green leaf area index” which is highly relevant from application and vegetation function point of view for photosynthesis, evapotranspiration and carbon balance studies. The green LAI can be accurately estimated from the image as shown in this study. The most important parameter for LAI determination was the saturated SVI value as the background SVI has very small effect for the final LAI. The saturated SVI value can be retrieved from the image particularly from large swath imagery during the growing season. SDVI was found out to be both scale and atmospheric correction invariant provided that the atmosphere all over the image scene is assumed to be constant. This is important attribute towards developing simple algorithm for LAI retrieval as the atmospheric correction can be rigorous process. MSAVI is also scale invariant. NDVI was found to be not appropriate index for large scale LAI inversion due to the sensitivity for atmospheric correction and scale. In the cases of well defined a priori information are lacking for radiative transfer inversion, or, sufficient ground measurements for robust empirical equations are not available, the method shown in this study can be utilised to retrieve LAI on operational bases particularly during the growing season. Acknowledgements The research was funded by the Academy of Finland through the TAITATOO-project and personal grant to Dr. Petri Pellikka (120190), and by Natural Sciences and Engineering Research Council of Canada funding to Dr. Douglas King. The SPOT image was obtained from OASIS program financed by the European Commission, DG Research. AG is financed partly by CIMO, Department of Geography of University of Helsinki, and the Finnish Graduate School of Geography. References Asner, G.P., Scurlock, J.M.O., Hicke, J.A., 2003. Global synthesis of leaf area index observations: implications for ecological and remote sensing studies. Global Ecology and Biogeography 12, 191–205. Abuelgasim, A.A., Fernandas, R.A., Leblanc, S.G., 2006. Evaluation of national and global LAI products derived from optical remote sensing instruments over Canada. IEEE Transactions on Geoscience and Remote Sensing 44, 1872–1883. Baret, F., Clevers, J.G.P.W., Steven, M.D., 1995. The robustness of canopy gap fraction estimations from red and near-infrared reflectances: a comparison of approaches. Remote Sensing of Environment 54, 141–151.

247

Baret, F., Guyot, G., Major, D., 1989. TSAVI: a vegetation index which minimizes soil brightness effects on LAI and APAR estimation. In: Proceedings of 12th Canadian Symposium on Remote Sensing and IGARSS’89, Vancouver, Canada, pp. 1355–1358. Bounoua, L., Masek, J., Tourre, Y.M., 2006. Sensitivity of surface climate to land surface parameters: a case study using the simple biosphere model SiB2. Journal of Geophysical Research: Atmosphere 111, D22S09. Bunting, P.J., Lucas, R.M., 2006. The delineation of tree crowns in Australian mixed species forests using hyperspectral Compact Airborne Spectrographic Imager (CASI) data. Remote Sensing of Environment 101, 230–248. Butson, C., King, D.J., 1999. Semivariance analysis of forest structure and remote sensing data to determine an optimal sample plot size. In: Proceedings of 4th International Airborne Remote Sensing Conference and Exhibition and 21st Canadian Symposium on Remote Sensing, Ottawa, Canada, pp. 155–162. Bowker, D.E., Davis, R.E., Myrick, D.L., Stacy, K., Jones, W.T., 1985. Spectral refectances of natural targets for use in remote sensing studies. NASA Reference Publication 1139, 188. Bréda, N.J.J., 2003. Ground-based measurements of leaf area index: a review of methods, instruments and current controversies. Journal of Experimental Botany 54, 2403–2417. Campbell, C, 1986. Extinction coefficients for radiation in plant canopies calculated using an ellipsoidal inclination angle distribution. Agriculture and Forest Meteorology 36, 317–321. Carlson, T.N., Ripley, D.A., 1997. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sensing of Environment 62, 241– 252. Chavez Jr., P.S., 1996. Image-based atmospheric corrections: revisited and improved. Photogrammetric Engineering and Remote Sensing 62, 1025–1036. Chen, J.M., Black, T.A., 1992. Defining leaf area index for non-flat leaves. Plant, Cell & Environment 15, 421–429. Chen, J.M., Liu, J., Leblanc, S.G., Lacaze, R., Roujean, J.-L., 2003. Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption. Remote Sensing of Environment 84, 516–525. Chen, J.M., Menges, C.H., Leblanc, S.G., 2005. Global mapping of foliage clumping index using multi-angular satellite data. Remote Sensing of Environment 97, 447–457. Choudhury, B.J., Ahmed, N.U., Idso, S.B., Reginato, R.J., Daughtry, C.S.T., 1994. Relations between evaporation coefficients and vegetation indices studied by model simulations. Remote Sensing of Environment 50, 1–17. Combal, B., Baret, F., Weiss, M., Trubuil, A., Macé, D., Pragnère, A., Myneni, R., Knyazikhin, Y., Wang, L., 2002. Retrieval of canopy biophysical variables from bidirectional reflectance using prior information to solve the ill-posed inverse problem. Remote Sensing of Environment 84, 1–15. Daniel, T.W., Helms, J.A., Baker, F.S., 1979. Principles of Silviculture, 2nd edition. McGraw-Hill, USA. Demarez, V., Duthoit, S., Baret, F., Weiss, M., Dedieu, G., 2008. Estimation of leaf area and clumping indexes of crops with hemispherical photographs. Agricultural and Forest Meteorology 148, 644–655. Fang, H., Liang, S., 2005. A hybrid inversion method for mapping leaf area index from MODIS data: experiments and application to broadleaf and needleleaf canopies. Remote Sensing of Environment 94, 405–424. Ganguly, S., Samanta, A., Schull, M.A., Shabanov, N.V., Milesi, C., Nemani, R., Knyazikhin, Y., Myneni, R., 2008. Generating vegetation leaf area index Earth system data record from multiple sensors. Part 2: Implementation, analysis and validation. Remote Sensing of Environment 112, 4318– 4332. Garrigues, S., Shabanov, N.V., Swanson, K., Morisette, J.T., Baret, F., Myneni, R., 2008. Intercomparison and sensitivity analysis of Leaf Area Index retrievals from LAI-2000. AccuPAR, and digital hemispherical photography over croplands. Agricultural and Forest Meteorology 148, 1193–1209. GCOS, 2006. Satellite-Based Products for Climate, Supplemental Details to the Satellite-Based Component of the Implementation Plan for the Global Observing System for Climate in Support of the UNFCCC, GCOS-107 (WMO/TD No. 1338), p. 90. GCOS, 2008. Future Climate Change Research and Observations: Global Climate Observing System (GCOS), World Climate Research Programme (WCRP) and International Geosphere-Biosphere Programme (IGBP) Learning from the Intergovernmental Panel on Climate Change (IPCC). Fourth Assessment Report—Workshop and Survey Report, Sydney, GCOS-117, p. 68. Gillies, R.R., Carlson, T.N., 1995. Thermal remote sensing of surface soil water content with partial vegetation cover for incorporation into climate models. Journal of Applied Meteorology 34, 745–756. Gobron, N., Verstraete, M.N., 2008. ECV T11: Leaf Area Index (LAI). Assessment of the Status of the Development of Standards for the Terrestrial Essential Climate Variables. Global Terrestrial Observing System, Rome, p. 24. Gonsamo Gosa, A., 2009. Remote sensing of leaf area index: enhanced retrieval from close-range and remotely sensed optical observations. PhD Dissertation, Department of Geography, University of Helsinki, Helsinki, Finland. Report Series A 147, p. 192. Gonsamo, A., Pellikka, P., 2008. Methodology comparison for slope correction in canopy leaf area index estimation using hemispherical photography. Forest Ecology and Management 256, 749–759. Gonsamo, A., Pellikka, P., 2009. The computation of foliage clumping index using hemispherical photography. Agricultural and Forest Meteorology 149, 1781–1787.

248

A. Gonsamo / International Journal of Applied Earth Observation and Geoinformation 12 (2010) 233–248

Gonsamo, A., Pellikka, P., King, D.J., in press. Large scale leaf area index inversion algorithms from high resolution airborne imagery. International Journal of Remote Sensing. Gutman, G., Ignatov, A., 1998. The derivation of the green vegetation fraction from NOAA/AVHRR data for use in numerical weather prediction models. International Journal of Remote sensing 19, 1533–1543. Heinsch, F.A., Zhao, M., Running, S.W., Kimball, J.S., Nemani, R.R., Davis, K.J., Bolstad, P.V., Cook, B.D., Desai, A.R., Ricciuto, D.M., Lwa, B.E., Oechel, W.C., Kwon, H.J., Luo, H., Wofsy, S.C., Dunn, A.L., Munger, J.W., Baldocchi, D.D., Xu, L., Hollinger, D.Y., Richardson, A.D., Stoy, P.C., Siqueira, M.B.S., Monson, R.K., Burns, S.P., Flanagan, L.B., 2006. Evaluation of remote sensing based terrestrial productivity from MODIS using regional tower eddy flux network observations. IEEE Transactions on Geoscience and Remote Sensing 44, 1908–1925. Houborg, R., Boegh, E., 2008. Mapping leaf chlorophyll and leaf area index using inverse and forward canopy reflectance modeling and SPOT reflectance data. Remote Sensing of Environment 112, 186–202. Hu, Z., Islam, S., 1997. A frame work for analyzing and designing scale invariant remote sensing algorithms. IEEE Transactions on Geoscience and Remote Sensing 35, 747–755. Huete, A., 1988. A soil adjusted vegetation index (SAVI). Remote Sensing of Environment 25, 295–309. Huete, A.R., Jackson, R.D., Post, D.F., 1985. Spectral response of a plant canopy with different soil backgrounds. Remote Sensing of Environment 17, 37–53. Jacquemoud, S., Bacour, C., Poilvé, H., Frange, J.-P., 2000. Comparison of four radiative transfer models to simulate plant canopies reflectance: direct and inverse mode. Remote Sensing of Environment 74, 471–481. Jacquemoud, S., Baret, F., 1990. PROSPECT: a model of leaf optical properties spectra. Remote Sensing of Environment 34, 75–91. Jacquemoud, S., Verhoef, W., Baret, F., Bacour, C., Zarco-Tejada, P.J., Asner, G.P., Franc¸ois, C., Ustin, S.L., 2009. PROSPECT + SAIL models: a review of use for vegetation characterization. Remote Sensing of Environment 113, S56–S66. Jensen, J.R., 1996. Introductory Digital Image Processing: A Remote Sensing Perspective. Prentice-Hall, Englewood Cliffs, NJ. Jiang, Z., Huete, A.R., Chen, J., Chen, Y., Li, J., Yan, G., Zhang, X., 2006. Analysis of NDVI and scaled difference vegetation index retrievals of vegetation fraction. Remote Sensing of Environment 101, 366–378. Lacaze, R., Chen, J.M., Roujean, J.-L., Leblanc, S.G., 2002. Retrieval of vegetation clumping index using hot spot signatures measured by POLDER instrument. Remote Sensing of Environment 79, 84–95. Lang, A.R.G., Xiang, Y., 1986. Estimation of leaf area index from transmission of direct sunlight in discontinuous canopies. Agricultural and Forest Meteorology 35, 229–243. Leblanc, S.G., Chen, J.M., White, H.P., Latifvic, R., 2005. Canada-wide foliage clumping index mapping from multi-angular POLDER measurements. Canadian Journal of Remote Sensing 31, 364–376. Moran, M.S., Jackson, R.D., Slater, P.N., Teillet, P.M., 1992. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output. Remote Sensing of Environment 41, 169–184. Morisette, J.T., Baret, F., Privette, J.L., Myneni, R.B., Nickeson, J.E., Garrigues, S., Shabanov, N.V., Weiss, M., Fernandes, R., Leblanc, S.G., Kalacska, M., SanchezAzofeifa, G.A., Chubey, M., Rivard, B., Stenberg, P., Rautiainen, M., Voipio, P., Manninen, T., Pilant, A.N., Lewis, T.E., Iiames, J.S., Colombo, R., Meroni, M., Busetto, L., Cohen, W.B., Turner, D.P., Warner, E.D., Petersen, G.W., Seufert, G., Cook, R., 2006. Validation of global moderate-resolution LAI products: a framework proposed within the CEOS land product validation subgroup. IEEE Transactions on Geoscience and Remote Sensing 44, 1804–1817. Myneni, R.B., Hall, F.G., Sellers, P.J., Marshak, A.L., 1995. The interpretation of spectral vegetation indexes. IEEE Transactions on Geoscience and Remote Sensing 33, 481–486. Myneni, R.B., Hoffman, S., Knyazikhin, Y., Privette, J.L., Glassy, J., Tian, Y., Wang, Y., Song, X., Zhang, Y., Smith, G.R., 2002. Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment 83, 214–231. National Capital Commission (NCC), 2005. Gatineau Park Master Plan., http://www.ncc-ccn.gc.ca/data/2/rec docs/1768 Master Plan e.pdf. Nilson, T., 1971. A theoretical analysis of the frequency of gaps in plant stands. Agriculture and Forest Meteorology 8, 25–38.

Pisek, J., Chen, J.M., 2007. Comparison and validation of MODIS and VEGETATION global LAI products over four BigFoot sites in North America. Remote Sensing of Environment 109, 81–94. Price, J., 1992. Estimating vegetation amount from visible and near-infrared reflectance. Remote Sensing of Environment 41, 29–34. Qi, J., Chehbouni, A., Huete, A., Kerr, Y., Sorooshian, S., 1994. A modified soil adjusted vegetation index. Remote Sensing of Environment 48, 119–126. Richardson, A.J., Wiegand, C.L., 1977. Distinguishing vegetation from soil background information. Photogrammetric Engineering and Remote Sensing 43, 1541–1552. Rouse, J.W., Haas, R.H., Schell, J.A., Deering, D.W., Harlan, J.C., 1974. Monitoring the Vernal Advancements and Retrogradation (Greenwave Effect) of Nature Vegetation. NASA/GSFC Type-III Final Report, Greenbelt, MD, USA, p. 164. Sellers, P.J., Los, S.O., Tucker, C.J., Justice, C.O., Dazlich, D.A., Collatz, D.J., Randall, D.A., 1996. A revised land surface parameterization (SiB2) for atmospheric GCMs. Part 2: The generation of global fields of terrestrial biophysical parameters from satellite data. Journal of Climate 9, 706–737. Song, C., Woodcock, C.E., Seto, K.C., Pax Lenney, M., Macomber, S.A., 2001. Classification and change detection using Landsat TM data: when and how to correct atmospheric effects? Remote Sensing of Environment 75, 230– 244. Sprintsin, M., Karnieli, A., Berliner, P., Rotenberg, E., Yakir, D., Cohen, S., 2007. The effect of spatial resolution on the accuracy of leaf area index estimation for a forest planted in the desert transition zone. Remote Sensing of Environment 109, 416–428. Stern, W.R., Donald, C.M., 1961. Relationship of Radiation, Leaf Area Index and Crop Growth-Rate. Nature 189, 597–598. Tang, S., Chen, J.M., Zhu, Q., Li, X., Chen, M., Sun, R., Zhou, Y., Deng, F., Xie, D., 2007. LAI inversion algorithm based on directional reflectance kernels. Journal of Environmental Management 85, 638–648. Verhoef, W., 1984. Light scattering by leaf layers with application to canopy reflectance modelling: the SAIL model. Remote Sensing of Environment 16, 125–141. Vermote, E.F., Tanré, D., Deuzé, J.L., Herman, M., Morcrette, J.J., 1997. Second simulation of the satellite signal in the solar spectrum, 6S: an overview. IEEE Transactions on Geoscience and Remote Sensing 35, 675–686. Vermote, E.F., Vermeulen, A., 1999. Atmospheric Correction Algorithm: SPECTRAL REFLECTANCES (MOD09). Algorithm Technical Background Document (ATBD)., http://modis.gsfc.nasa.gov/data/atbd/atbd mod08.pdf. Verstraete, M.M., Pinty, B., 1991. The potential contribution of satellite remote sensing to the understanding of arid lands processes. Vegetation 91, 59– 72. Vohland, M., Mader, S., Dorigo, W., 2010. Applying different inversion techniques to retrieve stand variables of summer barley with PROSPECT + SAIL. International Journal of Applied Earth Observation and Geoinformation 12, 71–80. Walter, J.-M., Fournier, R.A., Soudani, K., Meyer, E., 2003. Integrating clumping effects in forest canopy structure: an assessment through hemispherical photographs. Canadian Journal of Remote Sensing 29, 388–410. Wang, Y.J., Woodcock, C.E., Buermann, W., Stenberg, P., Voipio, P., Smolander, H., Hame, T., Tian, Y.H., Hu, J.N., Knyazikhin, Y., Myneni, R.B., 2004. Evaluation of the MODIS LAI algorithm at a coniferous forest site in Finland. Remote Sensing of Environment 91, 114–127. Weiss, M., Baret, F., Smith, G.J., Jonckheere, I., Coppin, P., 2004. Review of methods for in situ leaf area index (LAI) determination. Part II. Estimation of LAI, errors and sampling. Agricultural and Forest Meteorology 121, 37–53. Xiao, J., Moody, A., 2005. A comparison of methods for estimating fractional green vegetation cover within a desert-to-upland transition zone in central New Mexico, USA. Remote Sensing of Environment 98, 237–250. Yang, W.Z., Tan, B., Huang, D., Rautiainen, M., Shabanov, N.V., Wang, Y., Privette, J.L., Huemmrich, K.F., Fensholt, R., Sandholt, I., Weiss, M., Ahl, D.E., Gower, S.T., Nemani, R.R., Knyazikhin, Y., Myneni, R.B., 2006. MODIS leaf area index products: from validation to algorithm improvement. IEEE Transactions on Geoscience and Remote Sensing 44, 1885–1899. Zhou, X., Guan, H., Xie, H., Wilson, J.L., 2009. Analysis and optimization of NDVI definitions and areal fraction models in remote sensing of vegetation. International Journal of Remote Sensing 30, 721–751.