Coupled retrieval of aerosol optical thickness, columnar water vapor and surface reflectance maps from ENVISAT/MERIS data over land

An algorithm for the derivation of atmospheric parameters and surface reflectance data from MEdium Resolution Imaging Specrometer Instrument (MERIS) on board ENVIronmental SATellite (ENVISAT) images has been developed. Geo-rectified aerosol optical thickness (AOT), columnar water vapor (CWV) and spectral surface reflectance maps are generated from MERIS Level-1b data over land. The algorithm has been implemented so that AOT, CWV and reflectance products are provided on an operational manner, making no use of ancillary parameters apart from those attached to MERIS products. For this reason, it has been named Self-Contained Atmospheric Parameters Estimation from MERIS data (SCAPE-M). The fundamental basis of the algorithm and applicable error figures are presented in the first part of this paper. In particular, errors of ± 0.03, ± 4% and ± 8% have been estimated for AOT, CWV and surface reflectance retrievals, respectively, by means of a sensitivity analysis based on a synthetic data set simulated under a usual MERIS scene configuration over land targets. The assumption of a fixed aerosol model, the coarse spatial resolution of the AOT product and the neglection of surface reflectance directional effects were also identified as limitations of SCAPE-M. Validation results are detailed in the second part of the paper. Comparison of SCAPE-M AOT retrievals with data from AErosol RObotic NETwork (AERONET) stations showed an average Root Mean Square Error (RMSE) of 0.05, and an average correlation coefficient R² of about 0.7-0.8. R² values grew up to more than 0.9 in the case of CWV after comparison with the same stations. A good correlation is also found with the MERIS Level-2 ESA CWV product. Retrieved surface reflectance maps have been successfully compared with reflectance data derived from the Compact High Resolution Imaging Spectrometer (CHRIS) on board the PRoject for On-Board Autonomy (PROBA) in the first place. Reflectance retrievals have also been compared with reflectance data derived from MERIS images by the Bremen AErosol Retrieval (BAER) method. A good correlation in the red and near-infrared bands was found, although a considerably higher proportion of pixels was successfully processed by SCAPE-M.     

An algorithm for the retrieval of 30-m snow-free albedo from Landsat surface reflectance and MODIS BRDF

We present a new methodology to generate 30-m resolution land surface albedo using Landsat surface reflectance and anisotropy information from concurrent MODIS 500-m observations. Albedo information at fine spatial resolution is particularly useful for quantifying climate impacts associated with land use change and ecosystem disturbance. The derived white-sky and black-sky spectral albedos may be used to estimate actual spectral albedos by taking into account the proportion of direct and diffuse solar radiation arriving at the ground. A further spectral-to-broadband conversion based on extensive radiative transfer simulations is applied to produce the broadband albedos at visible, near infrared, and shortwave regimes. The accuracy of this approach has been evaluated using 270 Landsat scenes covering six field stations supported by the SURFace RADiation Budget Network (SURFRAD) and Atmospheric Radiation Measurement Southern Great Plains (ARM/SGP) network. Comparison with field measurements shows that Landsat 30-m snow-free shortwave albedos from all seasons generally achieve an absolute accuracy of ±0.02-0.05 for these validation sites during available clear days in 2003-2005, with a root mean square error less than 0.03 and a bias less than 0.02. This level of accuracy has been regarded as sufficient for driving global and regional climate models. The Landsat-based retrievals have also been compared to the operational 16-day MODIS albedo produced every 8-days from MODIS on Terra and Aqua (MCD43A). The Landsat albedo provides more detailed landscape texture, and achieves better agreement (correlation and dynamic range) with in-situ data at the validation stations, particularly when the stations include a heterogeneous mix of surface covers.     

Atmospheric correction of optical imagery from MODIS and Reanalysis atmospheric products

