Land cover mapping applications with MODIS: a literature review

Abstract

Land use/land cover monitoring and mapping is crucial to efficient management of the land and its resources. Since the late 1980s increased attention has been paid to the use of coarse resolution optical data. The Moderate Resolution Imaging Spectroradiometer (MODIS) has features, which make it particularly suitable to earth characterization purposes. MODIS has 10 products dedicated mainly to land cover characterization and provides three kinds of data: angular, spectral and temporal. MODIS data also includes information about the data quality through the ‘Quality Assessment’ product. In this paper, we review how MODIS data are used to map land cover including the preferred MODIS products, the preprocessing and classification approaches, the accuracy assessment, and the results obtained.

Keywords:

• MODIS
• mapping
• image processing
• classification, accuracy
• land use/land cover

1. Introduction

According to projections by the United Nations, the world population is expected to surpass 9 billion people by 2050. Most of the additional 2.3 billion people will add to the population of developing countries, which are projected to rise from 5.6 billion in 2009 to 7.9 billion in 2050. Available studies indicate that approximately 16% of the total arable land area has been degraded by human activity during the past half century (FAO Progress Report, 1997). To manage the daunting social challenges ahead, land use/land cover monitoring and mapping is crucial to efficient management of the land and its resources.

Reliable observations of the terrestrial environment are of crucial importance to understand climate change and its impacts, to sustain economic development, to properly manage natural resources, to promote conservation, to preserve biodiversity and to improve the scientific understanding of ecosystems (Herold et al. 2006). Satellite remote sensing has long been considered an ideal technology for land use/land cover monitoring and mapping due to its ability to provide synoptic and repetitive observations of the vegetation cover (Franklin and Wulder 2002). DeFries (1995) noted that the functional controls on exchanges of water, energy, and carbon dioxide between the atmosphere and biosphere are well enough understood that it is feasible to model these exchanges using a small number of vegetation characteristics. Many of these characteristics can be obtained in a cost effective manner through remote sensing coupled with some ground-based information.

Since the late 1980s, increased attention has been paid to the use of coarse resolution optical data, represented primarily by National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) images. The AVHRR were initially available at 8 km resolution and later, at the nominal resolution of 1 km for all land areas of the globe. The cost of data has gradually decreased (especially data for research purposes), and some data are available through direct broadcast at no charge. The launch of new satellite sensors such as SPOT 4 VEGETATION (VGT) (Saint 1992), Moderate Resolution Imaging Spectroradiometer (MODIS) (Salomonson et al. 1989, Running et al. 1994, Barnes et al. 1998, Wooster 2007), Medium Resolution Imaging Spectrometer (MERIS) with a systematic global acquisition strategy is inaugurating a new era in land remote sensing during which (1) high quality data-sets are available globally for land cover mapping applications; and (2) the remote sensing research community is developing sound methodologies for generating land cover information products (Cihlar 2000).

Among these satellite sensors, MODIS has features that make it particularly suitable to earth characterization. Early studies (Hansen 2002) indicated that MODIS data would provide a substantial improvement over AVHRR in mapping tree cover. The spatial detail available in MODIS imagery is unprecedented for moderate resolution satellites; however, preserving the finest spatial detail within the compositing process requires new approaches.

MODIS has two nearly identical sensors on-board two separate satellites (Terra and Aqua), MODIS surveys the entire Earth's surface every 1–2 days, acquiring data in 36 spectral bands between 0.405 and 14.385 µm, with the first seven bands focused on land observation (the so-called ‘land’ bands provide two spatial resolutions of 250 and 500 m in visible and infrared regions) and has a viewing swath width of 2330 km. Many of the data products derived from MODIS observations describe features of the land, the oceans and the atmosphere that can be used for studies of processes and trends on regional to global scales. Users with an appropriate X-band receiving system may capture regional data directly from the satellite using the MODIS direct broadcast signal.

This paper is a review of the prior studies using MODIS data in land use/land cover mapping. It is organized as follows: the first section outlines the MODIS land oriented products and their main characteristics and the way they are used in land cover mapping efforts; the second section presents preprocessing procedures; and the third section discusses the classification approaches. The fourth section is an examination of accuracy assessment including the data used as reference. The fifth section outlines the results obtained including the classification scheme, the number of categories and the accuracy achieved. Finally, the paper closes with a discussion of some of the drawbacks encountered by users when utilizing MODIS and also the positive outcomes resulting from MODIS development and its use in the remote sensing community.

2. MODIS data products for land use/land cover mapping

As noted above, many of the data products derived from MODIS observations describe features of the land, the oceans and the atmosphere that can be used for studies of processes and trends on local to global scales. The MODIS land products (MODLAND products) are generated using data from Terra, Aqua or a combination of both in a hierarchy of processing levels: retrieved geophysical parameters at the same location as the MODIS instrument data (Level 2), earth-gridded geophysical parameters (Level 2G and Level 3), and earth-gridded model outputs (Level 4). The smallest unit of MODIS Land data processed at any one time is defined at Level 2 as a granule, and at Levels 2G, 3, and 4 as a tile. A granule corresponds to 5 minutes of MODIS sensing such that there are approximately 2340×2330 km. Tiles are 10 degrees by 10 degrees in a Sinusoidal Grid (SIN) at the equator. Additionally, the Climate Modelling Grid (CMG) images are gridded on a simple 0.05 degree (5600 m) geographic projection. The MODIS land products are produced at four nominal spatial resolutions 250 m, 500 m, 1000 m, and 0.05 degrees (NASA MODIS Web-Components of MODIS 2008).

According to Morisette et al. (2002) the products focused on land characterization can be divided into three main groups including:

•	Radiation budget variables: (1) surface reflectance (MOD09); (2) land surface temperature (LST) and emissivity (MOD11); (3) surface reflectance, bidirectional reflectance distribution function (BRDF)/Albedo parameter (MOD43); and (4) snow cover (MOD10).
•	Ecosystem variables: (1) gridded vegetation indices (MOD13), (2) leaf area index (LAI), and fractional photosynthetically active radiation (FPAR) (MOD15), and (3) vegetation production, net primary productivity (NPP) (MOD17).
•	Land cover characteristics: (1) land cover/land cover change (MOD12), (2) thermal anomalies, fires and biomass burning (MOD14), (3) vegetation cover conversion (MOD44), and (4) burn area (MOD45).

