Spatial Enhancement of MODIS-based Images of Leaf Area Index: Application to the Boreal Forest Region of Northern Alberta, Canada

Abstract

Leaf area index (LAI) is one of the most commonly used ecological variables in describing forests. Since 2000, 1-km resolution Moderate Resolution Imaging Spectroradiometer (MODIS)-based 8-day composites of LAI have been operationally available from the National Aeronautics and Space Administration (NASA), USA, at no cost to the user. In this paper, we present a simple protocol to enhance the spatial resolution of NASA-produced LAI composites to 250-m resolution. This is done by fusing MODIS-based estimates of enhanced vegetation index (EVI), consisting of 16-day 250-m resolution composites (also from NASA), with estimates of LAI. We apply the protocol to derive 250-m resolution maps of LAI for the boreal forest region of northern Alberta, Canada. Data fusion was possible in this study because of the inherent linear correlation that exists between EVI and LAI for the April to October growing period of 2005–2008, producing r²-values of 0.85–0.95 and p-values < 0.0001. Comparison of MODIS-based LAI with field-based measurements using the Tracing Radiation and Architecture of Canopies (TRAC) sensor and LAI-2000 Plant Canopy Analyzer showed reasonable agreement across values; statistical comparison of LAI data points produced an r²-value of 0.71 and a p-value <0.0001. Seventy one percent of MODIS-based LAI were within ±20% of field estimates.

1. Introduction

Leaf area is an important factor in terrestrial ecosystems and represents the area involved in: (i) the plant-to-air and air-to-plant transfer of gases and particulate matter, (ii) energy transfer and exchange, (iii) eco-physiological processes of evapotranspiration, photosynthesis, and net primary production, and (iv) forest landscape growth and development. LAI is commonly defined as ½ of the total green area of forests, accounting for all sides of the leaves, per unit horizontal ground-surface area [1] and represents the vertical integration of leaf area along the length of forest canopies. LAI is dynamic and varies according to tree species composition, stand developmental stage, site conditions, time of year, and management regime [2].

Determination of spatiotemporal dynamics of LAI across large areas (e.g., covering hundreds to thousands of hectares) is only possible with remote-sensing-based techniques (e.g., [3,4,5,6]). These techniques typically rely on establishing empirical relationships between normalised difference vegetation index (NDVI) and LAI [7,8,9,10,11]. However, NDVI tends to saturate over dense vegetation (LAI > 3–4; [12]) or during the height of the growing season [13], despite LAI being dynamic during this time. Since 2000, Moderate Resolution Imaging Spectroradiometer (MODIS)-based 8-day composites of LAI have been made available to the scientific community at 1-km resolution by the National Aeronautics and Space Administration (NASA), USA, at no cost to the user. These products are primarily generated by solving a 3-dimensional radiative-transfer model using atmospherically-corrected MODIS spectral surface reflectance and biome identification [3,14,15]. In understanding ecosystem dynamics at large spatial scales, the scientific community would greatly benefit from simple techniques that could be used to improve the spatial resolution of existing products of LAI.

Enhanced vegetation index (EVI, at 250-m resolution [4]) brings with it two clear advantages over NDVI concerning the depiction of LAI, namely: (i) EVI is sensitive to variations in LAI (NDVI, in contrast, responds more to plant greenness [4]); and (ii) EVI, unlike NDVI, does not saturate at peak foliage concentrations [4,13,16]. Preliminary results from Huete et al. [4] reveal the presence of a linear relationship between EVI and LAI for a number of biomes. Here we attempt to confirm Huete et al.’s observation of linearity for the boreal forest of northern Alberta, Canada, for the April–October period of 2005–2008. Evidence of a strong linear correlation between EVI and LAI would provide basis for augmenting coarse-resolution images of LAI with data fusion methods.

The data fusion technique by Hassan et al. [17], which we implement in this study, was first developed to enhance images of growing degree days (GDD; a temperature-related index developed from MODIS-based land surface temperature [13,17]) from 1-km to 250-m resolution using images of EVI (at 250-m resolution) as basis for fusion. Data fusion was possible in this study because of the strong linear correlation between GDD and EVI (r² = 0.87; [13]). The method by Hassan et al. [17] differs from other data fusion techniques due to its dependence on the relationship between ecologically-correlated variables rather than on spectral similarities between reflectance bands [e.g., 18,19]. Following the development of spatially-enhanced images of LAI for the Athabasca Oil Sands Region of northern Alberta (Figure 1) for June 2009, we validate the capabilities of Hassan et al.’s method [17] by comparing point-estimates of LAI extracted from the enhanced images with corresponding field-based estimates.