In this paper we analyze the differences obtained in the atmospheric correction of optical imagery covering bands located in the Visible and Near Infra-Red (VNIR), Short-Wave Infra-Red (SWIR) and Themal-Infrared (TIR) spectral regions when atmospheric profiles extracted from different sources are used. In particular, three sensors were used, Compact High Resolution Imaging Spectrometer (CHRIS), Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER) and Landsat5 Thematic Mapper (TM), whereas four atmospheric profiles sources were considered: i) local soundings launched near the sensor overpass time, ii) Moderate Resolution Radiometer (MODIS) atmospheric profiles product (MOD07), iii) Atmospheric Correction Parameter Calculator (ACPC) generated by the National Center for Environmental Prediction (NCEP) and iv) Modified Atmospheric Profiles from Reanalysis Information (MAPRI), which includes data from NCEP and National Center of Atmospheric Research (NCAR) Reanalysis project but interpolated to 34 atmospheric levels and resampled to 0.5° × 0.5°. MODIS aerosol product (MOD04) was also used to extract Aerosol Optical Thickness (AOT) values at 550 nm. Analysis was performed for three test dates (12th July 2003, 18th July 2004 and 13th July 2005) over an agricultural area in Spain. Results showed that air temperature vertical profiles were similar for the four sources, whereas dew point temperature profiles showed significant differences at some particular levels. Atmospheric profiles were used as input to MODTRAN4 radiative transfer code in order to compute atmospheric parameters involved in atmospheric correction, with the aim of retrieving surface reflectances in the case of VNIR and SWIR regions, and Land Surface Temperature (LST) in the case of the TIR region. For the VNIR and SWIR region, significant differences depending on the atmospheric profile used were not found, particularly in the Visible region in which the AOT content is the main parameter involved in the atmospheric correction. In the case of TIR, differences depending on the atmospheric profile used were appreciable, since in this case the main parameter involved in the atmospheric correction is the water vapor content, which depends on the vertical profile. In terms of LST retrieval from ASTER data (2004 test case), all profiles provided satisfactory results compared to the ones obtained when using a local sounding, with errors of 0.3 K for ACPC and MAPRI cases and 0.7 K for MOD07. When retrieving LST from TM data (2005 test case), errors for MOD07 and MAPRI were 0.6 and 0.9 K respectively, whereas ACPC provided an error of 2 K. The results presented in this paper show that the different atmospheric profile sources are useful for accurate atmospheric correction when local soundings are not available. In particular, MOD07 product provides atmospheric information at the highest spatial resolution, 5 km, although its use is limited from 2000 to present, whereas MAPRI provides historical information from 1970 to present, but at lower spatial resolution.     

Comparison of atmospheric aerosol climatologies over southwestern Spain derived from AERONET and MODIS

The southwestern area of Spain, by its geographical and climatological conditions, is a key location for the characterization of atmospheric aerosol properties. The present study is aimed at evaluating the reliability of satellite-based aerosol climatologies, as inferred from level 2 standard aerosol products such as the Terra-MODIS (Moderate Resolution Imaging Radiometer) MOD04 aerosol product, with an application over this region during the period 2000-2008.This evaluation is carried out by means of comparison with ground-based data from the AERONET station of El Arenosillo (Spain, 37.1N, 6.7W), which has been providing continuous data since 2000. The focus of this paper is the climatology of two aerosol optical parameters: the aerosol optical depth (AOD) and the Ångström exponent.AERONET ground-based measurements give an annual mean value of 0.16 ± 0.12 and a median of 0.12 for the AOD, and a mean value of 1.20 ± 0.47 for the Ångström exponent. The seasonal pattern is characterized by two maxima, the most important maximum occurs in summer months, and the other one in late-winter/early-spring. Lowest values appear in fall and winter, however, a local minimum is observed in July which is only detected with the long-term data series.The mean climatological AOD based on AERONET exhibits complex seasonal patterns (i.e. with multiple local extrema), which are not always captured by MODIS-based climatology. MODIS only reproduces low values of the AOD in winter and high values in summer, as well as the local minimum of July which is sharper when using over-land retrievals. The time series of the AOD retrieved from MODIS both over land and ocean are in relatively good agreement with the ground-based measurements, with a monthly overestimation of about 30% on average, and higher differences in spring. Seasonal patterns from MODIS are better reproduced over land than over ocean. The agreement between daily AERONET and MODIS, as assessed by linear regression, gives correlation coefficients above 80% and an intercept bias below 0.03.     