MOD09 is an estimate of the surface spectral reflectance for the seven MODIS ‘Land’ bands in visible and near infrared as it would have been measured at ground level if there were no atmospheric scattering or absorption. MOD43 provides the coefficients for mathematical functions that describe the BRDF of each pixel in these seven bands, nadir BRDF-adjusted reflectance (NBAR) and albedo measures. MOD11 contains LST during the day and the night and the emissivity in several spectral bands. MOD10 products provide fractional snow cover, and snow albedo.

MOD13 includes the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI). MOD15 presents the FPAR, which is the fraction of the incoming solar radiation that a plant canopy absorbs for photosynthesis and growth in the 0.4–0.7 nm spectral range. MOD15 also includes the LAI, which is the biomass equivalent of FPAR and is defined as the ratio of total upper leaf surface of vegetation divided by the surface area of the land on which the vegetation grows. Both are biophysical variables which describe canopy structure and are related to functional process rates of energy and mass exchange. MOD17 provides both NPP and Gross Primary Production (GPP). GPP is the rate at which vegetation captures and stores biomass during photosynthesis. Some fraction of this fixed energy is used for cellular respiration. The remaining fixed energy is referred to as NPP. GPP/NPP provides information on several vegetation characteristics that have important social and economic impacts, including crop yields, range land forage, and forest production.

MOD12 consists of two suites of data-sets MOD12Q1 (global land cover) and MOD12Q2 (global land cover dynamics). MOD12Q1 includes five layers depicting different global land cover classifications while MOD12Q2 characterizes seasonal time scale dynamics (phenology). MOD14 includes fire occurrence (day/night), fire location, and an energy calculation for each fire. The MOD44 suite provides a database of land cover and change: MOD44A vegetative cover conversion provides information on land cover change. The MOD44B vegetation continuous fields product provides sub-pixel estimates of percent of woody, herbaceous, and bare cover. MOD45 provides information on burned areas (date and extension).

Each product may include various sub-products with different spatial or temporal resolutions and features. Table 1 shows the number of sub-products, processing level, spatial and temporal resolutions, and sensor for each product. It is worth noting that MODIS products are undergoing constant improvement and often present differences between versions. Most of these products have the potential to be used as input data to generate land use/land cover cartography, while MOD12, MOD44, and MOD45 are finished products that can be displayed as thematic maps.

Table 1. MODLAND products main characteristic.

Acronym	Name	Sub-products	Level	Spatial representation	Spatial resolution	Temporal resolution	Sensor^a
MOD09	Surface reflectance	5	2G, 3	Tile, CMG	250, 500, 1000 m0.05°	Day, 8 days	T, A
MOD11	Surface temperature and emissivity	7	2, 3	Granule, tile, CMG	1000, 6000 m0.05°	Day, 8 days, month	T, A
MOD12	Land cover/land cover change	4	3	Tile, CMG	500, 1000 m0.05°	Year	T, C
MOD13	Vegetation indices	6	3	Tile, CMG	250, 500, 1000 m0.05°	16 days, month	T, A
MOD14	Thermal anomalies, fires and biomass burning	3	2, 3	Granule, tile	1000 m	Day, 8 days	T, A
MOD15	LAI and FPAR	1	4	Tile	1000 m	8 days	T, A, C
MOD17	Vegetation production, NPP	1	4	Tile	1000 m	8 days	T, A
MOD43	Surface reflectance/albedo	12	3	Tile, CMG	500, 1000 m0.05°	16 days	C
MOD44	Vegetation cover conversion	2	3	Tile	250, 500 m	1 year	T
MOD45	Burn area	1	3	Tile	500 m	1 month	C
^aT, Terra; A, Aqua; C, combined.

The MODIS input data provides three different dimensions of sensor information:

1.	Spectral. The spectral information in a remotely sensed image is an estimate of the surface reflectance, as it would have been measured at ground level if there were not atmospheric scattering or absorption. This information includes daily images and multi-date composites. The most direct product is MOD09 daily images, but there are also MOD09 8-day composites. The vegetation indices are the most common transformations of the spectral information with the purpose of enhancing the green vegetation.
2.	Temporal. The high temporal coverage of the MODIS land products is advantageous for many studies since it allows for analysis over time. For example, analysing the vegetation phenology in a time series allows us to more readily distinguish classes with similar spectral response rather than through the evaluation of a single date.
3.	Angular. Due to the wide field of view provided by the MODIS sensors, it is possible to have records of the same ground pixel sampled under different geometries: object–sensor–sun. Because the reflectance distribution of many land covers (e.g. vegetation) is strongly anisotropic, multiview angle observations contains additional information beyond that acquired with nadir or single angle spectral measurements alone (Asner et al. 1998).

Additional information included in MODIS data are the data quality through the ‘Quality assessment’ (QA) product, made available on a routine basis and stored as product metadata and as per pixel information (Roy et al. 2002). At the per pixel level, the QA product is coded in 16-bit format where specific bit ranges provide a specific quality parameter. The most common data quality parameters indicate atmospheric conditions, for example the presence of clouds or give a general quality index (usefulness index). Other indices take into account parameters such as solar, viewing geometry, etc.

2.1 Studies based on vegetation indices

Despite the fact that there are a high number of available products that can be used in land cover/land use mapping, most of the reviewed authors use NDVI time series. NDVI is used more frequently than other indices largely to continue the AVHRR heritage and to facilitate the comparison with other studies (Gitelson and Kaufman 1998).

Vegetation types can be characterized using their seasonal or phenological variations. Annual seasonal parameters of reflectance data and NDVI, and metrics such as minimum, maximum, and amplitude, have the potential to improve land cover separability (Townshend et al. 1991, DeFries et al. 1995, Ji and Peters 2007). The effectiveness of NDVI time series has been demonstrated in several studies on crop discrimination of North America, Europe, and Asia (Mingwei and Qingbo 2008).

Some authors extract the seasonal parameters by using more than one year of NDVI records to retain sufficient data after eliminating noise (Tottrup 2007). Other authors use partial time series, restricting the data to a particular time of interest; for example, the growing season (Bagan et al. 2005, Sivanpillai and Latchininsky 2007), or a snow-free period (Heiskanen and Kivinen 2008).