Figure 1. Extent of study area. A landcover map, derived from an annual composite of 2004 MODIS images, appears in the background. The black polygon outlines the Province of Alberta. The red box identifies the Athabasca Oil Sands Region, an area where field-based estimates of LAI were acquired to validate estimates of LAI generated from MODIS data.

Figure 1. Extent of study area. A landcover map, derived from an annual composite of 2004 MODIS images, appears in the background. The black polygon outlines the Province of Alberta. The red box identifies the Athabasca Oil Sands Region, an area where field-based estimates of LAI were acquired to validate estimates of LAI generated from MODIS data.

2. Methods

2.1. General Description of the Study Area and Data Requirements

Figure 1 shows the extent of the study area and landcover for northern Alberta (MOD12Q1 v. 4; annual 1-km resolution composite of 2004 MODIS images of landcover). The study area extends to the northern boundary of the Province of Alberta in the north, ~200 km inside the Province of Saskatchewan to the east, just south of the City of Edmonton to the south, and ~25 km inside the Province of British Columbia to the west (Figure 1). The majority of the study area is characterised by (i) a boreal plains, forest-dominated landscape consisting of 71.4% needleleaf and 0.9% broadleaf forests, and (ii) a cold-dry climate [20]. Tree species common to the area include trembling aspen (Populus tremuloides), balsam poplar (Populus balsamifera), white spruce (Picea glauca), black spruce (Picea mariana), and jack pine (Pinus banksiana). In poorly drained areas (i.e., fens and bogs), tamarack (Larix laricina) and black spruce often dominate. For Fort McMurray, near the centre of the Athabasca Oil Sands Region (56°39.0'N latitude, 111°13.2'W longitude; Figure 1), the mean annual temperature and annual precipitation for the 1971-2000 normal period are 0.7 °C and 455.5 mm, respectively. The January and July mean temperatures for the same period are −18.8 °C and 16.8 °C [21].

MODIS-based products used in this study included (i) fifty two 8-day and twenty six 16-day composites of LAI and EVI (i.e., MOD15A2 v. 5 and MOD13Q1 v. 5 images at 1-km and 250-m resolution, respectively), for the April–October period of 2005–2008 for assessing the presence of correlation between LAI and EVI, and (ii) a second set of eight LAI and four EVI composites for June 2009 to validate spatially-enhanced estimates of LAI (250-m resolution). The first set of images was collected over the April–October period because (i) this period coincided with the growing season, and (ii) LAI varied most over this time [20]. To validate MODIS-based LAI for June 2009 after image enhancement, field-based LAI for June 2009 were acquired at 75 sites across the Athabasca Oil Sands Region of northern Alberta (Figure 1; see Section 2.4).

2.2. Image Pre-Processing

LAI and EVI images for the 2005–2008 and June 2009 periods were re-projected to UTM (NAD83), Zone 12. In order to generate images for the entire study area, it required that two adjacent scenes of both LAI and EVI images be mosiacked.

Based on a preliminary cell-by-cell evaluation of average LAI (at 1-km resolution) and EVI (at 250 m) for the growing season of 2006 (Figure 2), it was observed that acceptable correlation existed within 80% of the data for needleleaf forests (r² = 0.57).

Figure 2. Scatterplot of cell-to-cell comparison of growing season (2006) averages for LAI and EVI for needleleaf forests; F-statistic = 2,120.31 and p < 0.0001.

Figure 2. Scatterplot of cell-to-cell comparison of growing season (2006) averages for LAI and EVI for needleleaf forests; F-statistic = 2,120.31 and p < 0.0001.

High noise to signal ratio in the data (Figure 2) was most probably associated with (i) inconsistency in temporal and spatial resolutions in the MODIS-based products used in this study (8–16 day averages over 250-m to 1-km resolution), and (ii) spatial variation in atmospheric transmissivity and other extrinsic factors that affect LAI calculations based on satellite data. To reinforce the signal between LAI and EVI, values were averaged across the entire study area (irrespective of biome) for each 8-day period. Similar averaging schemes have been used by other researchers, e.g., [16], to reduce variability inherent in satellite data. Temporal behaviour of 2005-2008 mean values was then used to determine the degree of correlation between “study area” LAI and EVI (see Section 3, for more details).