Assessing the coupling between surface albedo derived from MODIS and the fraction of diffuse skylight over spatially-characterized landscapes

In this effort, the MODerate Resolution Imaging Spectroradiometer (MODIS) (Collection V005) Bidirectional Reflectance Distribution Function (BRDF)/Albedo algorithm is used to retrieve instantaneous surface albedo at a point in time and under specific atmospheric conditions. These retrievals are then used to study the role that the fraction of diffuse skylight plays under realistic scenarios of anisotropic diffuse illumination and multiple scattering between the surface and atmosphere. Simulations of the sky radiance using the MODTRAN®5.1 radiative transfer model were performed under different aerosol optical properties, illumination conditions, and surface characteristics to describe these effects on surface albedo retrievals from MODIS. This technique was examined using a validation scheme over four measurement sites with varied aerosol levels and landscapes, ranging from croplands to tundra ecosystems, and over extended time periods. Furthermore, a series of geostatistical analyses were performed to examine the types of spatial patterns observed at each measurement site. In particular, Enhanced Thematic Mapper Plus (ETM+) retrievals of surface albedo were acquired to analyze the change in variogram model parameters as a function of increased window-size. Results were then used to assess the degree to which a given point measurement is able to capture the intrinsic variability at the scale of MODIS observations. Assessments of MODIS instantaneous albedos that account for anisotropic multiple scattering, over snow-free and snow-covered lands and at all diurnal solar zenith angles, show a slight improvement over the albedo formulations that treat the downwelling diffuse radiation as isotropic. Comparisons with field measurements show biases improving by 0.004-0.013 absolute units (root-mean-squared error) or 0.1%-2.0% relative error.     

A hybrid inversion method for mapping leaf area index from MODIS data: experiments and application to broadleaf and needleleaf canopies

Leaf area index (LAI) is an important variable needed by various land surface process models. It has been produced operationally from the Moderate Resolution Imaging Spectroradiometer (MODIS) data using a look-up table (LUT) method, but the inversion accuracy still needs significant improvements. We propose an alternative method in this study that integrates both the radiative transfer (RT) simulation and nonparametric regression methods. Two nonparametric regression methods (i.e., the neural network [NN] and the projection pursuit regression [PPR]) were examined. An integrated database was constructed from radiative transfer simulations tuned for two broad biome categories (broadleaf and needleleaf vegetations). A new soil reflectance index (SRI) and analytically simulated leaf optical properties were used in the parameterization process. This algorithm was tested in two sites, one at Maryland, USA, a middle latitude temperate agricultural area, and the other at Canada, a boreal forest site, and LAI was accurately estimated. The derived LAI maps were also compared with those from MODIS science team and ETM+ data. The MODIS standard LAI products were found consistent with our results for broadleaf crops, needleleaf forest, and other cover types, but overestimated broadleaf forest by 2.0-3.0 due to the complex biome types.     

Separating surface emissivity and temperature using two-channel spectral indices and emissivity composites and comparison with a vegetation fraction method