Even though the use of NDVI is preferred by most authors, the use of other indices has been increasing. One of the most popular indices is the EVI provided in the MOD13 product. EVI was developed to optimize the vegetation signal with improved sensitivity in high biomass regions and improved vegetation monitoring through a decoupling of the canopy background signal and a reduction in atmospheric influences (Huete et al. 2002). Several authors indicate that EVI is more sensitive to variations in green vegetation cover when compared with NDVI (Justice et al. 1998, Ferreira et al. 2003, Xavier et al. 2006). NDVI exhibits scaling problems, asymptotic (saturated) signals over high biomass conditions, particularly strong with high canopy background brightness (Huete 1988). Xavier et al. (2006) observed that EVI is sensitive to variations in land use cover due to phenology and land management practices.

Other less commonly used vegetation indices are more suited to particular monitoring requirements. Gitelson et al. (2005) suggested the use of the Green Index (GI) defined as: GI = R_nir/R_green (where R_nir is the near infrared band and R_green is the green band). In the green spectrum, absorption of light is high enough to provide high sensitivity of the GI to chlorophyll content but much lower than in the blue and red to avoid saturation (Gitelson et al. 2003). Ozdogan and Gutman (2008) proved that the (GI) has a better performance than NDVI and EVI in the subpixel monitoring of irrigated crops. This can be explained by the fact that irrigated crops with very little or no soil moisture stress exhibit larger chlorophyll content than non-irrigated crops with potential moisture stress.

2.2 Studies based on surface reflectance data

Carrão et al. (2008) used multi-date composite surface reflectance products. Although preliminary investigations showed that the variable view geometry of the composites caused adjacent pixels to have large radiometric differences for the same cover type, leading to classification errors (Cihlar 2000). To avoid these difficulties, Wessels et al. (2004), Chen and Rao (2009) and Shimabukuro et al. (2009) used single date near cloud free imagery to map regional land cover. Sedano et al. (2005) used five single date images, as only five cloud-free single-day images corresponding to different periods of the dry season were available for southern Africa in this study.

2.3 Utilizing more than one kind of MODIS data

Some researchers find it useful to incorporate different kinds of MODIS data. For example, in order to generate the land cover and land cover dynamics products (MOD12), Friedl et al. (2010) added to bands 1–7 the information provided by the EVI and the LST. To compute the EVI they used NBAR.

Using more than one kind of data often enhances the spectral separability of the desired classes. Westra and De Wulf (2007) used information from medium infrared, near infrared, and an index based on these two bands, called the Normalized Difference Water Index combined with NDVI and EVI (Gao 1996). This method demonstrated a better performance than using only one variable by itself.

Another approach is the combination of more than one kind of sensor data. Braswell et al. (2003) took advantage of the nearly coincident imaging of the Multi-angle Imaging Spectroradiometer (MISR) and MODIS (using multiple shortwave infrared bands), and merged these sets of data to map sub-pixel land cover fractions in the Brazilian Amazon region. They obtained an apparent synergistic effect providing an increase of 20 and 30%, respectively, for the correlation value and the average accuracy of the derived maps.

Yang (2006) employed a data mining technique to evaluate the degree to which the accuracy of land cover classification can be increased using multi-spectral, multi-temporal, and multi-angle remote sensing data in a study area in South central New Mexico, USA. Data used for this study included EOS MISR surface BRF and MODIS 16-day NDVI composites acquired from 2002 to 2003. Eight land cover types were classified using a decision tree algorithm with multiple classifications obtained to evaluate classification accuracy using different input data (MODIS data only, MISR data only, and MODIS plus MISR data). The largest increase in accuracy was observed for open oak woodland, coniferous woodland, and woody wetland using MISR and MODIS data, and for irrigated cropland and barren land using MISR data only.

Acerbi-Junior et al. (2006) performed the fusion between MODIS and Landsat images using the wavelet transform. Their results indicated that the fusion of MODIS and TM images using the Fourier space wavelet transform provided the best quality measures for the fused images. This methodology can be beneficial in areas where there are gaps in the data series. Classification results showed that the fused images could be used for their study area with an accuracy level comparable to the Landsat TM image.

Hansen et al. (2008) used MODIS derived forest cover products to calibrate Landsat data allowing a standard processing stream that increases the internal consistency of the regional-scale cover characterizations. They found a standard error of the global MODIS vegetation continuous fields product for sites tested of 11.5%, according to the authors this accuracy is sufficient for identifying the broad tree cover strata required to map forest cover and change within the Congo Basin. They concluded that the multi-resolution methodology is portable, but requires study areas where the moderate and high spatial resolution thematic classes are the same and few in number, and where the land cover classes are spatially homogenous at scales larger than several moderate resolution pixels.

2.4 Three dimensions of data

A few authors compare the respective contributions of the three dimensions of data (temporal, spectral, and angular). Carrão et al. (2008) found that spectral information is more useful in discriminating land cover classes than temporal information. They demonstrated that the multi-temporal factor (choosing multiple dates) has a significant effect when coupled with combinations of few spectral bands, but the temporal contribution disappeared as soon as the full spectral information (seven bands) is exploited. In contrast, even with a full year of MODIS data, the results are strengthened when continuing to use no less than three spectral bands.

Heiskanen and Kivinen (2008) demonstrated that multi-temporal and angular variables can increase the accuracy of the cover estimates and forest mapping. Their results also suggest that seasonality affects the model performance; the late-spring and early-summer data providing superior results compared to mid- and late-summer data.

2.5 Incorporation of ancillary data

The spectral similarity of many land surface types presents a challenge in differentiating these classes by using only spectral information. This problem grows as the sensor resolution decreases. Many authors have tried to use other sources of information in order to improve the final accuracy of the maps. For instance, to separate vegetation types with a similar spectral response but located in different ecological zones, it is useful to take into account the ecological conditions. The most frequently used ancillary data are elevation because it is highly related to ecological factors and it is easily obtained. Slope can be readily derived from the elevation model, and temperature is another element used in analyses. Some authors chose to use a combination of these ancillary data instead of only one.