2.3. Data Fusion

Data fusion, after Hassan et al. [17], involves the generation of an artificial image (AI), which takes into account statistical properties of EVI within a 3 cell × 3 cell moving window, i.e.,:

$$AI=\frac{EV{I}_{ins}}{EV{I}_{mean}},$$

(1)

where EVI_ins is the instantaneous point-value of EVI and EVI_mean is the mean of all EVI values inside the moving window. In theory, AI can be viewed as an index that describes the relation of an instantaneous value of EVI to the mean value of EVI of surrounding pixels, and functions as a weight in the calculation of LAI at 250-m resolution, i.e.,:

$$LA{I}_{250m}=AI\text{ x }LA{I}_{1km},$$

(2)

where LAI_{1 km} represents the initial NASA-produced MODIS image of LAI at 1-km resolution.

2.4. Field-based LAI Data Acquisition and Processing

Validation of MODIS-based LAI was based on a June-2009 field-collection of LAI at 75 sites across the Athabasca Oil Sands Region of northern Alberta (Figure 1). To address variation in light conditions during the measurement period, LAI was determined with both the Tracing Radiation and Architecture of Canopies (TRAC) sensor and LAI-2000 Plant Canopy Analyzer because of their differences in lighting requirements. The TRAC is optimised for direct light and the LAI-2000, for diffused light conditions. Final selection of LAI from the two sensors is based on recorded sky conditions and the time when measurements were taken. For partially cloudy skies in late morning (09:30 AM LST) to mid-afternoon (05:30 PM LST), an average LAI is calculated from the individual LAI measurements determined with the TRAC and LAI-2000 sensors. At low sun elevations and increased diffused light conditions during early morning (sunrise to 09:30 AM LST) and late afternoon (05:30 PM LST to sundown), final LAI is based on the LAI-2000 sensor reading, regardless of sky conditions.

In general, light measurements were taken along transects of 250 m, unless site conditions (e.g., changes in stand composition, physical barriers) prevented the use of long transects. Both sensors required that light measurements be taken in open light conditions (in fields, alongside roads) prior to measuring light conditions in tall vegetation to define a light reference use in the calculation of LAI. Calculation of LAI is performed after each measurement session with the assistance of sensor-specific processing software available with the purchase of the equipment.

3. Results and Discussion

Figure 3 shows the temporal behaviour of LAI and EVI for the April–October period of 2005–2008. It revealed that both LAI and EVI followed near parallel tracks with a distinctive peak each year at around 20 July [day of year (DOY): 201; Figure 3a–e]. Polynomial fits to both LAI and EVI revealed strong correlations with DOY, i.e., for LAI, r²-values ranged from 0.92–0.95 for individual years and for EVI, r²-values ranged from 0.97–0.99; for all cases considered, p-values were <0.0001 (see Table 1, for additional statistics).

Figure 3. Temporal behaviour of LAI and EVI for the April–October period (DOY 97–297) of (a) 2005, (b) 2006, (c) 2007 and (d) 2008, and for (e) all years (2005–2008).

Figure 3. Temporal behaviour of LAI and EVI for the April–October period (DOY 97–297) of (a) 2005, (b) 2006, (c) 2007 and (d) 2008, and for (e) all years (2005–2008).

Table 1. Polynomial fit (i.e., y = a₀x⁴ − a₁x³ + a₂x² − a₃x + a₄) and associated statistics; x and y denote DOY and either LAI or EVI. In all cases considered, p-values were <0.0001.

Variable	Year	a0	a1	a2	a3	a4	r2	F-stat
LAI	2005	3e-08	2e-05	0.0065	0.70	26.82	0.93	70.03
EVI		3e-09	3e-06	0.0007	0.08	3.40	0.99	189.24
LAI	2006	3e-08	2e-05	0.0063	0.67	25.50	0. 92	59.93
EVI		2e-09	2e-06	0.0004	0.04	1.53	0. 98	149.15
LAI	2007	4e-08	3e-05	0.0075	0.81	31.00	0. 95	100.82
EVI		3e-09	2e-06	0.0006	0.06	2.59	0. 98	122.07
LAI	2008	4e-08	3e-05	0.0084	0.94	36.71	0. 94	76.98
EVI		3e-09	3e-06	0.0007	0.08	3.60	0. 98	110.84
LAI	All years	3e-08	3e-05	0.0072	0.78	30.00	0. 92	279.71
EVI		3e-09	2e-06	0.0006	0.07	2.78	0. 97	351.15