The temperature-independent thermal infrared spectral indices (TISI) method is employed for the separation of land surface temperature (LST) and emissivity from surface radiances (atmospherically corrected satellite data). The daytime reflected solar irradiance and the surface emission at ∼3.8 μm have comparable magnitudes. Using surface radiances and a combination of day-night 2-channel TISI ratios, the ∼3.8 μm reflectivity is derived. For implementing the TISI method, coefficients for NOAA 9-16 AVHRR channels are obtained. A numerical analysis with simulated surface radiances shows that for most surface types (showing nearly Lambertian behavior) the achievable accuracy is ∼0.005 for emissivity (AVHRR channel-5) and ∼1.5 K for LST. Data from the European Centre for Medium-Range Weather Forecasts (ECMWF) is used for calculation of atmospheric attenuation. Comparisons are made over a part of central Europe on two different dates (seasons). Clouds pose a major problem to surface observations; hence, monthly emissivity composites are derived. Additionally, using TISI-based monthly composites of emissivities, a normalized difference vegetation index (NDVI)-based method is tuned to the particular study area and the results are intercompared. Once the coefficients are known, the NDVI method is easily implemented but holds well only for vegetated areas. The error of the NDVI-based emissivities (with respect to the TISI results) ranges between −0.038 and 0.032, but for vegetated areas the peak of the error-histogram is at ∼0.002. The algorithm for retrieving emissivity via TISI was validated with synthetic data. Due to the different spatial scales of satellite and surface measurements and the lack of homogeneous areas, which are representative for low-resolution pixels and ground measurements, ground-validation is a daunting task. However, for operational products ground-truth validation is necessary. Therefore, also an approach to identify suitable validation sites for meteorological satellite products in Europe is described.     

Models for surface reflection of radiance and polarized radiance: Comparison with airborne multi-angle photopolarimetric measurements and implications for modeling top-of-atmosphere measurements

In this paper, we investigate the surface-atmosphere radiative interaction in application to the problem of aerosol satellite remote sensing over land. First, we test different models of the Bidirectional Reflectance and Polarization Distribution Function (BRDF and BPDF) for bare soil and vegetation surfaces using multi-angle, multi-spectral photopolarimetric airborne measurements of the Research Scanning Polarimeter (RSP). Then, we investigate the performance of different models of BRDF and BPDF for modeling top-of-atmosphere measurements. We have found that different BRDF models can describe the RSP measurements equally well. However, for soil surfaces, the different BRDF models show a different dependence on illumination geometry (solar zenith and azimuth angles), as well as a different dependence on viewing angle outside the range of RSP measurements. This implies that different models describe the surface-atmosphere interaction differently, leading for soil surfaces to differences in the top-of-atmosphere reflectance up to 4-5%, whereas at surface level the models agree within 2% for RSP illumination and measurement geometry. For vegetation, the different BRDF models show more similar dependence on illumination geometry, meaning that, in general, the differences in top-of-atmosphere reflectances are smaller than the differences in surface total reflectances. For the BPDF, we compare the empirical model of Nadal and Breon (1999) and the model developed by Maignan et al. (2009) with a newly developed model. The latter model compares better with RSP measurements. It was shown that, though all models have essentially different angular profiles at different illumination and viewing geometries, the difference of the top-of-atmosphere degree of linear polarization is less or is of the same order as the degree of linear polarization difference at the surface level taken at RSP illumination and measurement geometry. For the considered models, it can be up to 0.015 but is mostly below 0.005.     

Evaluation of MODIS land surface temperature data to estimate air temperature in different ecosystems over Africa