Miettinen et al. (2008) used elevation data as ancillary data to separate 12 classes in South-east Asia. Friedl et al. (2010) used a combination of MODIS products and maps of prior probabilities to create the MOD12 collection five including reflectance, LST, BRDF, and EVI. García-Mora and Mas (2010) used MODIS data along with an elevation model and a map of soils to map land cover in a highly diverse region of Mexico. The incorporation of these ancillary data provided only a small increase in accuracy which was attributed to the overlap of the distribution of various land cover types with respect to these variables. Zhang et al. (2008a, 2008b) used slope together with EVI to map land cover in Northern China. The overall accuracy obtained significantly exceeded the accuracy of the MOD12 product previously produced by NASA and USGS. Ozdogan and Gutman (2008) used an effective irrigation potential map derived from climate and radiation data to map irrigated areas from multi-date MODIS data. Table 2 summarizes the different mapping efforts reviewed in this study.

Table 2. Mapping efforts with MODIS data.

Author(s) (date)	Kind of mapping	MODIS products	Ancillary data	Study area location
Mingwai and Qingbo (2008)	Crop discrimination	NDVI time series	No	Northern china
Tottrup (2007)	Fractional forest cover	NDVI time series	Elevation data	Southeast Asia
Bagan et al. (2005)	Land cover	EVI growing season time series	No	Arid and semi-arid area of Northern China
Heiskanen and Kivinen (2008)	Forestry	MOD12	No	Finlandia
Xavier et al. (2006)	Sugarcane crops	EVI time series	No	Sao Paulo, Brazil
Ozdogan and Gutman (2008)	Mapping irrigated areas	Green Index	Climate-based indices of surface moisture status map	USA
Carrão et al. (2008)	Land cover	Surface Reflectance composites	No	Portuguese mainland
Sedano et al. (2005)	Land cover	5 single-day Surface reflectance bands 1–7	Land cover maps, Elevation	Miombo woodlands, Zambezia, Mozambique.
Friedl et al. (2010)	Land cover and L/c dynamics	EVI, NDVI, bands 1–7	LST BRDF NBAR	Global
Westra and De Wulf (2007)	Land cover/land use	NDVI, EVI	Normalized Water Index	Semi-arid zone of Africa
Zhang et al. (2008a)	Land cover	Bands 1–7, NDVI snow cover, LST BRDF	Terrain elevation Texture/Water mask	North China plain
Braswell et al. (2003)	Land cover	VIS NIR SWIR	MISR	Brazilian Amazon region
Matsuoka et al. (2007)	Land cover	MOD09Q1	MOD10A2	East Asia
Miettinen et al. (2008)	Land cover		Elevation information and peatland maps	Malaysian peninsula and the major islands
Sivanpillai and Latchininsky (2007)	Mapping locust habitats	Multi-temporal MODIS imagery 250M		Amudarya River Delta, Uzbekistan
Wessels et al. (2004)	Mapping regional land cover	Single date near-cloud free 250 and 500 m		Greater Yellowstone
Shimabukuro et al. (2009)	Burnt area mapping	Daily reflectance (MOD09) for 5 days	No	Western Brazilian Amazonia

3. Preprocessing

Due to its high temporal resolution and large number of bands and products, the MODIS data sets used to map land cover tend to present a high dimensionality and are therefore prone to suffer the Hughes Effect which is a decrease in classification accuracy when the number of features exceeds a given limit, for a fixed training-sample size. A reduction in the number of features overcomes this problem, thus improving classification accuracy (Hughes 1968, Landgrebe 2003). To minimize the effects of correlation and high dimensionality, some form of dimension reduction is a desirable preprocessing step. Feature selection (selecting the most relevant attributes) and feature extraction (combining attributes into a new reduced set of features) are two approaches to dimension reduction.

A second objective of preprocessing is to reduce noise. The process of creating temporal 8- or 16-day composites can be seen as a noise reduction process where pixels with cloud contamination or highly oblique view angles are eliminated. However, even after noise-reduction, data are often still strongly contaminated particularly in cloudy areas or during specific periods of the year. Additional procedures such as time series filtering may then be applied in order to further reduce noise. In fact, many preprocessing methods cannot be strictly divided into either dimension reduction or noise reduction because many methods accomplish both purposes. In addition, methods presented in the time-series analysis section (Section 3.3) such as Fourier and wavelet transforms can also be considered as feature extraction methods.

3.1 Feature selection

A simple way to select input data and reduce dimensionality is taking advantage of the MODIS quality assessment (QA) data to carry out data quality filtering. Since quality information is provided for both the entire product and for each individual pixel, the selection can serve to discard the entire data-set (e.g. the entire tile) or specific pixels or groups of pixels inside the image. For example, Sedano et al. (2005) used the QA per pixel information to select images without clouds and smoke reducing the number of bands from an entire year to five cloud/smoke-free single day images. García-Mora and Mas (2010) extracted the usefulness index of the vegetation index and the MODLAND QA index for reflectance data to discard entire images by taking into account two criteria: (1) the quality value average of the image and (2) the proportion of pixels reported with bad quality. There were set threshold values, which were customized for each study area in order to provide flexibility. For example, in tropical areas which are cloudy during significant parts of the year, the selection criteria for each image had to be adjusted in order to retain sufficient images for evaluation. Within the reduced data-set, a visual examination of the selected images was still necessary since bad quality images remained. A certain amount of poor quality images were not labelled appropriately in the QA index. Sivanpillai and Latchininsky (2007) also used quality control data provided by the LPDAAC in combination with visual analyses to screen for clouds and other anomalies such as missing data and image striping.

Some authors chose to build masks in order to remove specific pixels under a pre-established threshold and thus, do not base quality filtering on QA information alone. For instance, Xiao (2005) observed that thick clouds still remained in pixels which were not labelled as clouds by the cloud quality flag. Xiao applied an additional filter to the pixels by setting a blue reflectance threshold of 0.2. Though this results in an image with data gaps; these gaps can later be filled by methods such as time filtering.

When the input data are highly redundant, feature selection can be accomplished using separability measures. These measures allow the evaluation of the statistical distance between samples in different land cover categories in a space defined by a given set of bands. Computing the statistical distances between each of the pairs of categories for each of the possible subsets of bands allow the determination of the optimal number of input bands (the minimum number of bands on a subset allowing the maximum or a near maximum separability) and the optimal input band combination. Popular measure indices are the Bhattacharyya, the Mahalanobis and the Jeffreys-Matusita distances along with the divergence and the transformed divergence (Landgrebe 2003).