Table 1. Polynomial fit (i.e., y = a₀x⁴ − a₁x³ + a₂x² − a₃x + a₄) and associated statistics; x and y denote DOY and either LAI or EVI. In all cases considered, p-values were <0.0001.

Variable	Year	a0	a1	a2	a3	a4	r2	F-stat
LAI	2005	3e-08	2e-05	0.0065	0.70	26.82	0.93	70.03
EVI		3e-09	3e-06	0.0007	0.08	3.40	0.99	189.24
LAI	2006	3e-08	2e-05	0.0063	0.67	25.50	0. 92	59.93
EVI		2e-09	2e-06	0.0004	0.04	1.53	0. 98	149.15
LAI	2007	4e-08	3e-05	0.0075	0.81	31.00	0. 95	100.82
EVI		3e-09	2e-06	0.0006	0.06	2.59	0. 98	122.07
LAI	2008	4e-08	3e-05	0.0084	0.94	36.71	0. 94	76.98
EVI		3e-09	3e-06	0.0007	0.08	3.60	0. 98	110.84
LAI	All years	3e-08	3e-05	0.0072	0.78	30.00	0. 92	279.71
EVI		3e-09	2e-06	0.0006	0.07	2.78	0. 97	351.15

Figure 4 shows the relation between the “study area” (Figure 1) mean values of LAI and EVI. Strong linear relationships are observed, producing r²-values ranging from 0.85–0.95 and p-values <0.0001. Results for northern Alberta confirm Huete et al.’s [4] observation of linearity in other biomes. However, more investigation is needed for widespread application across a greater number of biomes globally.

Figure 4. Data comparison of LAI and EVI for the April–October period (DOY 97–297) of (a) 2005 (F-statistics = 87.69), (b) 2006 (F-statistics = 225.08), (c) 2007 (F-statistics = 63.68) and (d) 2008 (F-statistics = 119.65), and for (e) all years (2005–2008; F-statistics = 416.98). For all cases considered, the p-values were <0.0001.

Figure 4. Data comparison of LAI and EVI for the April–October period (DOY 97–297) of (a) 2005 (F-statistics = 87.69), (b) 2006 (F-statistics = 225.08), (c) 2007 (F-statistics = 63.68) and (d) 2008 (F-statistics = 119.65), and for (e) all years (2005–2008; F-statistics = 416.98). For all cases considered, the p-values were <0.0001.

In this study, we assume that the 1-km resolution images of LAI have been previously validated by NASA [3,14,15]. As a result, examination of the level of agreement between mean values of LAI at 250-m resolution and the original values at 1-km resolution serves as an appropriate means of validation (Figure 5). For each of the landcover types identified in Figure 1, Figure 5 provides some sample calculations for six arbitrary 1-km pixels. It revealed that the mean LAI at 250-m resolution compared to the original LAI at 1-km resolution were nearly similar for two landcover categories, i.e., (i) 3.6 vs. 3.7 for needleleaf forests [Figure 5a], and (ii) 4.6 vs. 4.7 for grasses/cereal crops [Figure 5b], and exactly the same for the remaining four categories, i.e., (i) 1.7 for shrubs [Figure 5c]; (ii) 3.0 for broadleaf crops [Figure 5d]; (iii) 2.7 for savannah [Figure 5e]; and (iv) 5.4 for broadleaf forests [Figure 5f]. A study-area-wide comparison of the original 1-km and enhanced 250-m image of LAI for an 8-day period in June 2009 (DOY 169–176), further illustrates the extent the two products match at 1-km resolution (r² = 0.87, with a slope = 1.02; Figure 6). Eighty percent of the data points are distributed along the 1:1 line (blue symbols in Figure 6). Mean absolute difference between LAI at 1-km and 250-m resolution (Figure 6) is roughly 0.59.