The estimation of near surface air temperature (Ta) is useful for a wide range of applications such as agriculture, climate related diseases and climate change studies. Air temperature is commonly obtained from synoptic measurements in weather stations. In Africa, the spatial distribution of weather stations is often limited and the dissemination of temperature data is variable, therefore limiting their use for real-time applications. Compensation for this paucity of information may be obtained by using satellite-based methods. However, the derivation of near surface air temperature (Ta), from the land surface temperature (Ts) derived from satellite is far from straight forward. Some studies have tried to derive maximum Ta from satellites through regression analysis but the accuracy obtained is quite variable according to the study. The main objective of this study was to explore the possibility of retrieving high-resolution Ta data from the Moderate Resolution Imaging Spectroradiometer (MODIS) Ts products over different ecosystems in Africa. First, comparisons between night MODIS Ts data with minimum Ta showed that MODIS nighttime products provide a good estimation of minimum Ta over different ecosystems (with (ΔTs − Ta) centered at 0 °C, a mean absolute error (MAE) = 1.73 °C and a standard deviation = 2.4 °C). Secondly, comparisons between day MODIS Ts data with maximum Ta showed that (ΔTs − Ta) strongly varies according to the seasonality, the ecosystems, the solar radiation, and cloud-cover. Two factors proposed in the literature to retrieve maximum Ta from Ts, i.e. the Normalized Difference Vegetation Index (NDVI) and the Solar Zenith Angle (SZA), were analyzed. No strong relationship between (ΔTs − Ta) and (i) NDVI and (ii) SZA was observed, therefore requiring further research on robust methods to retrieve maximum Ta.     

Variability in surface BRDF at different spatial scales (30 m-500 m) over a mixed agricultural landscape as retrieved from airborne and satellite spectral measurements

Over the past decade, the role of multiangle remote sensing has been central to the development of algorithms for the retrieval of global land surface properties including models of the bidirectional reflectance distribution function (BRDF), albedo, land cover/dynamics, burned area extent, as well as other key surface biophysical quantities impacted by the anisotropic reflectance characteristics of vegetation. In this study, a new retrieval strategy for fine-to-moderate resolution multiangle observations was developed, based on the operational sequence used to retrieve the Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 5 reflectance and BRDF/albedo products. The algorithm makes use of a semiempirical kernel-driven bidirectional reflectance model to provide estimates of intrinsic albedo (i.e., directional-hemispherical reflectance and bihemispherical reflectance), model parameters describing the BRDF, and extensive quality assurance information. The new retrieval strategy was applied to NASA's Cloud Absorption Radiometer (CAR) data acquired during the 2007 Cloud and Land Surface Interaction Campaign (CLASIC) over the well-instrumented Atmospheric Radiation Measurement Program (ARM) Southern Great Plains (SGP) Cloud and Radiation Testbed (CART) site in Oklahoma, USA. For the case analyzed, we obtained ~ 1.6 million individual surface bidirectional reflectance factor (BRF) retrievals, from nadir to 75° off-nadir, and at spatial resolutions ranging from 3 m to 500 m. This unique dataset was used to examine the interaction of the spatial and angular characteristics of a mixed agricultural landscape; and provided the basis for detailed assessments of: (1) the use of a land cover type-specific a priori knowledge in kernel-driven BRDF model inversions; (2) the interaction between surface reflectance anisotropy and instrument spatial resolution; and (3) the uncertainties that arise when sub-pixel differences in the BRDF are aggregated to a moderate resolution satellite pixel. Results offer empirical evidence concerning the influence of scale and spatial heterogeneity in kernel-driven BRDF models; providing potential new insights into the behavior and characteristics of different surface radiative properties related to land/use cover change and vegetation structure.     

Effects of spectral response function on surface reflectance and NDVI measured with moderate resolution satellite sensors: Extension to AVHRR NOAA-17, 18 and METOP-A