Carrão et al. (2008) used the median of the Mahalanobis interclass distances as a sorting criterion to measure the temporal and the spectral features in their land cover classification study. García-Mora and Mas (2011) used the Transformed Divergence (TD) since according to Maussel et al. (1990), it has been shown to correlate well with the results of the maximum likelihood (ML) classifier. However, García-Mora and Mas observed that separability measurement analyses tend to select noisy images since noise was interpreted as new uncorrelated information which increased separability. These authors carried out a two-step selection procedure first selecting noise free products using the QA information together with visual inspection before computing separability to select an optimal subset.

3.2 Feature extraction

When the input data are too large to be processed efficiently and are suspected to be highly redundant, then feature extraction should be applied to the input data to produce a reduced representation of the feature set. If the extracted features are carefully chosen it is expected that the feature set will extract the relevant information from the input data in order to perform the desired task using this reduced representation instead of the full size input.

The most common feature extraction method in remote sensing is the principal component analysis (PCA) which is an image transformation method that allows the concentration of useful information into the first principal components, which are uncorrelated. Ferreira et al. (2004) applied a PCA to MODIS-like data. When the first three components, which accounted for 99.5% of the total variance, are considered simultaneously, approximately 91% of the total data-set was correctly classified. Chen (2005) applied three PCA to reflectance data with 250, 500, and 1000 m spatial resolution, respectively, and found that more than 95% of variance of the two 250 m bands was in the first component and that 500 and 1000 m reflectance can be reduced to two or three components with a loss under 10% of the total variance at the different sites he evaluated. Sedano et al. (2005) used PCA in a different way; they applied PCA as a change detection technique to enhance the seasonal changes of the different vegetation types in the study area.

Another image transformation method is the Tasseled Cap Transformation (TCT). This method, initially developed for Landsat MSS images, was inspired by the PCA although TCT shows important advantages over PCA since it provides concrete physical and easily interpretable indices (brightness, greenness, and wetness). The coefficients are set in an empirical way and are valid for only one sensor. TCT was adapted to MODIS by Lobser and Cohen (2007) who showed that MODIS TCT has great potential for the analysis of land cover characteristics. Shimabukuro et al. (2009) produced shade fraction images using a linear spectral mixture model to map burnt areas.

3.3 Time series filtering and analysis

Due to the sensor's temporal resolution, an interesting feature of MODIS is its ability to produce time series, especially for vegetation index time series. Time series data has a natural temporal ordering which makes their analysis distinct from other common data-sets, in which there is no natural ordering of the observations (e.g. spectral bands). Therefore, times series analysis tools can be employed, most of which have been previously applied to AVHRR times series data-sets. It should be noted that Fourier and wavelet transforms are procedures based mainly on time filtering.

Time series filtering consists of smoothing aberrant values at a certain date by taking into account the previous and following values in the time series. Smoothing can be done by a simple moving average or through more sophisticated fitting methods and can be applied to all the observations or only those detected as being noisy. The most commonly used smoothing functions applied to time-series data include: the moving-average filter, the asymmetric Gaussian, the double logistic, and the adaptive Savitzky-Golay filter. The moving-average filter calculates a smooth value of an observation averaging its value with previous and posterior values. The size of the time-window (number of previous and posterior observations taken into account for each average) determines the smoothing effect. An alternative is to use the median instead of the average. The adaptive Savitzky-Golay filtering approach uses local polynomial functions in fitting. It is able to capture subtle and rapid changes in the time series but this also makes it sensitive to noise. Both asymmetric Gaussian and double logistic approaches use semi-local methods. These are less sensitive to noise and can give a better description at the beginning and end of the seasons (Jonsson and Eklundh 2002).

Boschetti et al. (2009) used the Savitzky-Golay approach to smooth NDVI data. Zhang et al. (2008a, 2008b) also applied the Savitzky-Golay method to remove pixels labelled as clouds in the QA. Even though this filter compensates the loss of information in the gaps of the masked image to a certain degree, this is not equivalent to having the original information for the entire image. To avoid biases due to time-specific contaminations, Van Dijk et al. (1987) suggested that a time-sliding median filter window outperforms other smoothing techniques, such as linear filters or polynomial best fits, by removing aberrant measures from the time series. This approach was also used by Carrão et al. (2008).

Fourier and wavelet analysis are commonly used for noise removing and curve fitting in MODIS vegetation indices data-sets ( Wessels et al. 2004, Bruce et al. 2006, Colditz et al. 2008). Fourier transform enables the decomposition of temporal data to the frequency domain, whereby the frequency information is represented as constituent sine and cosine waves. These two waves can then be combined into a single cosine wave, with characteristic amplitude (size of the wave) and phase angle (offset of the wave) (Morettin and Toloi 2006, Westra and De Wulf 2007). If the original data are discrete rather than continuous, the discrete Fourier transform (DFT), which requires regular spacing of samples within the temporal domains, should be applied.

In wavelet transform, the function is reconstructed as a linear combination of wavelet and scaling function and their corresponding coefficients. Wavelets cut up data into different frequency components, where each component can be studied with a resolution matched to its scale. Wavelets have advantages over traditional Fourier methods in situations where the signal contains discontinuities and sharp spikes.

Fourier analysis provides a new representation of the time-series of images, which allows the examination of the vegetation phenology using only the amplitude and phase of the most important periodic components. For instance, images of amplitude and phase values were used as attributes to vegetation mapping of Southern Africa. Amplitude and phase angle images were correlated with crop type information and used to discriminate crops (Jakubauskas et al. 2001). Westra and De Wulf (2007) used amplitude, amplitude variance, and phase images to map flooded areas.

Bruce et al. (2006) proposed a novel wavelet-based feature extraction method which outperformed Fourier transform when classifying two types of vegetation from an NDVI times series. These authors found that the wavelet-based features measure more detailed characteristics than the Fourier transform, and the method is more sensitive to noisy spikes in the MODIS temporal signatures also performed better at discriminating the two vegetation classes. Sakamoto et al. (2005) and Oliveira et al. (2009) also found that the wavelet transform performed better than the Fourier transform to filter vegetation indices time series.