Figure 5. Comparison of LAI values at 1-km resolution without data fusion (left panel) and at 250-m resolution with data fusion (right panel) for various landcover types identified in Figure 1, including (a) needleleaf forests (comprising 71.4% of the study area), (b) grasses/cereal crops (15.3%), (c) shrubs (3.7%), (d) broadleaf crops (1.8%), (e) savannah (1.6%), and (f) broadleaf forests (0.9%).

Figure 5. Comparison of LAI values at 1-km resolution without data fusion (left panel) and at 250-m resolution with data fusion (right panel) for various landcover types identified in Figure 1, including (a) needleleaf forests (comprising 71.4% of the study area), (b) grasses/cereal crops (15.3%), (c) shrubs (3.7%), (d) broadleaf crops (1.8%), (e) savannah (1.6%), and (f) broadleaf forests (0.9%).

Figure 6. Study-area-wide comparison of LAI from the original 1-km and enhanced 250-m image for an 8-day period in June 2009 (DOY 169–176); F-statistic=10164 and p < 0.0001.

Figure 6. Study-area-wide comparison of LAI from the original 1-km and enhanced 250-m image for an 8-day period in June 2009 (DOY 169–176); F-statistic=10164 and p < 0.0001.

Figure 7 provides an example display of the spatial distribution of LAI at 250-m resolution for an 8-day period in June 2009 (DOY 169–176) based on the application of data fusion to the Athabasca Oil Sands Region (Figure 1). Greatest LAI (5–8) is observed to occur in the river valley, where an abundance of hardwood trees are found. Lowest LAI (0–1) occurs on a major plain formation (east of the road network; Figure 7), characterised by coarse sandy soils and a predominance of low-productivity jack pine stands in the dryer areas and an abundance of wetlands (fens and bogs) in local depressions [22].

Figure 7. Spatial distribution of LAI for an 8-day period in June 2009 (DOY 169–176) at 250-m resolution for the Athabasca Oil Sands Region (red box in Figure 1). Open circles in the centre, along the road network (green lines) and elsewhere in the image, denote the sites where field-based measurements of LAI were acquired.

Figure 7. Spatial distribution of LAI for an 8-day period in June 2009 (DOY 169–176) at 250-m resolution for the Athabasca Oil Sands Region (red box in Figure 1). Open circles in the centre, along the road network (green lines) and elsewhere in the image, denote the sites where field-based measurements of LAI were acquired.

To provide further validation of Hassan et al.’s method [17], we have also compared the LAI-values acquired from the field with those developed by data fusion (Figure 8). It revealed the presence of a strong linear relationship between estimates for 72 of the field sites (i.e., r² = 0.71 and slope = 0.70; Figure 8). LAI for three sites near open water were removed from the analysis, because at 1-km resolution the MODIS sensor was unable to discriminate the LAI-site characteristics from the adjacent water body—a feature that could not be corrected with data fusion. Of the 72 MODIS-based LAI used in the analysis 71% of the values fell within ±20% of the field-based values. Twenty percent deviation between satellite-based and field-based estimates is common when ecological systems are concerned [23,24]. This discrepancy is most likely a characteristic of the fact:

(i)

Field-based measurements are collected at a particular time of a day, whereas MODIS-based estimates are generated from 8- and 16-day composites, and

(ii)

Field-based measurements are collected by sensing the canopy from below, and as a result the measure of LAI cannot account for the understorey component. In contrast, RS-based methods measure LAI from sensing the top of the canopy. Under sub-optimal overstorey densities, understorey canopy effects are most likely to influence sensor readings and subsequent calculation of LAI.

Figure 8. Comparison of 72 field-based measurements with their corresponding MODIS-based LAI for June 2009 at 250-m resolution; F-statistic = 171.23 and p-value < 0.0001.

Figure 8. Comparison of 72 field-based measurements with their corresponding MODIS-based LAI for June 2009 at 250-m resolution; F-statistic = 171.23 and p-value < 0.0001.

4. Concluding Remarks

In this paper we demonstrate a relatively strong linear correlation between LAI and EVI. Because of this linear relationship, we can use an existing data fusion technique to enhance existing LAI images (at 1-km resolution) with EVI (at 250-m resolution). Reasonable agreement was obtained when point-estimates from spatially-enhanced LAI images for June 2009 were compared with their field-based equivalents collected at the locations where the field measurements were taken. Having the ability to enhance existing LAI products is critical to our understanding of earth-system processes at moderate spatial resolution over large areas.