This work extends the previous study of Trishchenko et al. [Trishchenko, A. P., Cihlar, J., & Li, Z. (2002). Effects of spectral response function on surface reflectance and NDVI measured with moderate resolution satellite sensors. Remote Sensing of Environment 81 (1), 1-18] that analyzed the spectral response function (SRF) effect for the Advanced Very High Resolution Radiometer (AVHRR) onboard the NOAA satellites NOAA-6 to NOAA-16 as well as the Moderate Resolution Imaging Spectroradiometer (MODIS), the VEGETATION sensor (VGT) and the Global Imager (GLI). The developed approach is now applied to cover three new AVHRR sensors launched in recent years on NOAA-17, 18, and METOP-A platforms. As in the previous study, the results are provided relative to the reference sensor AVHRR NOAA-9. The differences in reflectance among these three radiometers relative to the AVHRR NOAA-9 are similar to each other and range from − 0.015 to 0.015 (− 20% to + 2% relative) for visible (red) channel, and from − 0.03 to 0.02 (− 5% to 5%) for the near infrared (NIR) channel. The absolute change in the Normalized Difference Vegetation Index (NDVI) ranged from − 0.03 to + 0.06. Due to systematic biases of the visible channels toward smaller values and the NIR channels toward slightly larger values, the overall systematic biases for NDVI are positive. The polynomial approximations are provided for the bulk spectral correction with respect to the AVHRR NOAA-9 for consistency with previous study. Analysis was also conducted for the SRF effect only among the AVHRR-3 type of radiometer on NOAA-15, 16, 17, 18 and METOP-A using AVHRR NOAA-18 as a reference. The results show more consistency between sensors with typical correction being under 5% (or 0.01 in absolute values). The AVHRR METOP-A reveals the most different behavior among the AVHRR-3 group with generally positive bias for visible channel (up to + 5%, relative), slightly negative bias for the NIR channel (1%-2% relative), and negative NDVI bias (− 0.02 to + 0.005). Polynomial corrections are also suggested for normalization of AVHRR on NOAA-15, 16, 17 and METOP-A to AVHRR NOAA-18.     

Land surface emissivity retrieval from combined mid-infrared and thermal infrared data of MSG-SEVIRI

This work addressed the retrieval of Land Surface Emissivity (LSE) from combined mid-infrared and thermal infrared data of Spinning Enhanced Visible and Infra-Red Imager (SEVIRI) onboard the geostationary satellite—Meteosat Second Generation (MSG). To correct for the atmospheric effects in satellite measurements, a new atmospheric correction scheme was developed for both Middle Infra-Red (MIR) and Thermal Infra-Red (TIR) channels. For the MIR channel, because it is less sensitive to the change of water vapor content, the clear-sky and time-nearest European Centre for Median-range Weather Forecast (ECMWF) atmospheric data were used for the images where no atmospheric data are available. For TIR channels, a modified model of Diurnal Temperature Cycle (DTC) used by Göttsche and Olesen [Göttsche, F. M., and Olesen, F. S. (2001). Modeling of diurnal cycles of brightness temperature extracted from METEOSAT data. Remote Sensing of Environment, 76, 337-348.] and Schädlich et al. [Schädlich, S., Göttsche, F. M., and Olesen, F. S. (2001). Influence of land surface parameters and atmosphere on METEOSAT brightness Temperatures and generation of land surface temperature maps by temporally and spatially interpolating atmospheric correction. Remote Sensing of Environment, 75, 39-46.] was adopted. The separation of Land Surface Temperature (LST) and LSE is based on the concept of the Temperature Independent Spectral Indices (TISI) [Becker, F., and Li, Z. L. (1990a). Temperature independent spectral indices in thermal infrared bands. Remote Sensing of Environment, 32, 17-33.] constructed with one channel in MIR and one channel in TIR. The results of two different combinations (combination of channels 4 and 9 and of channels 4 and 10) and two successive days at six specific locations over North Africa show that the retrievals are consistent. The range of emissivity in MSG-SEVIRI channel 4 goes from 0.5 for bare areas to 0.96 for densely vegetated areas, whereas the emissivities in MSG-SEVIRI channels 9 and 10 are usually from 0.9 to 0.95 for bare areas and from 0.95 to 1.0 for vegetated areas. For densely vegetated areas, the emissivities in MSG-SEVIRI channel 9 are larger than the ones in channel 10, whereas the opposite is observed over bare areas. The rms differences between two combinations over the whole studied region are 0.017 for emissivity in channel 4, 0.008 for emissivity in channel 9 and 0.007 for emissivity in channel 10.     

An assessment of photosynthetic light use efficiency from space: Modeling the atmospheric and directional impacts on PRI reflectance