There are several very useful programs for analysing and processing time series images. The TIMESAT program (Jonsson and Eklundh 2002, 2004) provides three different smoothing functions to fit the time-series data: asymmetric Gaussian, double logistic, and adaptive Savitzky-Golay filtering. The program allows the user to develop their own algorithms. For instance, Nightingale et al. (2007) used a modified version of the TIMESAT software to smooth and spatially gap-fill FPAR time series data. This approach provides a weighting mechanism based on the MODIS quality assessment (QA) layers by maximizing the use of high-quality retrievals to fit an annual curve. The curve is used to replace missing or poor-quality observations. If large gaps exist in the time-series, TIMESAT does not fit a curve and a separate spatial gap-filling procedure is used. Another program is the time-series generator developed by Colditz et al. (2008), which analyses the pixel-level quality-assurance science data-sets of all gridded MODIS land (MODLand) products suitable for time-series generation. According to user-defined settings, the tool visualizes the spatial and temporal data by generating two indices; the number of invalid pixels and the maximum gap length. Quality settings can be modified spatially and temporally to account for regional and seasonal variations of data quality. The user compares several quality settings and masks or interpolates the data gaps. The program Harmonic Analysis of Time-Series (HANTS) deal with time-series of irregularly spaced observations and identifies cloud contaminated observations (Roerink et al. 2000). It performs two tasks: (1) screening and removal of cloud affected observations; and (2) temporal interpolation of the remaining observations to reconstruct gap-less images at a prescribed time. It calculates a Fourier series to model a time series of pixelwise observations, while at the same time identifying outliers relative to the model of the time series. The algorithm eliminates such outliers and replaces them with the value given by the Fourier series. Atkinson et al. (2009) compared HANTS and TIMESAT programs.

4. Methods of classification

In general terms there are two broad approaches to land cover classification, the hard or crisp approach is the most used and consists of the assignation of one pixel to one class. The fuzzy or soft approach consists of associating various categories to one pixel through membership values, which typically range from 0 (no membership) to 1 (full membership). Fuzzy classifications are useful to express classification ambiguities in transitional areas (ecotone) or pixel mixing. In this last case, the membership values represent the proportions of pixel covered by the different categories, and overcome the problems related to coarse spatial resolution of data such as MODIS. Ozdogan and Gutman (2008) used the fuzzy approach to estimate proportions of irrigated crops, Heiskanen (2008) proportions of tree covers, Braswell et al. (2003) and Tottrup (2007) proportions of mature forest, secondary forest, and non-forest.

Over the last few decades, a large number of classification methods have been developed. Depending on the way classification parameters are handled, all are able to produce hard or fuzzy classifications. A great variety of classification methods have been used to process MODIS data. A few authors used the ML algorithm as did Zhang et al. (2008a, 2008b) to map different types of crops. Many authors used the ML as a benchmark when testing a new classification method (for example, Bagan et al. 2005, Westra and De Wulf 2007, Chen and Rao 2009).

A variety of neural network methods have been applied to process MODIS data such as self-organizing map (SOM) (Bagan et al. 2005), Fuzzy ARTMAP (Borak and Strahler 1999, Sedano et al. 2005), and LNNS, an in-house developed neural network simulator (Westra and De Wulf 2007).Borak and Strahler (1999) and Bagan et al. (2005) reported superior accuracy using neural network than ML.

Decision tree (DT) is also a commonly used algorithm. For instance, Matsuoka et al. (2007) used this approach to map 11 land cover types in East Asia. The DT was defined manually based on the experiment because of insufficient training data and ease of tuning by visual interpretation. Tottrup (2007) used a regression tree algorithm which produces a rule-based model for predicting a single continuous response variable (i.e. fractional land cover) from multi-date MODIS reflectance data. Regression trees are built through an iterative process of splitting the data into subsets. The single explanatory variable that maximizes the reduction in the residual sum of squares is used to split the data, and the process is repeated on the two descendant nodes. Generally, DT allows better results compared to other methodologies. For instance, Chen and Rao (2009) reported that the DT approach outperformed the ML (overall accuracy of 64.5 vs. 48.9%) in classifying MODIS reflectance data to map land cover. Stern et al. (2009) also found that DT obtained higher accuracy than ISODATA in land use classification.

Some authors use less common approaches such as Heiskanen and Kivinen (2008) who used a binomial generalized regression model (GLM) to estimate tree cover, Shimabukuro et al. (2009) who used segmentation and clustering unsupervised classification, Gao et al. (2009) which applied object-based classification, Gonçalves et al. (2006) and Carrão et al. (2008) who used Support Vector Machine.

5. Accuracy assessment

Accuracy assessment is an essential task implemented in the final phase of the land cover classification. However, it is strongly constrained by available resources, and it remains a particular challenge, especially for global or regional products (DeFries et al. 1998, Cihlar 2000). It is important to rely on a robust reference data-set to calculate the accuracy of the final map. The most commonly used approach for accuracy assessment of MODIS classifications is the comparison against a classification obtained from higher spectral resolution imagery (Friedl et al. 2002, Latifovic et al. 2004), other land cover maps (Liu et al. 2003, Boles et al. 2004, Han et al. 2004, Giri et al. 2005) or a subset of the training data (cross validation) (Friedl et al. 2010). The comparisons with census and field work are less common due to the lack of availability of this kind of data. It is important to handle the reference data taking into account the resolution capabilities of MODIS. Some authors resample the reference data in order to match the resolution of the coarser image. In the case of field work a fuzzy approach can be applied taking into account more than one class when the study area is fragmented or heterogeneous. A fuzzy approach adapted from Woodcock and Gopal (2000) was used in Mexico to prepare the reference data on the accuracy measurements (García-Mora and Mas 2011).

Another issue is to insure the independence between training data and reference data, especially in case of cross-validation. To do that, Friedl et al. (2010) stratified the training site database into 10 unique subsets using sites (not individual pixels) as the sampling unit in order to avoid spatial correlation in training and test data. Using this approach, 10 classifications were performed, each based on a unique combination of nine subsets to train the data, and using the remaining subset as reference data. In this manner, every pixel in the reference data were classified based on an independent training set and cross-validation is likely to have modestly lower predictive accuracy than the classifications used in operational generation of the product. However, as training sites consist of a polygon, delineated on higher resolution imagery, where the land cover is uniform and representative of a single IGBP class, accuracy assessment is also biased toward homogeneous areas, more likely correctly classified (Jung et al. 2006).

