- Email: [email protected]
Optimal interpolation analysis of leaf area index using MODIS data
Remote Sensing of Environment 104 (2006) 283 – 296 www.elsevier.com/locate/rse
Optimal interpolation analysis of leaf area index using MODIS data Yingxin Gu ⁎, Stéphane Bélair, Jean-François Mahfouf, Godelieve Deblonde Meteorological Research Branch, Meteorological Service of Canada, Dorval, Québec, Canada Received 28 October 2005; received in revised form 28 April 2006; accepted 29 April 2006
Abstract A simple data analysis technique for vegetation leaf area index (LAI) using Moderate Resolution Imaging Spectroradiometer (MODIS) data is presented. The objective is to generate LAI data that is appropriate for numerical weather prediction. A series of techniques and procedures which includes data quality control, time-series data smoothing, and simple data analysis is applied. The LAI analysis is an optimal combination of the MODIS observations and derived climatology, depending on their associated errors σo and σc. The “best estimate” LAI is derived from a simple three-point smoothing technique combined with a selection of maximum LAI (after data quality control) values to ensure a higher quality. The LAI climatology is a time smoothed mean value of the “best estimate” LAI during the years of 2002–2004. The observation error is obtained by comparing the MODIS observed LAI with the “best estimate” of the LAI, and the climatological error is obtained by comparing the “best estimate” of LAI with the climatological LAI value. The LAI analysis is the result of a weighting between these two errors. Demonstration of the method described in this paper is presented for the 15-km grid of Meteorological Service of Canada (MSC)'s regional version of the numerical weather prediction model. The final LAI analyses have a relatively smooth temporal evolution, which makes them more appropriate for environmental prediction than the original MODIS LAI observation data. They are also more realistic than the LAI data currently used operationally at the MSC which is based on land-cover databases. © 2006 Elsevier Inc. All rights reserved. Keywords: Leaf area index; Data analysis; MODIS; Vegetation; Atmospheric and environmental models; Remote sensing
1. Introduction The representation of land surface processes is an important component of atmospheric and environmental models since it has a major impact on precipitation and on the evolution of the atmospheric boundary layer. To make sure that land surface processes have an optimal impact in these models, surface characteristics such as soil moisture, vegetation, and snow cover have to be specified as accurately as possible. In the current Meteorological Service of Canada (MSC) forecasting and assimilation systems, vegetation characteristics are specified from the United States Geological Survey (USGS) global land-cover characteristics (GLCC) database that is processed from numerous thematic maps (Loveland et al., 2000). The vegetation ⁎ Corresponding author. SAIC contractor to the USGS Center for Earth Resources Observation and Science (EROS), Sioux Falls, SD 57198, USA. Tel.: +1 605 594 6576; fax: +1 605 594 6529. E-mail address: [email protected] (Y. Gu). 0034-4257/$ - see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2006.04.021
property (LAI) is obtained from the USGS GLCC 1-km database, whereas the seasonal variation is specified from a pre-determined (arbitrary) look-up table (Bélair et al., 2003; Giard & Bazile, 2000) without using any observation data. Remote sensing has been successfully used to detect and monitor the coverage and density of global vegetation (e.g. Cohen et al., 2003; Hansen et al., 2002; Myneni et al., 2002; Zeng et al., 2002). Green leaves commonly have larger reflectances in the near infrared (NIR) than in the visible range, a fact that is used in the normalized difference vegetation index (NDVI), which is the normalized reflectance difference between the NIR and red bands (Rouse et al., 1974). Higher NDVI values are usually associated with greater density and greenness of the plant canopy. Several relationships have been proposed to link NDVI to leaf area index (LAI) which is a measure of the foliage density directly related to plant photosynthesis and evapotranspiration (Bonan, 1993). These relationships between NDVI and LAI were obtained either empirically through analysis of field measurements for various vegetation types (e.g. Nemani &
284
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Running, 1995; Shippert et al., 1995) or theoretically through radiative transfer models (Myneni et al., 1997). Many satellite sensors, such as AVHRR (Advanced Very High Resolution Radiometer), MODIS (Moderate Resolution Imaging Spectroradiometer), Landsat TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper Plus), and SPOT-VEGETATION, are currently used to derive vegetation products (NDVI and LAI). High-resolution MODIS vegetation products for the global scale are easily obtained with a latency of 3 to 5 days from the Earth Observing System (EOS) data gateway (see EOS Data Gateway in Web References) making these products suitable for the Canadian forecasting and assimilation systems. Our main objective in this study is to present a simple method to generate time varying vegetation LAI appropriate for environmental forecasting systems. This new method involves several operations upon MODIS data such as data quality control, time-series data smoothing, and a simple data analysis method. The performance of this new method is evaluated for a 15-km grid regional version of the Global Environmental Multiscale (GEM) model currently used for regional weather forecasting at MSC (Mailhot et al., 2006). 2. Data analysis strategy 2.1. MODIS LAI product The MODIS LAI product is defined as the one-sided green leaf area per unit ground area (Myneni et al., 2002; Privette et al., 2002). Its production algorithm is based on three-dimensional radiative transfer theory (Myneni et al., 2002) and is developed for inversion using a look-up-table approach (Knyazikhin et al., 1998a,b; Tian et al., 2004). When this main algorithm fails, a back-up method based on the relationship between NDVI and LAI (Knyazikhin et al., 1998a; Myneni et al., 1997) is employed (Myneni et al., 2002). The 1-km MODIS LAI products are generated over an 8-day compositing period (MOD15A2, version4) which is based on the maximum fraction of absorbed photosynthetically active radiation (FPAR) (Myneni et al., 2002). The products are available to the public and can be ordered through the EOS Data Gateway or directly via FTP through Data Pool at the Land Processes Distributed Active Archive Center. The LAI products are reprojected on a Sinusoidal 10° grid with 36 × 18 tiles spanning the globe (Myneni et al., 2002) and are distributed in HDF (Hierarchical Data Format)-EOS format. Extensive quality control information regarding cloud and data processing conditions is also included in the data. The accuracy of MODIS LAI (version 4) for needleleaf forest is within 50% (Wang et al., 2004) and within 30% for agricultural area (Tan et al., 2005) from ground base validations. 2.2. Optimal estimation of LAI Data assimilation techniques are widely used for initializing numerical weather prediction (NWP) models in national weather centers (Daley, 1997; Lorenc, 1986). The main purpose of data assimilation is to combine data from several sources of information to provide optimal initial conditions for modeling and
forecasting systems. A simple one-dimensional data assimilation (i.e., data analysis) method (Mahfouf, 1991) is used in this paper (as described in this section) to produce an analysis of LAI from the MODIS data. If we assume that three pieces of information on LAI are available at a given time t: (i) an observation LAIo derived from satellite with an observation error σo, (ii) a background value LAIb with an associated error σb (i.e., first-guess from a vegetation model), and (iii) a climatological value LAIc with an associated error σc, then the analysis LAIa is obtained by combining in an optimal manner these three values. This optimal combination can be found by minimizing a quadratic cost-function representing the departures from these three LAI values (i.e., least-square estimate): J ðLAIa Þ ¼
LAIa −LAIo 2 LAIa −LAIb 2 LAIa −LAIc 2 þ þ ro rb rc
ð1Þ
It can be shown (Lorenc, 1986) that, when the associated errors are Gaussian, the choice of J is such that the value of LAI at its minimum corresponds to the best linear unbiased estimate, for which the variance of analysis error is minimum. The optimal value of LAIa is simply found by setting the gradient of this costfunction to zero: LAIa −LAIo LAIa −LAIb LAIa −LAIc ¼ 0 ð2Þ jJ ¼ 2 þ þ r2o r2c r2b From which the analysis LAIa is given: LAIo LAIb LAIc 1 1 1 −1 þ 2 þ 2 þ þ LAIa ¼ r2o r2b r2c r2o rc rb
ð3Þ
The equation can be rewritten as: LAIa ¼ ao LAIo þ ab LAIb þ ac LAIc
ð4Þ
with ao ¼ ab ¼
r2b r2c r2o r2c þ r2b r2o þ r2c r2b r2o r2c
r2o r2c þ r2b r2o þ r2c r2b
ac ¼
r2o r2c
r2b r2o þ r2b r2o þ r2c r2b
ð5Þ
where αo, αb, and αc are the weighting factors for the observations, the background state, and the climatology. In numerical weather prediction, the background information is provided by a short-range forecast. Since in this study we do not have a prognostic model to evolve LAI in time, the background information is not used in this paper. Only the observational and climatological parts are considered, and the LAI analysis is given by: ð6Þ
LAIa ¼ ao LAIo þ ac LAIc where ao ¼
r2c ; r2o þ r2c
and
ac ¼
r2o r2o þ r2c
ð7Þ
Additional quality control is done in the data analysis process, based on the difference between the current observation LAIo and the climatological information LAIc. If both errors are uncorrelated and follow Gaussian statistics with standard deviations σo and σc,
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
the difference shouldpalso follow a Gaussian statistics with a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi standard deviation of r2o þ r2c . Therefore, when qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð8Þ jLAIo −LAIc j > 2 r2o þ r2c which means that the difference is two standard deviations away from the mean of the distribution, the observation data is believed to be suspicious and is not used in the analysis (i.e., the analysis is set to the climatological value). In summary, this simple technique for assimilating MODIS LAI data requires four elements of information: the observed and climatological LAI (i.e., LAIo and LAIc) as well as their error statistics σo and σc. The method used to specify these four elements is described in the two following subsections. 2.3. LAI observation and climatological data The LAI 1-km data (MOD15A2, version 4) for 2002 to 2004 is used in this study to generate both the LAI observation and climatology. The process is schematically shown in Fig. 1 and a description follows. For the LAI observation (LAIo), an initial quality control is performed to remove bad values due to the presence of cloud or
285
due to atmospheric effects (see FPAR, LAI User's Guide in Web References). The pixels with bad values are identified from the quality control (QC) information included in the MOD15A2 file and are not included in the averaging process to the 15-km GEM model grid on which the assimilation is done. These 15-km LAI observations (LAIo), in addition to being directly fed into the data analysis system, are also used to produce the climatology LAIc by following the steps described below. A first operation (done on the 15-km LAIo) remedies to deficiencies of the MODIS LAI algorithm over snow-covered regions, especially those covered by mid-to-high-latitude needleleaf forests (Tian et al., 2004; Zhang et al., 2004). Over these regions, MODIS LAI values are too small (near zero in fact), and do not capture the evergreen trees that exist in reality (Tian et al., 2004). To solve this problem and make sure that LAI observations used to derive the LAI climatology are realistic for every season and every region of North America, the LAI values over regions covered with evergreen vegetation (e.g., evergreen needleleaf trees, mixed shrubs, and mixed wood forests) are replaced by LAI values derived from the summer season over the same regions. The replaced winter LAI will be overestimated if the understory of evergreen vegetation is deciduous, as the deciduous understory contribute to LAI in summer, but not in winter LAI. The identification of the problematic
Fig. 1. Flowchart of MODIS LAI data processing.
286
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
regions (non-computed pixels of LAI product) is given by the MSC's GEM LAI database, and the corrective operation can be expressed as: LAIoV¼ MAX½LAIo ; LAIo
summer * ðVF1
þ VF2 þ VF3Þ
ð9Þ
where LAIo′ is the LAI values averaged on the 15-km grid (after quality control), corrected for mid- and high-latitude deficiencies due to snow; LAIo_summer represents LAIo for July 4–July 12 8-day period of each year (on the 15-km grid); and VF1, VF2, and VF3 are the fractional vegetation cover (from USGS land-cover database) for evergreen needleleaf trees, mixed shrubs, and mixed wood
forests, respectively. There are some deciduous forests within mixed wood forests, and their LAI will be zero during winter. Considering this factor, the winter LAI values for mixed wood forest and mixed shrubs regions will be overestimated in Eq. (9). The second step in the process of generating LAIc is to apply a temporal smoother to LAIo′. Indeed, LAI should have a relatively smooth temporal evolution that follows seasonal and annual cycles. Furthermore, since noise signal caused by clouds and poor atmospheric conditions is negatively biased (i.e., systematic underestimation of LAI) (Chen et al., 2004), an “upper LAI envelope” smoothing strategy is required. Following this, a simple three-point
Fig. 2. Examples of “upper envelope” three-point smoothing method for MODIS LAI time series data (for year 2004). “Smooth 1”, “Smooth 2”, and “Smooth 3” refer to the number of processes that were used for the smoother. (a) Forest (point “a” in Fig. 3), (b) crops (point “c” in Fig. 3).
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
287
Fig. 3. Climatological LAI (m2 m− 2) derived from MODIS data for the years 2002 to 2004.
smoothing technique is combined with a selection of maximum LAI values to provide what is considered here as the “best estimate” for LAI (i.e., LAIo_smooth): LAIo
smooth ðtÞ
¼ MAX½LAIo V; ð0:5*LAIo VðtÞ
ð10Þ
þ0:25*ðLAIo Vðt−1Þ þ LAIo Vðt þ 1ÞÞÞ where t− 1 represents the data at the previous time-step and t +1 represents the data at the following time-step. Note that, in order to remove the noise more efficiently, the smoothing formula (10) is repeated three times in LAIo_smooth processing in this study. An example of the LAI evolution obtained by using this threepoint “upper envelope” smoothing method is given in Fig. 2. It shows that this smoothing method works well when the perturbations (errors) are only for one data point. For a longer sequence of contaminated data (e.g., see data during month “11” in Fig. 2a), the noise will not be removed as effectively. The two final steps for the LAI climatology is the averaging for the 3 years (N = 3) available in the context of this study (i.e., 2002–2004): LAImean ðt Þ ¼
X 1 2002;2004 LAIo N year
smooth ðt; year Þ
ð11Þ
and a final temporal smoothing: LAIc ðtÞ ¼ 0:5*LAImean ðtÞ þ 0:25*ðLAImean ðt−1Þ þ LAImean ðt þ 1ÞÞ
ð12Þ
This method ensures that the climatology of LAI is smooth and of a higher quality. Maps of LAIc are shown in Fig. 3 for particular times during each of the four seasons. They indicate that the LAI climatology captures the significant changes that occur during the year, as well as the slow evolution of boreal forests, even though MODIS observations are near zero for these high-latitude regions during the winter. 2.4. Climatological and observation errors (σc and σo) The variance of the climatological LAI represents the yearly variation of the LAI “best estimate” and is simply defined as follows: X h 1 2002;2004 r2c ðt Þ ¼ LAIo smooth ðt; yearÞ−LAIc ðt Þ2 ð13Þ N year The observation error (σo), on the other hand, is given by the departure of the LAI observation (prior to the modification done for evergreen forest in the northern regions) with respect to the “best estimate” of LAI: X h 1 2002;2004 r2o ðt Þ ¼ LAIo ðt; yearÞ−LAIo smooth ðt; yearÞ2 ð14Þ N year where LAIo (t, year) represents the observation data that will be assimilated to produce LAI for the atmospheric model, i.e., the observation data for each year without blending and smoothing.
288
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Fig. 4. Observed (left panels) and analyzed (right panels) LAI (m2 m− 2) from MODIS.
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
289
Because an “upper envelope” method is used to smooth the LAI data (Eq. (10)), it is possible to have LAIo(t) = LAIo_smooth(t) which leads to σo(t) values of zero. For the few points where this happens, an empirical method is used to specify the observation error σo [i.e., σo(t) = 0.1 × LAIo(t)]. The final observation error is expressed according to:
Finally, the σc and σo are temporally smoothed by using a threepoint smoothing method (without the “upper-envelope” aspect though).
ro ðtÞ ¼ MAXfro ðtÞ; ð0:1 LAIo
ð15Þ
3.1. Comparison of LAI analyses and LAI MODIS observations
Since the error variances σo2 and σc2 appear in the denominator (can not be zero) of (3), the minima for σo and σc are set to 0.01 m2 m− 2.
The assimilation of MODIS LAI data is done for 2004 following the strategy described above. LAI maps from the
smooth ðtÞÞg
3. Results and discussion
Fig. 5. Upper panels: time series of observed LAI (for year 2004), analyzed LAI (for year 2004), climatological LAI. Lower panels: time series of observation error, climatological error, and the difference between the analysis and climatology (Laia_c). (a). Deciduous broadleaf trees in western West Virginia, (b) deciduous broadleaf trees in southern Missouri, (c) crops in Iowa, (d) crops in southern Saskatchewan, (e) deciduous shrubs in northern Manitoba, (f) evergreen needleleaf trees in Ontario, (g) long grasses in southern Kansas, and (h) short grasses in western Oklahoma. The locations of the above 8 points are shown in Fig. 3.
290
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Fig. 5 (continued ).
analysis and MODIS data (prior to blending step in Fig. 2) are provided in Fig. 4. The comparison is done for four particular days: one day in each season. For winter and spring LAI maps (Fig. 4a and b), LAI values from MODIS are very small over northern regions, due to the strong snow effects on the LAI retrievals. The LAI analysis, on the other hand, seems more realistic over these regions and is more consistent with the land-cover databases for the regions with evergreen needleleaf trees, mixed shrubs, and mixed wood forests. For summer LAI maps (Fig. 4c), the observed and analyzed LAI
are in better agreement as expected (i.e., less problems associated with snow). During fall (Fig. 4d), differences between the analyses and MODIS observations are seen over northern Canada. Because of contamination by clouds or other processes, there are also differences in northern Alabama and Georgia (see red circle in left panel of Fig. 4d). Over these regions, MODIS LAI values are near zero, which is apparently not correct as the region should be covered by vegetation at that time of the year. The LAI analysis (see red circle in right panel of Fig. 4d) seems more realistic, showing why LAI MODIS observations can not be directly used
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
291
Fig. 5 (continued ).
in environmental modeling applications; and why a data analysis technique such as the one proposed in this study is necessary. The LAI time evolution during 2004 is shown for several points in the upper panels of Fig. 5 for the observation, the climatology, and the analysis (see Fig. 3 for the locations of the time series). One should first note that LAI analyses exhibit a smooth temporal evolution in which the noise signal, present in the LAI MODIS observations and caused by cloud contamination and atmospheric variability (Chen et al., 2004; Myneni et al., 2002), is significantly reduced. Another feature worthy of mention is related to the differences that can be seen at the start and end of the growing season. At these times, the LAI analyses have more realistic seasonal changes compared to the climatological LAI
(e.g., Fig. 5a fall season, Fig. 5d crops growing season, and Fig. 5e late spring season), due to the impact of LAIo in the analysis. Indeed, as indicated in (6) and (7), the LAI analysis is an optimal combination of the MODIS observations and the climatology, depending on their associated errors σo and σc. Therefore, climatological and observation errors (σo and σc) play key roles in the data analysis process. For example, for large σo and small σc, the LAI analysis is mainly given by the climatology. Conversely, for small σo with large σc, the analysis is mainly determined by the observations. Time series of σo and σc at particular points are shown in the lower panels of Fig. 5, as well as the difference between the analysis and climatology (LAIa_c). The annual variation of the
292
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Fig. 5 (continued ).
climatological errors σc mainly depends on the type of land cover. For deciduous broad leaf trees and deciduous shrubs (e. g., Fig. 5a, b, and e), σc in spring and fall are higher than during the other seasons, indicating that the highest yearly variation of LAI is in the leaf growing and falling seasons. During summer, the climatological errors σc are small and thus the analysis LAIa values are mainly determined by the climatology. For crop regions (e.g., Fig. 5c and d), the greatest yearly variation of LAI is found during summer, with σc values that are three times larger than the corresponding σo. Therefore, the summer time LAI analysis at these locations is mainly determined by the observation LAI. In regions of boreal evergreen needleleaf forests (Fig. 5f), because of the strong snow effect discussed before, σo values are much higher in
the winter season and the LAIa are nearly totally determined by the LAIc values. For grass regions (Fig. 5g and h), high σo values occur during April and October, indicating that the highest yearly variation of LAI is in the warm season. One can notice that, although significantly reduced, some noise signal still exist in the LAI analysis time series (Fig. 5). This is related to the fact that the analysis is done in a “real-time” fashion (i.e., with no knowledge of future LAI observations). This is different from a “smoothing” technique which benefits from knowledge of the LAI values after the time of analysis. Using background information (i.e., from model forecasts) in the data analysis process would likely reduce the noise signal more efficiently. We are therefore considering the possibility of using a prognostic model to
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
evolve LAI in time. Moreover, the short time period (3 years) that was used to produce the climatological LAI data is probably not sufficient. A longer time period (at least 5 years) would be necessary for a better climatological LAI database. Such a climatology will be built for data analysis as more data becomes available. For some extreme and unnatural situations (i.e., fires), it will take some time, maybe a few years, before this effect will be represented appropriately in the analysis (through climatology). This is also a limitation of the approach presented here. It should be noted that the “upper envelop” smoothing method may introduce positive biases of LAI, for instance in situations for which surface reflectances are over-corrected for atmospheric conditions. Also, the method developed in this study is only applied for North America, but it should be general enough to be applied
293
over other regions of the world. Some regions, like those covered by rainforests, may require some modifications of the algorithms. 3.2. Comparison of analysis MODIS LAI data and LAI data currently used for GEM model Since the main objective of the assimilation strategy proposed in this study is to provide LAI analyses to the environmental forecasting systems, it is important to examine how different these new analyses are compared to what is currently used in operation at MSC. As shown in the upper panel of Fig. 6 for a particular summertime day (4 August 2004), the MODIS LAI analyses are found to be significantly different from the operational LAI. As expected, LAI fields used in the GEM operational model (GEM
Fig. 6. Upper panel: MODIS LAI analysis (m2 m− 2) and GEM operational LAI (m2 m− 2), valid for 4 August 2004. Lower panel: The LAI (m2 m− 2) map obtained from CCRS [Fernandes et al., 2003] over Canada, valid for the period of 1–10 August 2004.
294
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
LAI) are much smoother than those obtained by assimilating MODIS data. Indeed, the vegetation classification database used in GEM exhibits less spatial variability than the LAI analysis, especially in summertime when this analysis is strongly influenced by current MODIS LAI observation data. It can also be noticed that canopy density in GEM seems to be severely overestimated for large regions in Saskatchewan, British Columbia, and southeastern USA, while they seem to be underestimated in other regions such as upper Michigan. In order to verify our new MODIS LAI analysis, it is compared with a LAI map (lower panel of Fig. 6) from the Canada Centre for Remote Sensing (CCRS) (Fernandes et al., 2003), validated over
the Canadian territory (Abuelgasim et al., 2006). The comparison result illustrates that the MODIS LAI analysis is generally in good agreement with the CCRS LAI map. It should first be noted that both LAI products have similar horizontal variability, while GEM LAI is smooth in comparison. There are many regions where the MODIS LAI analysis is in better agreement with the CCRS map than the GEM operational LAI is. For instance, both MODIS LAI analysis and CCRS LAI map show low LAI values (1–2 m2 m− 2) in northern Saskatchewan, while the GEM operational LAI is on the order of 4–5 m2 m− 2 in the same region. Also in the southern boundary region between Newfoundland and Quebec, the LAI values from both the MODIS analysis and the CCRS LAI map
Fig. 7. Time series of analyzed LAI (m2 m− 2) and of GEM operational LAI for forest (point “a” in Fig. 6) and crops (point “b” in Fig. 6) regions, valid for 2004.
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
have low values (2–4 m2 m− 2), while the LAI values from the GEM operational map have high values (4–6 m2 m− 2) in the same region. The LAI time series for points over forests and crops (Fig. 7) reveal that the current representation of vegetation in GEM could be quite far from that observed. Over the forests point (upper panel of Fig. 7), large differences are found in the leaf growing and falling seasons. Over the crop point (lower panel of Fig. 7), the simple correspondence tables used to specify vegetation in GEM overestimates LAI by a factor of 2 to 3 during summertime. These differences will certainly have an influence on the surface energy and water budgets over a large portion of North America (and the world) and it is likely that this could have a non-negligible impact on numerical weather prediction. An evaluation of this impact will be done in the near future. 4. Conclusions A simple data analysis technique for MODIS LAI data is presented in this paper. MODIS LAI data is obtained from the EOS Data Gateway. Data quality control, time-series data smoothing, and a simple data analysis method are applied in this work. Demonstration of the method described in this paper is presented for the 15-km grid of MSC's regional version of the GEM model. The LAI analysis is an optimal combination of the MODIS observations and derived climatology, depending on their associated errors σo and σc. The LAI climatology is a time smoothed mean value of the “best estimate” LAI (Eq. (10)) during the years of 2002–2004, which is of high quality. The computed LAI climatology captures the significant changes that occur during the year, as well as the slow evolution of boreal forests, even though MODIS observations are near zero for these high-latitude regions during the winter. The climatological and observation errors (σo and σc) play key roles in the data analysis process. The observation error is obtained by comparing the MODIS observed LAI with the “best estimate” of the LAI, and the climatological error is obtained by comparing the “best estimate” of LAI with the climatological LAI value. The LAI analysis is the result of a weighting between these two errors. For large σo and small σc, the LAI analysis is mainly given by the climatology (e.g., LAI analyses for evergreen needleleaf trees in winter). Conversely, for small σo with large σc, the LAI analysis is mainly determined by the MODIS observation (e.g., LAI analyses for crops in summer). Results show the final LAI analyses are of higher quality, have smoother temporal evolution, and are more consistent with land-cover databases. There are significant differences between the LAI analysis and the fields currently used operationally in MSC's atmospheric models, and the LAI analysis data is more realistic and appropriate for environmental forecasting than either the old operational LAI data used in the MSC's GEM model or the original MODIS LAI observation data. These differences should have a notable influence on the surface energy and water budget over a large portion of North America, implying that important impacts could be generated in the numerical weather forecasts. Testing the potential improvement of using this new analysis LAI data for GEM model will be
295
described in a future study. Refining the LAI analysis over Canada by using CCRS LAI products will be done in the future. Acknowledgement The authors would like to acknowledge Richard Fernandes for providing the CCRS LAI map. References Abuelgasim, A. A., Fernandes, R. A., & Leblanc, S. G. (2006). Evaluation of national and global LAI products derived from optical remote sensing instruments over Canada. Accepted by IEEE TGARSS Special Issue on Land validation. Bélair, S., Crevier, L. -P., Mailhot, J., Bilodeau, B., & Delage, Y. (2003). Operational implementation of the ISBA land surface scheme in the Canadian regional weather forecast model: Part I. Warm season results. Journal of Hydrometeorology, 4, 352−370. Bonan, G. (1993). Importance of leaf area index and forest type when estimating photosynthesis in boreal forests. Remote Sensing of Environment, 43, 303−314. Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., & Eklundh, L. (2004). A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sensing of Environment, 91, 332−344. Cohen, W. B., Maierpserger, T. K., Gower, S. T., & Turner, D. P. (2003). An improved strategy for regression of biophysical variables and Landsat ETM+ data. Remote Sensing of Environment, 84, 561−571. Daley, R. (1997). Atmospheric data assimilation. Journal of the Meteorological Society of Japan, 75(1B), 319−329. EOS Data Gateway: http://edcimswww.cr.usgs.gov/pub/imswelcome/ Fernandes, R. A., Butson, C., Leblanc, S. G., & Latifovic, R. (2003). Landsat-5 and Landsat-7 ETM+ based accuracy assessment of leaf area index products for Canada derived from SPOT-4 VEGETATION data. Canadian Journal of Remote Sensing, 29(2), 241−258. FPAR, LAI User's Guide: http://cybele.bu.edu/modismisr/products/modis/ userguide.pdf Giard, D., & Bazile, E. (2000). Implementation of a new assimilation scheme for soil and surface variables in a global NWP model. Monthly Weather Review, 128, 997−1015. Hansen, M. C., DeFries, S., Townshend, J. R. G., Sohlberg, A., Dimiceli, C., & Carroll, M. (2002). Towards an operational MODIS continuous field of percent tree cover algorithm: Examples using AVHRR and MODIS data. Remote Sensing of Environment, 83, 303−319. Knyazikhin, Y., Martonchik, J. V, Myneni, R. B., Diner, D. J, & Running, S. W (1998). Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from MODIS and MISR data. Journal of Geophysical Research, 103, 32,257−32,275. Knyazikhin, Y., Martonchik, J. V., Diner, D. J., Myneni, R. B., Verstraete, M. M, Pinty, B., et al. (1998). Estimation of vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from atmospherecorrected MISR data. Journal of Geophysical Research, 103, 32,239−32,256. Lorenc, A. C. (1986). Analysis methods for numerical weather prediction. Quarterly Journal of the Royal Meteorological Society, 112, 1177−1194. Loveland, T. R., Reed, B. C., Brown, J. F., Ohlen, D. O., Zhu, Z., Yang, L., et al. (2000). Development of a global land cover characteristics database and IGBP DISCover from 1 km AVHRR data. International Journal of Remote Sensing, 21(6/7), 1303−1330. Mahfouf, J. F. (1991). Analysis of soil moisture from near-surface parameters: A feasibility study. Journal of Applied Meteorology, 30, 1534−1547. Mailhot, J., Bélair, S., Lefaivre, L., Bilodeau, B., Desgagné, M., Girard, C., Glazer, A., Leduc, A. -M., Methot, A., Patoine, A., Plante, A., Rahill, A., Robinson, T., Talbot, D., Tremblay, A., Vaillancourt, P., & Zadra, A. (2006). The 15-km version of the Canadian regional forecast system. Atmosphere-Ocean, 44(2), 133−149. Myneni, R. B., Hoffman, S., Knyazikhin, Y., Privette, J. L, Glassy, J., Tian, Y., et al. (2002). Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment, 83, 214−231.
296
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Myneni, R. B., Nemani, R. R., & Running, S. W. (1997). Estimation of global leaf area index and absorbed PAR using radiative transfer model. IEEE Transactions on Geoscience and Remote Sensing, 35, 1380−1393. Nemani, R. R., & Running, S. W. (1995). Satellite monitoring of global land cover changes and their impact on climate. Climatic Change, 31, 395−413. Privette, J. L., Myneni, B., Knyazikhin, Y., Mukelabai, M., Roberts, G., Tian, Y., et al. (2002). Early spatial and temporal validation of MODIS LAI product in Africa. Remote Sensing of Environment, 83, 232−243. Rouse, J. W., Jr., Haas, R. H., Deering, D. W., Schell, J. A., & Harlan, J. C. (1974). Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation. NASA/GSFC Type III Final Report, Greenbelt, MD (pp. 371). Shippert, M. M., Walker, D. A., Auerbach, N. A., & Lewis, B. E. (1995). Biomass and leaf-area index maps derived from SPOT images for Toolik Lake and Imnavait Creek areas, Alaska. Polar Record, 31, 147−154. Tan, B., Hu, J., Zhang, P., Huang, D., & Shabanov, N. (2005). Validation of Moderate Resolution Imaging Spectroradiometer leaf area index product in
croplands of Alpilles, France. Journal of Geophysical Research, 110, D01107, doi:10.1029/2004JD004860 Tian, Y., Dickinson, R. E., Zhou, L., Zeng, X., Dai, Y., Myneni, R. B., et al. (2004). Comparison of seasonal and spatial variations of LAI/FPAR from MODIS and common land model. Journal of Geophysical Research, 109 (D1), D01103, doi:10.1029/2003JD003777 Wang, Y., Woodcock, C. E., Buermann, W., Stenberg, P., Voipio, P., Smolander, H., et al. (2004). Evaluation of the MODIS LAI algorithm at a coniferous forest site in Finland. Remote Sensing of Environment, 91, 114−127. Zeng, X., Shaikh, M., Dai, Y., Dickinson, R. E., & Myneni, R. (2002). Coupling of the common land model to the NCAR community climate model. Journal of Climate, 15, 1832−1854. Zhang, P., Anderson, B., Barlow, M., Tan, B., & Myneni, R. B. (2004). Climate related vegetation characteristics derived from MODIS LAI and NDVI. Journal of Geophysical Research, 109, D20105, doi:10.1029/2004JD004720
Optimal interpolation analysis of leaf area index using MODIS data Yingxin Gu ⁎, Stéphane Bélair, Jean-François Mahfouf, Godelieve Deblonde Meteorological Research Branch, Meteorological Service of Canada, Dorval, Québec, Canada Received 28 October 2005; received in revised form 28 April 2006; accepted 29 April 2006
Abstract A simple data analysis technique for vegetation leaf area index (LAI) using Moderate Resolution Imaging Spectroradiometer (MODIS) data is presented. The objective is to generate LAI data that is appropriate for numerical weather prediction. A series of techniques and procedures which includes data quality control, time-series data smoothing, and simple data analysis is applied. The LAI analysis is an optimal combination of the MODIS observations and derived climatology, depending on their associated errors σo and σc. The “best estimate” LAI is derived from a simple three-point smoothing technique combined with a selection of maximum LAI (after data quality control) values to ensure a higher quality. The LAI climatology is a time smoothed mean value of the “best estimate” LAI during the years of 2002–2004. The observation error is obtained by comparing the MODIS observed LAI with the “best estimate” of the LAI, and the climatological error is obtained by comparing the “best estimate” of LAI with the climatological LAI value. The LAI analysis is the result of a weighting between these two errors. Demonstration of the method described in this paper is presented for the 15-km grid of Meteorological Service of Canada (MSC)'s regional version of the numerical weather prediction model. The final LAI analyses have a relatively smooth temporal evolution, which makes them more appropriate for environmental prediction than the original MODIS LAI observation data. They are also more realistic than the LAI data currently used operationally at the MSC which is based on land-cover databases. © 2006 Elsevier Inc. All rights reserved. Keywords: Leaf area index; Data analysis; MODIS; Vegetation; Atmospheric and environmental models; Remote sensing
1. Introduction The representation of land surface processes is an important component of atmospheric and environmental models since it has a major impact on precipitation and on the evolution of the atmospheric boundary layer. To make sure that land surface processes have an optimal impact in these models, surface characteristics such as soil moisture, vegetation, and snow cover have to be specified as accurately as possible. In the current Meteorological Service of Canada (MSC) forecasting and assimilation systems, vegetation characteristics are specified from the United States Geological Survey (USGS) global land-cover characteristics (GLCC) database that is processed from numerous thematic maps (Loveland et al., 2000). The vegetation ⁎ Corresponding author. SAIC contractor to the USGS Center for Earth Resources Observation and Science (EROS), Sioux Falls, SD 57198, USA. Tel.: +1 605 594 6576; fax: +1 605 594 6529. E-mail address: [email protected] (Y. Gu). 0034-4257/$ - see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2006.04.021
property (LAI) is obtained from the USGS GLCC 1-km database, whereas the seasonal variation is specified from a pre-determined (arbitrary) look-up table (Bélair et al., 2003; Giard & Bazile, 2000) without using any observation data. Remote sensing has been successfully used to detect and monitor the coverage and density of global vegetation (e.g. Cohen et al., 2003; Hansen et al., 2002; Myneni et al., 2002; Zeng et al., 2002). Green leaves commonly have larger reflectances in the near infrared (NIR) than in the visible range, a fact that is used in the normalized difference vegetation index (NDVI), which is the normalized reflectance difference between the NIR and red bands (Rouse et al., 1974). Higher NDVI values are usually associated with greater density and greenness of the plant canopy. Several relationships have been proposed to link NDVI to leaf area index (LAI) which is a measure of the foliage density directly related to plant photosynthesis and evapotranspiration (Bonan, 1993). These relationships between NDVI and LAI were obtained either empirically through analysis of field measurements for various vegetation types (e.g. Nemani &
284
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Running, 1995; Shippert et al., 1995) or theoretically through radiative transfer models (Myneni et al., 1997). Many satellite sensors, such as AVHRR (Advanced Very High Resolution Radiometer), MODIS (Moderate Resolution Imaging Spectroradiometer), Landsat TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper Plus), and SPOT-VEGETATION, are currently used to derive vegetation products (NDVI and LAI). High-resolution MODIS vegetation products for the global scale are easily obtained with a latency of 3 to 5 days from the Earth Observing System (EOS) data gateway (see EOS Data Gateway in Web References) making these products suitable for the Canadian forecasting and assimilation systems. Our main objective in this study is to present a simple method to generate time varying vegetation LAI appropriate for environmental forecasting systems. This new method involves several operations upon MODIS data such as data quality control, time-series data smoothing, and a simple data analysis method. The performance of this new method is evaluated for a 15-km grid regional version of the Global Environmental Multiscale (GEM) model currently used for regional weather forecasting at MSC (Mailhot et al., 2006). 2. Data analysis strategy 2.1. MODIS LAI product The MODIS LAI product is defined as the one-sided green leaf area per unit ground area (Myneni et al., 2002; Privette et al., 2002). Its production algorithm is based on three-dimensional radiative transfer theory (Myneni et al., 2002) and is developed for inversion using a look-up-table approach (Knyazikhin et al., 1998a,b; Tian et al., 2004). When this main algorithm fails, a back-up method based on the relationship between NDVI and LAI (Knyazikhin et al., 1998a; Myneni et al., 1997) is employed (Myneni et al., 2002). The 1-km MODIS LAI products are generated over an 8-day compositing period (MOD15A2, version4) which is based on the maximum fraction of absorbed photosynthetically active radiation (FPAR) (Myneni et al., 2002). The products are available to the public and can be ordered through the EOS Data Gateway or directly via FTP through Data Pool at the Land Processes Distributed Active Archive Center. The LAI products are reprojected on a Sinusoidal 10° grid with 36 × 18 tiles spanning the globe (Myneni et al., 2002) and are distributed in HDF (Hierarchical Data Format)-EOS format. Extensive quality control information regarding cloud and data processing conditions is also included in the data. The accuracy of MODIS LAI (version 4) for needleleaf forest is within 50% (Wang et al., 2004) and within 30% for agricultural area (Tan et al., 2005) from ground base validations. 2.2. Optimal estimation of LAI Data assimilation techniques are widely used for initializing numerical weather prediction (NWP) models in national weather centers (Daley, 1997; Lorenc, 1986). The main purpose of data assimilation is to combine data from several sources of information to provide optimal initial conditions for modeling and
forecasting systems. A simple one-dimensional data assimilation (i.e., data analysis) method (Mahfouf, 1991) is used in this paper (as described in this section) to produce an analysis of LAI from the MODIS data. If we assume that three pieces of information on LAI are available at a given time t: (i) an observation LAIo derived from satellite with an observation error σo, (ii) a background value LAIb with an associated error σb (i.e., first-guess from a vegetation model), and (iii) a climatological value LAIc with an associated error σc, then the analysis LAIa is obtained by combining in an optimal manner these three values. This optimal combination can be found by minimizing a quadratic cost-function representing the departures from these three LAI values (i.e., least-square estimate): J ðLAIa Þ ¼
LAIa −LAIo 2 LAIa −LAIb 2 LAIa −LAIc 2 þ þ ro rb rc
ð1Þ
It can be shown (Lorenc, 1986) that, when the associated errors are Gaussian, the choice of J is such that the value of LAI at its minimum corresponds to the best linear unbiased estimate, for which the variance of analysis error is minimum. The optimal value of LAIa is simply found by setting the gradient of this costfunction to zero: LAIa −LAIo LAIa −LAIb LAIa −LAIc ¼ 0 ð2Þ jJ ¼ 2 þ þ r2o r2c r2b From which the analysis LAIa is given: LAIo LAIb LAIc 1 1 1 −1 þ 2 þ 2 þ þ LAIa ¼ r2o r2b r2c r2o rc rb
ð3Þ
The equation can be rewritten as: LAIa ¼ ao LAIo þ ab LAIb þ ac LAIc
ð4Þ
with ao ¼ ab ¼
r2b r2c r2o r2c þ r2b r2o þ r2c r2b r2o r2c
r2o r2c þ r2b r2o þ r2c r2b
ac ¼
r2o r2c
r2b r2o þ r2b r2o þ r2c r2b
ð5Þ
where αo, αb, and αc are the weighting factors for the observations, the background state, and the climatology. In numerical weather prediction, the background information is provided by a short-range forecast. Since in this study we do not have a prognostic model to evolve LAI in time, the background information is not used in this paper. Only the observational and climatological parts are considered, and the LAI analysis is given by: ð6Þ
LAIa ¼ ao LAIo þ ac LAIc where ao ¼
r2c ; r2o þ r2c
and
ac ¼
r2o r2o þ r2c
ð7Þ
Additional quality control is done in the data analysis process, based on the difference between the current observation LAIo and the climatological information LAIc. If both errors are uncorrelated and follow Gaussian statistics with standard deviations σo and σc,
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
the difference shouldpalso follow a Gaussian statistics with a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi standard deviation of r2o þ r2c . Therefore, when qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð8Þ jLAIo −LAIc j > 2 r2o þ r2c which means that the difference is two standard deviations away from the mean of the distribution, the observation data is believed to be suspicious and is not used in the analysis (i.e., the analysis is set to the climatological value). In summary, this simple technique for assimilating MODIS LAI data requires four elements of information: the observed and climatological LAI (i.e., LAIo and LAIc) as well as their error statistics σo and σc. The method used to specify these four elements is described in the two following subsections. 2.3. LAI observation and climatological data The LAI 1-km data (MOD15A2, version 4) for 2002 to 2004 is used in this study to generate both the LAI observation and climatology. The process is schematically shown in Fig. 1 and a description follows. For the LAI observation (LAIo), an initial quality control is performed to remove bad values due to the presence of cloud or
285
due to atmospheric effects (see FPAR, LAI User's Guide in Web References). The pixels with bad values are identified from the quality control (QC) information included in the MOD15A2 file and are not included in the averaging process to the 15-km GEM model grid on which the assimilation is done. These 15-km LAI observations (LAIo), in addition to being directly fed into the data analysis system, are also used to produce the climatology LAIc by following the steps described below. A first operation (done on the 15-km LAIo) remedies to deficiencies of the MODIS LAI algorithm over snow-covered regions, especially those covered by mid-to-high-latitude needleleaf forests (Tian et al., 2004; Zhang et al., 2004). Over these regions, MODIS LAI values are too small (near zero in fact), and do not capture the evergreen trees that exist in reality (Tian et al., 2004). To solve this problem and make sure that LAI observations used to derive the LAI climatology are realistic for every season and every region of North America, the LAI values over regions covered with evergreen vegetation (e.g., evergreen needleleaf trees, mixed shrubs, and mixed wood forests) are replaced by LAI values derived from the summer season over the same regions. The replaced winter LAI will be overestimated if the understory of evergreen vegetation is deciduous, as the deciduous understory contribute to LAI in summer, but not in winter LAI. The identification of the problematic
Fig. 1. Flowchart of MODIS LAI data processing.
286
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
regions (non-computed pixels of LAI product) is given by the MSC's GEM LAI database, and the corrective operation can be expressed as: LAIoV¼ MAX½LAIo ; LAIo
summer * ðVF1
þ VF2 þ VF3Þ
ð9Þ
where LAIo′ is the LAI values averaged on the 15-km grid (after quality control), corrected for mid- and high-latitude deficiencies due to snow; LAIo_summer represents LAIo for July 4–July 12 8-day period of each year (on the 15-km grid); and VF1, VF2, and VF3 are the fractional vegetation cover (from USGS land-cover database) for evergreen needleleaf trees, mixed shrubs, and mixed wood
forests, respectively. There are some deciduous forests within mixed wood forests, and their LAI will be zero during winter. Considering this factor, the winter LAI values for mixed wood forest and mixed shrubs regions will be overestimated in Eq. (9). The second step in the process of generating LAIc is to apply a temporal smoother to LAIo′. Indeed, LAI should have a relatively smooth temporal evolution that follows seasonal and annual cycles. Furthermore, since noise signal caused by clouds and poor atmospheric conditions is negatively biased (i.e., systematic underestimation of LAI) (Chen et al., 2004), an “upper LAI envelope” smoothing strategy is required. Following this, a simple three-point
Fig. 2. Examples of “upper envelope” three-point smoothing method for MODIS LAI time series data (for year 2004). “Smooth 1”, “Smooth 2”, and “Smooth 3” refer to the number of processes that were used for the smoother. (a) Forest (point “a” in Fig. 3), (b) crops (point “c” in Fig. 3).
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
287
Fig. 3. Climatological LAI (m2 m− 2) derived from MODIS data for the years 2002 to 2004.
smoothing technique is combined with a selection of maximum LAI values to provide what is considered here as the “best estimate” for LAI (i.e., LAIo_smooth): LAIo
smooth ðtÞ
¼ MAX½LAIo V; ð0:5*LAIo VðtÞ
ð10Þ
þ0:25*ðLAIo Vðt−1Þ þ LAIo Vðt þ 1ÞÞÞ where t− 1 represents the data at the previous time-step and t +1 represents the data at the following time-step. Note that, in order to remove the noise more efficiently, the smoothing formula (10) is repeated three times in LAIo_smooth processing in this study. An example of the LAI evolution obtained by using this threepoint “upper envelope” smoothing method is given in Fig. 2. It shows that this smoothing method works well when the perturbations (errors) are only for one data point. For a longer sequence of contaminated data (e.g., see data during month “11” in Fig. 2a), the noise will not be removed as effectively. The two final steps for the LAI climatology is the averaging for the 3 years (N = 3) available in the context of this study (i.e., 2002–2004): LAImean ðt Þ ¼
X 1 2002;2004 LAIo N year
smooth ðt; year Þ
ð11Þ
and a final temporal smoothing: LAIc ðtÞ ¼ 0:5*LAImean ðtÞ þ 0:25*ðLAImean ðt−1Þ þ LAImean ðt þ 1ÞÞ
ð12Þ
This method ensures that the climatology of LAI is smooth and of a higher quality. Maps of LAIc are shown in Fig. 3 for particular times during each of the four seasons. They indicate that the LAI climatology captures the significant changes that occur during the year, as well as the slow evolution of boreal forests, even though MODIS observations are near zero for these high-latitude regions during the winter. 2.4. Climatological and observation errors (σc and σo) The variance of the climatological LAI represents the yearly variation of the LAI “best estimate” and is simply defined as follows: X h 1 2002;2004 r2c ðt Þ ¼ LAIo smooth ðt; yearÞ−LAIc ðt Þ2 ð13Þ N year The observation error (σo), on the other hand, is given by the departure of the LAI observation (prior to the modification done for evergreen forest in the northern regions) with respect to the “best estimate” of LAI: X h 1 2002;2004 r2o ðt Þ ¼ LAIo ðt; yearÞ−LAIo smooth ðt; yearÞ2 ð14Þ N year where LAIo (t, year) represents the observation data that will be assimilated to produce LAI for the atmospheric model, i.e., the observation data for each year without blending and smoothing.
288
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Fig. 4. Observed (left panels) and analyzed (right panels) LAI (m2 m− 2) from MODIS.
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
289
Because an “upper envelope” method is used to smooth the LAI data (Eq. (10)), it is possible to have LAIo(t) = LAIo_smooth(t) which leads to σo(t) values of zero. For the few points where this happens, an empirical method is used to specify the observation error σo [i.e., σo(t) = 0.1 × LAIo(t)]. The final observation error is expressed according to:
Finally, the σc and σo are temporally smoothed by using a threepoint smoothing method (without the “upper-envelope” aspect though).
ro ðtÞ ¼ MAXfro ðtÞ; ð0:1 LAIo
ð15Þ
3.1. Comparison of LAI analyses and LAI MODIS observations
Since the error variances σo2 and σc2 appear in the denominator (can not be zero) of (3), the minima for σo and σc are set to 0.01 m2 m− 2.
The assimilation of MODIS LAI data is done for 2004 following the strategy described above. LAI maps from the
smooth ðtÞÞg
3. Results and discussion
Fig. 5. Upper panels: time series of observed LAI (for year 2004), analyzed LAI (for year 2004), climatological LAI. Lower panels: time series of observation error, climatological error, and the difference between the analysis and climatology (Laia_c). (a). Deciduous broadleaf trees in western West Virginia, (b) deciduous broadleaf trees in southern Missouri, (c) crops in Iowa, (d) crops in southern Saskatchewan, (e) deciduous shrubs in northern Manitoba, (f) evergreen needleleaf trees in Ontario, (g) long grasses in southern Kansas, and (h) short grasses in western Oklahoma. The locations of the above 8 points are shown in Fig. 3.
290
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Fig. 5 (continued ).
analysis and MODIS data (prior to blending step in Fig. 2) are provided in Fig. 4. The comparison is done for four particular days: one day in each season. For winter and spring LAI maps (Fig. 4a and b), LAI values from MODIS are very small over northern regions, due to the strong snow effects on the LAI retrievals. The LAI analysis, on the other hand, seems more realistic over these regions and is more consistent with the land-cover databases for the regions with evergreen needleleaf trees, mixed shrubs, and mixed wood forests. For summer LAI maps (Fig. 4c), the observed and analyzed LAI
are in better agreement as expected (i.e., less problems associated with snow). During fall (Fig. 4d), differences between the analyses and MODIS observations are seen over northern Canada. Because of contamination by clouds or other processes, there are also differences in northern Alabama and Georgia (see red circle in left panel of Fig. 4d). Over these regions, MODIS LAI values are near zero, which is apparently not correct as the region should be covered by vegetation at that time of the year. The LAI analysis (see red circle in right panel of Fig. 4d) seems more realistic, showing why LAI MODIS observations can not be directly used
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
291
Fig. 5 (continued ).
in environmental modeling applications; and why a data analysis technique such as the one proposed in this study is necessary. The LAI time evolution during 2004 is shown for several points in the upper panels of Fig. 5 for the observation, the climatology, and the analysis (see Fig. 3 for the locations of the time series). One should first note that LAI analyses exhibit a smooth temporal evolution in which the noise signal, present in the LAI MODIS observations and caused by cloud contamination and atmospheric variability (Chen et al., 2004; Myneni et al., 2002), is significantly reduced. Another feature worthy of mention is related to the differences that can be seen at the start and end of the growing season. At these times, the LAI analyses have more realistic seasonal changes compared to the climatological LAI
(e.g., Fig. 5a fall season, Fig. 5d crops growing season, and Fig. 5e late spring season), due to the impact of LAIo in the analysis. Indeed, as indicated in (6) and (7), the LAI analysis is an optimal combination of the MODIS observations and the climatology, depending on their associated errors σo and σc. Therefore, climatological and observation errors (σo and σc) play key roles in the data analysis process. For example, for large σo and small σc, the LAI analysis is mainly given by the climatology. Conversely, for small σo with large σc, the analysis is mainly determined by the observations. Time series of σo and σc at particular points are shown in the lower panels of Fig. 5, as well as the difference between the analysis and climatology (LAIa_c). The annual variation of the
292
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Fig. 5 (continued ).
climatological errors σc mainly depends on the type of land cover. For deciduous broad leaf trees and deciduous shrubs (e. g., Fig. 5a, b, and e), σc in spring and fall are higher than during the other seasons, indicating that the highest yearly variation of LAI is in the leaf growing and falling seasons. During summer, the climatological errors σc are small and thus the analysis LAIa values are mainly determined by the climatology. For crop regions (e.g., Fig. 5c and d), the greatest yearly variation of LAI is found during summer, with σc values that are three times larger than the corresponding σo. Therefore, the summer time LAI analysis at these locations is mainly determined by the observation LAI. In regions of boreal evergreen needleleaf forests (Fig. 5f), because of the strong snow effect discussed before, σo values are much higher in
the winter season and the LAIa are nearly totally determined by the LAIc values. For grass regions (Fig. 5g and h), high σo values occur during April and October, indicating that the highest yearly variation of LAI is in the warm season. One can notice that, although significantly reduced, some noise signal still exist in the LAI analysis time series (Fig. 5). This is related to the fact that the analysis is done in a “real-time” fashion (i.e., with no knowledge of future LAI observations). This is different from a “smoothing” technique which benefits from knowledge of the LAI values after the time of analysis. Using background information (i.e., from model forecasts) in the data analysis process would likely reduce the noise signal more efficiently. We are therefore considering the possibility of using a prognostic model to
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
evolve LAI in time. Moreover, the short time period (3 years) that was used to produce the climatological LAI data is probably not sufficient. A longer time period (at least 5 years) would be necessary for a better climatological LAI database. Such a climatology will be built for data analysis as more data becomes available. For some extreme and unnatural situations (i.e., fires), it will take some time, maybe a few years, before this effect will be represented appropriately in the analysis (through climatology). This is also a limitation of the approach presented here. It should be noted that the “upper envelop” smoothing method may introduce positive biases of LAI, for instance in situations for which surface reflectances are over-corrected for atmospheric conditions. Also, the method developed in this study is only applied for North America, but it should be general enough to be applied
293
over other regions of the world. Some regions, like those covered by rainforests, may require some modifications of the algorithms. 3.2. Comparison of analysis MODIS LAI data and LAI data currently used for GEM model Since the main objective of the assimilation strategy proposed in this study is to provide LAI analyses to the environmental forecasting systems, it is important to examine how different these new analyses are compared to what is currently used in operation at MSC. As shown in the upper panel of Fig. 6 for a particular summertime day (4 August 2004), the MODIS LAI analyses are found to be significantly different from the operational LAI. As expected, LAI fields used in the GEM operational model (GEM
Fig. 6. Upper panel: MODIS LAI analysis (m2 m− 2) and GEM operational LAI (m2 m− 2), valid for 4 August 2004. Lower panel: The LAI (m2 m− 2) map obtained from CCRS [Fernandes et al., 2003] over Canada, valid for the period of 1–10 August 2004.
294
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
LAI) are much smoother than those obtained by assimilating MODIS data. Indeed, the vegetation classification database used in GEM exhibits less spatial variability than the LAI analysis, especially in summertime when this analysis is strongly influenced by current MODIS LAI observation data. It can also be noticed that canopy density in GEM seems to be severely overestimated for large regions in Saskatchewan, British Columbia, and southeastern USA, while they seem to be underestimated in other regions such as upper Michigan. In order to verify our new MODIS LAI analysis, it is compared with a LAI map (lower panel of Fig. 6) from the Canada Centre for Remote Sensing (CCRS) (Fernandes et al., 2003), validated over
the Canadian territory (Abuelgasim et al., 2006). The comparison result illustrates that the MODIS LAI analysis is generally in good agreement with the CCRS LAI map. It should first be noted that both LAI products have similar horizontal variability, while GEM LAI is smooth in comparison. There are many regions where the MODIS LAI analysis is in better agreement with the CCRS map than the GEM operational LAI is. For instance, both MODIS LAI analysis and CCRS LAI map show low LAI values (1–2 m2 m− 2) in northern Saskatchewan, while the GEM operational LAI is on the order of 4–5 m2 m− 2 in the same region. Also in the southern boundary region between Newfoundland and Quebec, the LAI values from both the MODIS analysis and the CCRS LAI map
Fig. 7. Time series of analyzed LAI (m2 m− 2) and of GEM operational LAI for forest (point “a” in Fig. 6) and crops (point “b” in Fig. 6) regions, valid for 2004.
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
have low values (2–4 m2 m− 2), while the LAI values from the GEM operational map have high values (4–6 m2 m− 2) in the same region. The LAI time series for points over forests and crops (Fig. 7) reveal that the current representation of vegetation in GEM could be quite far from that observed. Over the forests point (upper panel of Fig. 7), large differences are found in the leaf growing and falling seasons. Over the crop point (lower panel of Fig. 7), the simple correspondence tables used to specify vegetation in GEM overestimates LAI by a factor of 2 to 3 during summertime. These differences will certainly have an influence on the surface energy and water budgets over a large portion of North America (and the world) and it is likely that this could have a non-negligible impact on numerical weather prediction. An evaluation of this impact will be done in the near future. 4. Conclusions A simple data analysis technique for MODIS LAI data is presented in this paper. MODIS LAI data is obtained from the EOS Data Gateway. Data quality control, time-series data smoothing, and a simple data analysis method are applied in this work. Demonstration of the method described in this paper is presented for the 15-km grid of MSC's regional version of the GEM model. The LAI analysis is an optimal combination of the MODIS observations and derived climatology, depending on their associated errors σo and σc. The LAI climatology is a time smoothed mean value of the “best estimate” LAI (Eq. (10)) during the years of 2002–2004, which is of high quality. The computed LAI climatology captures the significant changes that occur during the year, as well as the slow evolution of boreal forests, even though MODIS observations are near zero for these high-latitude regions during the winter. The climatological and observation errors (σo and σc) play key roles in the data analysis process. The observation error is obtained by comparing the MODIS observed LAI with the “best estimate” of the LAI, and the climatological error is obtained by comparing the “best estimate” of LAI with the climatological LAI value. The LAI analysis is the result of a weighting between these two errors. For large σo and small σc, the LAI analysis is mainly given by the climatology (e.g., LAI analyses for evergreen needleleaf trees in winter). Conversely, for small σo with large σc, the LAI analysis is mainly determined by the MODIS observation (e.g., LAI analyses for crops in summer). Results show the final LAI analyses are of higher quality, have smoother temporal evolution, and are more consistent with land-cover databases. There are significant differences between the LAI analysis and the fields currently used operationally in MSC's atmospheric models, and the LAI analysis data is more realistic and appropriate for environmental forecasting than either the old operational LAI data used in the MSC's GEM model or the original MODIS LAI observation data. These differences should have a notable influence on the surface energy and water budget over a large portion of North America, implying that important impacts could be generated in the numerical weather forecasts. Testing the potential improvement of using this new analysis LAI data for GEM model will be
295
described in a future study. Refining the LAI analysis over Canada by using CCRS LAI products will be done in the future. Acknowledgement The authors would like to acknowledge Richard Fernandes for providing the CCRS LAI map. References Abuelgasim, A. A., Fernandes, R. A., & Leblanc, S. G. (2006). Evaluation of national and global LAI products derived from optical remote sensing instruments over Canada. Accepted by IEEE TGARSS Special Issue on Land validation. Bélair, S., Crevier, L. -P., Mailhot, J., Bilodeau, B., & Delage, Y. (2003). Operational implementation of the ISBA land surface scheme in the Canadian regional weather forecast model: Part I. Warm season results. Journal of Hydrometeorology, 4, 352−370. Bonan, G. (1993). Importance of leaf area index and forest type when estimating photosynthesis in boreal forests. Remote Sensing of Environment, 43, 303−314. Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., & Eklundh, L. (2004). A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sensing of Environment, 91, 332−344. Cohen, W. B., Maierpserger, T. K., Gower, S. T., & Turner, D. P. (2003). An improved strategy for regression of biophysical variables and Landsat ETM+ data. Remote Sensing of Environment, 84, 561−571. Daley, R. (1997). Atmospheric data assimilation. Journal of the Meteorological Society of Japan, 75(1B), 319−329. EOS Data Gateway: http://edcimswww.cr.usgs.gov/pub/imswelcome/ Fernandes, R. A., Butson, C., Leblanc, S. G., & Latifovic, R. (2003). Landsat-5 and Landsat-7 ETM+ based accuracy assessment of leaf area index products for Canada derived from SPOT-4 VEGETATION data. Canadian Journal of Remote Sensing, 29(2), 241−258. FPAR, LAI User's Guide: http://cybele.bu.edu/modismisr/products/modis/ userguide.pdf Giard, D., & Bazile, E. (2000). Implementation of a new assimilation scheme for soil and surface variables in a global NWP model. Monthly Weather Review, 128, 997−1015. Hansen, M. C., DeFries, S., Townshend, J. R. G., Sohlberg, A., Dimiceli, C., & Carroll, M. (2002). Towards an operational MODIS continuous field of percent tree cover algorithm: Examples using AVHRR and MODIS data. Remote Sensing of Environment, 83, 303−319. Knyazikhin, Y., Martonchik, J. V, Myneni, R. B., Diner, D. J, & Running, S. W (1998). Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from MODIS and MISR data. Journal of Geophysical Research, 103, 32,257−32,275. Knyazikhin, Y., Martonchik, J. V., Diner, D. J., Myneni, R. B., Verstraete, M. M, Pinty, B., et al. (1998). Estimation of vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from atmospherecorrected MISR data. Journal of Geophysical Research, 103, 32,239−32,256. Lorenc, A. C. (1986). Analysis methods for numerical weather prediction. Quarterly Journal of the Royal Meteorological Society, 112, 1177−1194. Loveland, T. R., Reed, B. C., Brown, J. F., Ohlen, D. O., Zhu, Z., Yang, L., et al. (2000). Development of a global land cover characteristics database and IGBP DISCover from 1 km AVHRR data. International Journal of Remote Sensing, 21(6/7), 1303−1330. Mahfouf, J. F. (1991). Analysis of soil moisture from near-surface parameters: A feasibility study. Journal of Applied Meteorology, 30, 1534−1547. Mailhot, J., Bélair, S., Lefaivre, L., Bilodeau, B., Desgagné, M., Girard, C., Glazer, A., Leduc, A. -M., Methot, A., Patoine, A., Plante, A., Rahill, A., Robinson, T., Talbot, D., Tremblay, A., Vaillancourt, P., & Zadra, A. (2006). The 15-km version of the Canadian regional forecast system. Atmosphere-Ocean, 44(2), 133−149. Myneni, R. B., Hoffman, S., Knyazikhin, Y., Privette, J. L, Glassy, J., Tian, Y., et al. (2002). Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment, 83, 214−231.
296
Y. Gu et al. / Remote Sensing of Environment 104 (2006) 283–296
Myneni, R. B., Nemani, R. R., & Running, S. W. (1997). Estimation of global leaf area index and absorbed PAR using radiative transfer model. IEEE Transactions on Geoscience and Remote Sensing, 35, 1380−1393. Nemani, R. R., & Running, S. W. (1995). Satellite monitoring of global land cover changes and their impact on climate. Climatic Change, 31, 395−413. Privette, J. L., Myneni, B., Knyazikhin, Y., Mukelabai, M., Roberts, G., Tian, Y., et al. (2002). Early spatial and temporal validation of MODIS LAI product in Africa. Remote Sensing of Environment, 83, 232−243. Rouse, J. W., Jr., Haas, R. H., Deering, D. W., Schell, J. A., & Harlan, J. C. (1974). Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation. NASA/GSFC Type III Final Report, Greenbelt, MD (pp. 371). Shippert, M. M., Walker, D. A., Auerbach, N. A., & Lewis, B. E. (1995). Biomass and leaf-area index maps derived from SPOT images for Toolik Lake and Imnavait Creek areas, Alaska. Polar Record, 31, 147−154. Tan, B., Hu, J., Zhang, P., Huang, D., & Shabanov, N. (2005). Validation of Moderate Resolution Imaging Spectroradiometer leaf area index product in
croplands of Alpilles, France. Journal of Geophysical Research, 110, D01107, doi:10.1029/2004JD004860 Tian, Y., Dickinson, R. E., Zhou, L., Zeng, X., Dai, Y., Myneni, R. B., et al. (2004). Comparison of seasonal and spatial variations of LAI/FPAR from MODIS and common land model. Journal of Geophysical Research, 109 (D1), D01103, doi:10.1029/2003JD003777 Wang, Y., Woodcock, C. E., Buermann, W., Stenberg, P., Voipio, P., Smolander, H., et al. (2004). Evaluation of the MODIS LAI algorithm at a coniferous forest site in Finland. Remote Sensing of Environment, 91, 114−127. Zeng, X., Shaikh, M., Dai, Y., Dickinson, R. E., & Myneni, R. (2002). Coupling of the common land model to the NCAR community climate model. Journal of Climate, 15, 1832−1854. Zhang, P., Anderson, B., Barlow, M., Tan, B., & Myneni, R. B. (2004). Climate related vegetation characteristics derived from MODIS LAI and NDVI. Journal of Geophysical Research, 109, D20105, doi:10.1029/2004JD004720