Estimation of photosynthetic light use efficiency (ε) from satellite observations is an important component of climate change research. The photochemical reflectance index, a narrow waveband index based on the reflectance at 531 and 570 nm, allows sampling of the photosynthetic activity of leaves; upscaling of these measurements to landscape and global scales, however, remains challenging. Only a few studies have used spaceborne observations of PRI so far, and research has largely focused on the MODIS sensor. Its daily global coverage and the capacity to detect a narrow reflectance band at 531 nm make it the best available choice for sensing ε from space. Previous results however, have identified a number of key issues with MODIS-based observations of PRI. First, the differences between the footprint of eddy covariance (EC) measurements and the MODIS footprint, which is determined by the sensor's observation geometry make a direct comparison between both data sources challenging and second, the PRI reflectance bands are affected by atmospheric scattering effects confounding the existing physiological signal. In this study we introduce a new approach for upscaling EC based ε measurements to MODIS. First, EC-measured ε values were “translated” into a tower-level optical PRI signal using AMSPEC, an automated multi-angular, tower-based spectroradiometer instrument. AMSPEC enabled us to adjust tower-measured PRI values to the individual viewing geometry of each MODIS overpass. Second, MODIS data were atmospherically corrected using a Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm, which uses a time series approach and an image-based rather than pixel-based processing for simultaneous retrievals of atmospheric aerosol and surface bidirectional reflectance (BRDF). Using this approach, we found a strong relationship between tower-based and spaceborne reflectance measurements (r² = 0.74, p < 0.01) throughout the vegetation period of 2006. Swath (non-gridded) observations yielded stronger correlations than gridded data (r² = 0.58, p < 0.01) both of which included forward and backscatter observations. Spaceborne PRI values were strongly related to canopy shadow fractions and varied with different levels of ε. We conclude that MAIAC-corrected MODIS observations were able to track the site-level physiological changes from space throughout the observation period.     

A scalable approach to mapping annual land cover at 250 m using MODIS time series data: A case study in the Dry Chaco ecoregion of South America

Land use and land cover (LULC) maps from remote sensing are vital for monitoring, understanding and predicting the effects of complex human-nature interactions that span local, regional and global scales. We present a method to map annual LULC at a regional spatial scale with source data and processing techniques that permit scaling to broader spatial and temporal scales, while maintaining a consistent classification scheme and accuracy. Using the Dry Chaco ecoregion in Argentina, Bolivia and Paraguay as a test site, we derived a suite of predictor variables from 2001 to 2007 from the MODIS 250 m vegetation index product (MOD13Q1). These variables included: annual statistics of red, near infrared, and enhanced vegetation index (EVI), phenological metrics derived from EVI time series data, and slope and elevation. For reference data, we visually interpreted percent cover of eight classes at locations with high-resolution QuickBird imagery in Google Earth. An adjustable majority cover threshold was used to assign samples to a dominant class. When compared to field data, we found this imagery to have georeferencing error < 5% the length of a MODIS pixel, while most class interpretation error was related to confusion between agriculture and herbaceous vegetation. We used the Random Forests classifier to identify the best sets of predictor variables and percent cover thresholds for discriminating our LULC classes. The best variable set included all predictor variables and a cover threshold of 80%. This optimal Random Forests was used to map LULC for each year between 2001 and 2007, followed by a per-pixel, 3-year temporal filter to remove disallowed LULC transitions. Our sequence of maps had an overall accuracy of 79.3%, producer accuracy from 51.4% (plantation) to 95.8% (woody vegetation), and user accuracy from 58.9% (herbaceous vegetation) to 100.0% (water). We attributed map class confusion to limited spectral information, sub-pixel spectral mixing, georeferencing error and human error in interpreting reference samples. We used our maps to assess woody vegetation change in the Dry Chaco from 2002 to 2006, which was characterized by rapid deforestation related to soybean and planted pasture expansion. This method can be easily applied to other regions or continents to produce spatially and temporally consistent information on annual LULC.