For fuzzy classifications, map and reference are expressed as continuous variables (e.g. fuzzy values estimating cover proportions and true proportions), and the most commonly used indices to assess the differences between estimated and true values are the coefficient of correlation, the root mean square error, both used for example by Braswell et al. (2003), and the bias (Muukkonen and Heiskanen 2007).

For hard classifications, indices are overall accuracy and Kappa coefficient. A few authors like Cohen et al. (2003) and Heiskanen (2008) provide per class accuracy indices such as user and producer accuracy or confusion matrices. Per-class accuracy is important data since the range at class-specific accuracies is generally large (see for example Friedl et al. 2010). When the entire confusion matrix is presented, it allows other users to compute new accuracy indices (Stehman and Czaplewski 1998).

6. Obtained results

In this section we review some characteristics of maps derived from MODIS data including the classification scheme, the number of categories, and the obtained accuracy.

Some of the studies are focused on a specific problem which requires only a limited number of land cover categories. For instance, some maps are binary with only two categories, such as irrigated/non irrigated (Ozdogan and Gutman 2008) or forest/no forest (Westra and De Wulf, 2007). Many maps have a reduced number of categories and few have more than 10 categories. For example, Borak and Strahler (1999) and Matsuoka et al. (2007) use 11 categories. General purpose land cover maps are often based upon the IGBP classification scheme which has 17 categories at the global level (Borak and Strahler 1999, Bagan et al. 2005, Friedl et al. 2010).

The obtained accuracies of the studies using MODIS vary highly from one case to another depending on the classification method, the number of categories, and the accuracy evaluation performed (Table 3). In most cases when using a broad classification the obtained results varied from 65 to 91% with large accuracy variation between classes. These results seem to be sufficient for many land cover mapping purposes. However, since different databases and approaches are used, it must be emphasized that reported accuracy measures are not comparable and should not be regarded as truly robust quantitative estimates. While a few accuracy assessment procedures used design-based sample schemes, others are based on a cross-validation – i.e. using several subsets of the training data or/and are biased to homogeneous areas giving likely optimistic accuracy estimates.

Table 3. Map characteristics obtained from MODIS. In the case where various approaches are used, the accuracy reported on the table is the higher obtained within the study.

Author(s) (date)	Classification scheme (# of categories)	Categories	Reference data used for accuracy assessment	Accuracy	Study area
Bagan et al. (2005)	IGBP (9)	1. Evergreen Needleleaf, 2. Evergreen Broadleaf 3. Deciduous Broadleaf 4. Grasslands5. Permanent Wetlands 6. Cropland 7. Urban and Built-up 8. Barren/Sparse 9. Water	Land Cover Dataset generated from Landsat and ASTER data	OA = 91.1% K = 0.8998 // per class accuracies from 75 to 100%	Various regions of China
Borak and Strahler (1999)^a	IGBP (11)	Evergreen needleleaf forest, Evergreen broadleaf forest, Deciduous broadleaf forest, Mixed forest, Woody savanna, Grassland, Closed shrubland, Open shrubland, Cropland, Urban and built-up, Bare soil and rocks	Map from Gap Analysis Program (GAP)	Overall accuracy 78.6%, with per class accuracies from 28 to 100%	Semi-arid region of Arizona, USA
Braswell et al. (2003)	Designed for the purpose of the study (3 to 4)	forest, cleared area, secondary (in one case old and young secundary)	Landsat TM classification	RMSE from 0.01 to 0.27	Amazon Region, Brazil
Carrão et al. (2008)	CORINE Land cover	artificial areas, rainfed crops, irrigated crops, broadleved forest, needleleaved forest natural grasslands shrubland, barren, water bodies	CORINE Land cover map (strata) + Landsat and air photographs	OA 90%	Portugal
Chen and Rao (2009)	Designed for the purpose of the study (4)	Degraded grassland, slat affected grassland, cropland, water	Landsat interpretation validated by field survey	OA 64.5 K 0.45	North China
Matsuoka et al. (2007)	IGBP	Evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, open shrublands, grasslands, croplands (including paddy), croplands (non-paddy), double-cropping fields (including paddy), double-cropping fields (non-paddy), irrigated fields, urban and built-up lands, snow and ice, barren, water bodies	A pixel-based comparison with MOD12Q1province-based comparisons for China using an existing 1-km resolution land cover map and ground-based census for China		East Asia
Miettinen et al. (2008)	Designed for the purpose of the study (1)	Wáter, mangrove, peatswamp forest, lowland forest, lower montane forest, upper montane forest, plantation, montane mosaic, lowland open, montan open, urban	XX Stratified random sampling against spot	Overall accuracy of 82%/K0.77	Western part of insular Southeast Asia
Ozdogan and Gutman (2008)	Designed for the purpose of the study (1)	Irrigated/non irrigated	Irrigation database	Correlation 0.88 between estimated and reported irrigation areas	USA
Shimabukuro et al. (2009)	Designed for the purpose of the study (1)	Burnt area	Landsat images classification, filed observations	Visual agreement, no quantitative data provided	Western Brazilian Amazonia
Sivanpillai and Latchininsky (2007)	Designed for the purpose of the study (1)	Reeds, redes and shrub mix, shrubs, sparse vegetation, Sandy soil, water.	Digital photographs, 203 ground data, landsat, interviews,	Error matrix OA 74% K 0.67	Amurdaya River Delta Uzbekistan
Tottrup (2007)	Designed for the purpose of study (3–4)	Mature forest, secondary forest, and non-forest (Fractional cover)	SPOT, ASTER and IKONOS lassifications	Mean absolute deviation between 14.6 and 21.6%	Mainland Southeast Asia uplands of Vietnam, Laos, Cambodia, and Thailand
Westra and De Wulf (2007)	Designed for the purpose of the study	Dry land, flooded land, irrigated rice, cultivation, and permanent water bodies	Cross validation	K = 0.98	Sahelian floodplains, Cameroon
Zhang et al. (2008a, 2008b)	Designed for the purpose of the study (8)	Crop, dense grasslands, sparse grasslands, Needleleaf Forests, Broadleaf Forests, shrublands, urban, water	Landsat classification	OA 75.5% PCA 35.4–96.3%	North China Plain
^aStudies based on MODIS like data.

Sedano et al. (2005) found that the global land cover classification derived from MODIS does not capture the spatial land cover variability of the study area located in Mozambique, providing a very broad classification in which most of the region is covered by only two classes (savanna and woody savanna). These two classes account for 59 and 36%, respectively of the study area, misrepresenting the heterogeneous distribution pattern of the land cover and land use types. Sivanpillai and Latchininsky (2007) pointed out that high commission errors resulted from misclassification of reeds mixed with shrubs class and shrubs classed as reeds. This could have resulted in over-prediction of the area designated as reeds. Such over-prediction of certain classes is also reported by Couturier (in press).

Additional research is therefore needed in order to assess misclassifications and biases in mapping complex landscapes such as mosaic landscape with patches showing several degrees of perturbation or recovering. A promising alternative to assess accuracy is quantifying the thematic and positional fuzziness of accuracy through a fuzzy logic approach (Couturier et al. 2009). This approach was used to carry out the accuracy assessment of the MOD12 IGBP map of the forest areas of Mexico (~1,000,000 km²) with a sample of field data (more than 24,000 conglomerates of four observation sites) and a land cover map (Couturier, in press). Based upon the most strict (without fuzzy tolerance) and flexible fuzzy option, the map accuracy is 54.6 and 65%, respectively.

7. Discussion and conclusion

MODIS presents unique features not present in other sensor: it has an ample variety of products including reflectance data, vegetation indices, biophysical variables, and thematic products which describe land cover characteristics. Moreover, the quality assessment per pixel information is an innovative feature potentially very useful to apply preprocessing procedures in order to select, transform, and/or filter the data. However, such information is not commonly applied: we only observed its use in preprocessing procedures to screen data or interpolate data gaps (Colditz, 2008). The limitation in its use is likely due to the fact that the information contained in the QA is not totally trustworthy to identify pixels with low quality as reported by various authors (Xiao 2005, García-Mora and Mas 2011). The improvements in the QA information contained within version 5 of the MODIS products will likely increase the incorporation of the QA information in the preprocessing of the data. MODIS also presents the advantage of offering a high temporal resolution data-set that potentially characterizes the phenology of the vegetation. It should be noted that the benefit of the high temporal resolution is often reduced by the lack of good quality data during the rainy season (Wessels et al. 2004).

Due to the large quantity of data and the presence of noise, preprocessing is a critical issue and is carried out using a large variety of approaches such as feature selection, feature extraction, and time series filtering and analysis. When using the QA per pixel information to select noise free images, an additional visual inspection is highly recommended to insure that the images that are chosen for analysis meet quality requirements.

The wide array of methods used to classify MODIS data are a reflection of the large number of MODIS products and the variety of mapping applications for this data. Among these methods, it is not clear which are more appropriate for specific applications and which can be applied with satisfactory results in other regions. However, DT is a classification algorithm which is repeatedly reported as very efficient in classifying MODIS data. Novel classification procedures which take into account the QA information in order to evaluate classification certainty can be envisioned.

Accuracy assessment is often not carried out in a systematic and rigorous manner with a probability sampling design and clear protocol for determining the reference land-cover classification of a sampling unit. Many reports of accuracy are likely optimistically biased. This can be an issue, particularly in fragmented diverse landscapes, where, due to its relatively coarse spatial resolution and bias of training data toward homogeneous regions, MODIS tends to misclassify heterogeneous landscape. However, results reported in the literature show that MODIS data can be successfully used to classify land cover/land use with a reasonable accuracy (>70%).

Due to its coarse spatial resolution MODIS cannot replace medium spatial resolution data such as that provided by Landsat and other higher resolution sensors, but it can present a useful complement, for example, filling the time-gap between two mapping efforts based upon high resolution imagery. It also provides the opportunity of mapping large areas using a single approach, which leads to a uniform final product and therefore represents an attractive alternative to map land cover at regional to global level.

As MODIS nears the end of life, the many years of work on sensors such as AVHRR, MODIS, and SPOT-vegetation by the involved agencies, together with the positive outcomes resulting from this collaborative effort will provide insight for the systems that will inevitably follow.

Notes on contributors

Tzitziki Janik Garcia-Mora has a bachelor in Biology and a Master in Conservation and Management of Natural Resources by the Universidad Michoacana de San Nicolás de Hidalgo. She has a Ph.D. in Geography with emphasis in Remote Sensing by the Universidad Nacional Autónoma de México. Her research interests are focused in land use/land cover mapping using moderate resolution data, classification of digital data, and accuracy assessment methodologies.

Jean-François Mas holds a Ph.D. from the University Paul Sabatier, Toulouse, France and he is a researcher at the Centro de Investigaciones en Geografía Ambiental (CIGA) at the Universidad Nacional Autónoma de México (UNAM). His research interests are focused on remote sensing applied to land use/land cover mapping and monitoring, land use/land cover changes modeling, conservation, community forestry, and remote sensing images analysis.

Everett Hinkley is the National Remote Sensing Program Manager for the U.S.D.A. – Forest Service, Geospatial Management Office located in Arlington, Virginia. In this position, Hinkley is responsible for providing national remote sensing program guidance and coordination to Forest Service field units throughout the USA, and serves as the Forest Service liaison to other federal and state agencies. His areas of interest include: land cover mapping, forest change detection, wildfire mapping and evaluation of unmanned aerial vehicles, and new sensors for remote sensing applications. He holds a Master of Science degree in mapping and geographic information systems from Ohio State University and a Master of Science degree in Geology from Miami University of Ohio.

Acknowledgements

This review was carried out under the projects Un sistema de monitoreo de la deforestación en México (Fondo Sectorial de Investigación para la Educación SEP-CONACYT, clave 47198) and Evaluación del sensor MODIS para el monitoreo anual de la vegetación forestal de México (Fondo Sectorial para la Investigación, el Desarrollo y la Innovación Tecnológica Forestal CONACyT-CONAFOR, clave 14741).

The first author is grateful for the financial support provided by the PhD scholarship funded by CONACyT on the 2007 call. Final writing of the manuscript has been carried out during a sabbatical stay of the second author at the University of California-Santa Barbara with the support of Consejo Nacional de Ciencias y Tecnología (CONACyT) and the Dirección General de Asuntos del Personal Académico (DGAPA) at the Universidad Nacional Autónoma de México.