Characterizing temporal vegetation dynamics of land use in regional scale of Java Island, Indonesia

Abstract

Improving the understanding of land use and land cover is a major research challenge for the human-environmental sciences and is essential for many aspects of global environmental research. Considering seasonal vegetation dynamics or phenological dynamics in multi-year series leads to a broader view of land use and land cover. This study is based on the hypothesis that a pixel representing a complex but consistent land use has a typical, distinct and repeated temporal pattern of vegetation index inter-annually, which can be used as characteristic signatures for land use classification. Considering the seasonal events and climatic variability in Indonesia, we characterized the temporal vegetation dynamics of long-term land use by using multi-temporal Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index 16-day composite data from 2001 to 2007, and then generated a land use map using those characteristics. Accuracy assessment of the results showed the need to evaluate such methods for land use types that do not have a consistent yearly pattern. On the other hand, the identification of the intensive agriculture lands, such as paddy rice and upland, was satisfactory. Although the mixed pixel issue is quite problematic when using MODIS data, the results indicate that MODIS data offer great promise for characterizing seasonal as well as multi-year variation at large scales. Indeed, the methodology proposed in this research distinguished among many specific land use classes based on temporal land cover information properties. Characterization of temporal vegetation dynamics patterns would provide sufficient, significant and useful information regarding the patterns of land use; consequently it should be possible to consider the actual, subtle nature of inter-annual land use change as well as overall land use.

Keywords:

• temporal vegetation dynamics
• land use
• MODIS EVI
• Java Island

Introduction

Information on global land use and land cover (LULC) is essential for many aspects of global environmental research (Foley et al. 2005; Turner, Lambin, and Reenberg 2007). Earth observation from satellites provides spatial measurements of biophysical surface attributes, such as vegetation cover, biomass and landscape heterogeneity (Lambin, Geist, and Rindfuss 2006). Analysis of multi-year time series of those attributes and their seasonal evolutions has led to a broader view of LULC (Ramankutty et al. 2006; Philippon, Jarlan, Martiny, Camberlin, and Mougin 2007; Erasmi, Propastin, Kappas, and Panferov 2009).

Commonly, land surface is characterized by distinguishing among different land cover types (Jansen and Gregorio 2002). Land cover addresses the layer of soils and biomass, including natural vegetation, crops and human structures that cover the land surface. Land cover is thus directly observable, both in the field as well as from remote sensing images. Land use, in contrast, refers to the purposes for which humans exploit the land cover (Lambin et al. 2006) and includes the land management practices. In many cases, land cover may drastically change without any change in land use. The mixed use of land cover and land use in maps and models leads to inconsistencies in land change assessments (Verburg, Van De Steeg, Veldkamp, and Willemen 2009).

The monitoring of land surface continuously in space and time, at both seasonal and inter-annual scales, allows characterization of the temporal vegetation dynamics of land use (Lambin, Geist, and Lepers 2003; Verburg et al. 2009), in other words, distinguishing the actual type of land use based on land cover dynamics. In tropical regions, for example Java Island, Indonesia, such discrimination is an important issue due to the sequential land use patterns that exist in the field. An area classified as a paddy field might undergo a sequence of covers through the year, such as (1) paddy–bareland–secondary crops–bareland, (2) paddy–bareland–inundated–paddy–bareland–secondary crops and (3) paddy–bareland–secondary crops–bareland–inundated–paddy, where the sequence is repeated year after year following the seasons. These temporal covers are characteristic of the cropping systems in some regions of the tropics. In this example, attributes of the land surface are: vegetation (paddy), soil (bareland) and water (inundated), which is defined as land cover. On the other hand, the land use type, which describes the purpose of land, is the paddy rice field. Similar distinctions are important when considering plantation crops, forest and other types, even if there is a difference in time scale of the trajectory cycle (Verburg, Veldcamp, Willemen, Overmars, and Castella 2004).

Vegetation dynamics have attracted attention as a means for better understanding land characterization (Verburg et al. 2009). Understanding vegetation dynamics requires the characterization of vegetation changes at different temporal scales (Lunetta, Knight, Ediriwickrema, Lyon, and Worthy 2006; Galford et al. 2008; Martinez and Gilabert 2009). The appropriate temporal changes are characteristic of biologically complex ecosystems and are known as phenological dynamics of terrestrial ecosystems, reflecting the response of the earth's biosphere to annual dynamics of the earth's climate and hydrologic cycle (Zhang et al. 2003). Moreover, vegetation phenology, as used and studied by remote sensing, refers to the relationship between climate and periodic development of photosynthetic biomass (Li and Kafatos 2000; Jakubauskas, Legates, and Kastens 2002; Kang et al. 2003).

Characterization of vegetation dynamics has often been made by using vegetation index values (Zhang et al. 2003; Knight, Lunetta, Ediriwickrema, and Khorram 2006; Lunetta et al. 2006). The temporal dynamics of those index values are useful for distinguishing land surface conditions by differentiating among vegetation types and their distributions (DeFries, Hansen, and Townshend 1995; Hansen, DeFries, Townshend, and Sohlberg 2000; Lloyd 2000). The vegetation indices relate to the amount of green leaf biomass. They are computed by algebraic combination of two or more spectral bands and are designed to enhance the contribution of vegetation properties (Huete, Liu, Batchily, and Van Leeuwen 1997). The normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) are commonly used to measure reliable spatial and temporal inter-comparisons of terrestrial photosynthetic activity and canopy structural variations (Huete et al. 2002). Whereas the NDVI is chlorophyll sensitive, the EVI is more responsive to canopy structural variations, including leaf area index, canopy type, plant physiognomy and canopy architecture (Gao, Huete, Ni, and Miura 2000). The 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 atmosphere influences (Huete et al. 2002).

Some previous studies (Goward, Tucker, and Dye 1985; Justice, Townshend, Holben, and Tucker 1985; Loveland, Merchant, Ohlen, and Brown 1991; Schwartz 1999) emphasize that time-series images are needed to recognize the vegetation phenology. In addition, DeFries and Townshend (1994) also pointed out that multi-year data sets would improve the capability to (1) characterize portions of the study area, such as some sites with high inter-annual variability in vegetative activity and (2) ensure that a classification result and land cover change analysis were not the product of short-term fluctuations in vegetative activity. The uses of multi-year data sets are necessary to develop methodologies that utilize information on inter-annual variations to increase the accuracy of the land surface characterization. Initially this work relied on the use of coarse spatial resolution (>1.0 km²) Advanced Very High Resolution Radiometer data but, with the availability of the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on board Terra and Aqua Satellites after launch in December 1999 and May 2002, respectively, the time-series data analysis can be applied to 250 m pixels (Privette et al. 2002; Lobell and Asner 2004; Xiao, Liang, Wang, Song, and Wu 2009).

Extending the themes mentioned above, this study is based on the hypothesis that a pixel representing a complex but consistent land use has a typical, distinct and repeated temporal pattern of vegetation index inter-annually, which can be used as a characteristic signature for land use classification Consequently, a land use map could be drawn from the temporal pattern of the time-series data. The approach in this article is different from the previous approach used to produce the available land use/cover maps of Indonesia, since we consider a variability of temporal vegetation dynamics in land use classes. The consideration of temporal vegetation dynamics of land use could also be used to support the land use change analysis at a regional scale. This article will investigate temporal vegetation dynamics of long-term land use considering seasonal events and climatic variability in a tropical region. Then, we will generate a land use map using those temporal patterns.

Methods

Study area

This study encompassed all of the island of Java which is located on the southern rim of the Indonesian archipelago (upper left corner: 5°52′S 105°04′E; lower right corner: 8°47′S 114°36′E) and has an area of 132,792 km² (Figure 1). The island is administratively divided into six provinces, namely Banten, Jakarta (special capital region), West Java, Central Java, Yogyakarta (special region) and East Java.

Figure 1. The study area.

The climate of Java reflects the region's location in the tropics, with temperature and humidity remaining constant throughout the year, and rainfall varying between the seasons. Most areas receive 2000–4000 mm of rain annually. The driest areas, with less than 1500 mm, are along the north coast reaching a minimum of less than 1000 mm around Asembagus in the northeast. There are dry intermontane basins with less than 2000 mm around Bandung and Garut in West Java; around Surakarta and Yogyakarta in the centre; and around Madiun, Kediri and the Malang basin in East Java. The wettest areas of Java are mostly inland. Over 4000 mm fall on the mountains bordering the Sunda Strait and on or close to the upper slopes of most peaks in the Central Volcanic Mountains. These include the northern slopes of Ragajembangan, which is the wettest place in Java with more than 6000 mm (Land Resources Department/Bina Program 1990). In general, from November until March, the zone of heavy rainfall covers most of Java, while from May to September; the zone of rainfall is limited to the southeast, leaving the eastern parts of the island dry.

According to the BPS-Statistics of Indonesia (2003), about 70.62% of Java is considered to be agricultural, with land use as follows: paddy fields, mixed gardens, uplands/drylands, fishponds and plantations, with as much as 5.43% of the area covered by settlements. Furthermore, Ministry of Public Works (2002) stated that paddy rice accounts for about 2,807,770 ha or 21.14% of Java's total area, while uplands/drylands and non-irrigated areas cover about 6,352,850 ha or 47.84% of the total area.

Forests in Java cover 22.98% of Java's area, consisting 1,246,728 ha of production forest (9.39%), 550,849 ha of protection forest (4.15%) and 529,702 ha of conservation forest (3.99%) (Department of Forestry 2002). According to the Act of the Republic of Indonesia No. 41 of 1999 on forestry, a conservation forest is a forest area with specific characteristics, having the main function of preserving plant and animal diversity and its ecosystem. Meanwhile, a protection forest is a forest area having the main function of protecting life-supporting systems for hydrology, preventing floods, controlling erosion, preventing sea water intrusion and maintaining soil fertility (State Gazette of the Republic of Indonesia 1999).

Data

Satellite image

The MODIS product used in this study is the Vegetation Indices (VI) Composite 16-day Global 250 m SIN Grid V005 or MOD13Q1 product, which provided the needed vegetation phenology data. In addition, the product had already been systematically corrected for the effects of gaseous and aerosol scattering. The MODIS EVI is embedded in the MOD13Q1 product. The MODIS Land Discipline Group (MODLAND 2010) developed the EVI for use with MODIS data following this equation:

(1)

where, , and are the remote sensing reflectances in the near-infrared, red and blue, respectively, L is a soil adjustment factor and C ₁ and C ₂ describe the use of the blue band in correction of the red band for atmospheric aerosol scattering. The coefficients, C ₁, C ₂ and L, are empirically determined as 6.0, 7.5 and 1.0, respectively. G is a gain factor set to 2.5 (Huete et al. 1997).

The EVI data were developed in the above form (Equation (1)(1)) in order to optimize the vegetation signal with improved sensitivity in high biomass regions. The EVI also minimizes atmospheric influences with the ‘aerosol resistance’ term which uses the blue band to correct aerosols influence in band red (Huete et al. 1997). In this study we used the MODIS EVI data sets which were acquired from January 2001 to December 2007 and captured 161 time series with the interval time 16 days. The data were obtained at no cost from USGS Land Processes Distributed Active Archive Center (USGS LP DAAC 2009a). The study area is covered by two MODIS tiles: h28v09 and h29v09.

MODIS EVI data were extracted from the MODIS VI product (MOD13Q1) using the MODIS Reprojection Tool (USGS LP DAAC 2009b) and the selected output format was GeoTIFF and coordinate system was geographic coordinate systems on datum World Geodetic System of 1984. The tiled MODIS EVI data were then mosaicked and clipped to cover Java for each composite period and then sequentially stacked to produce the time-series data set.

To examine the details of the temporal pattern of the MODIS EVI, we used a finer spatial resolution Landsat 7 Enhanced Thematic Mapper Plus (ETM+) as well as reference data derived from the ground survey points. The study area is covered by 18 scenes of Landsat images, collected from 2001 to 2007.

Constraints of the MODIS EVI product

Although MODIS data have some advantages in providing basic information related to vegetation phenology, time series of these data inevitably contain disturbances caused by atmospheric variability (Huete and Liu 1994) and aerosol scattering (Xiao, Braswell, Zhang, Boles, and Frolking 2003). However, improving aerosol correction at the surface reflectance level (MOD09) and a new filtering scheme in the VI algorithm implemented in collection 5.0 have positively impacted the MODIS VI (Didan and Huete 2006).

Moreover, to get a greater percentage of clear-sky data, the maximum value composited method is applied to the MODIS VI and is combined with the MODIS BRDF (bidirectional reflectance distribution function) or MOD43 product to generate the 16-day composite MODIS VI (MOD13Q1 product) (Huete, Justice, and Van Leeuwen 1999). Nevertheless, if the composite period is too long, the land surface does not remain static; and if it is too short, the atmospheric disturbance cannot be removed effectively (Lu, Liu, and Liang 2007); consequently, there are some residual errors. Such noise degrades the data quality and introduces considerable uncertainty in temporal sequences, confusing the analysis of temporal images sequences by introducing significant variations in the EVI time-series data. Therefore, noise reduction or fitting a model to observed data is necessary before phenological stages can be determined (Sakamoto et al. 2005; Lu et al. 2007).

Various strategies for image preprocessing including smoothers have been applied, such as polynomial and median filters (Van Dijk, Callis, Sakamoto, and Decker 1987), a moving window to select the local maximum VI (Viovy, Arino, and Belward 1992), temporal window operation (Park, Tateishi, and Matsuoka 1999), logistic curve fitting (Zhang et al. 2003), the asymmetric Gaussian function fitting approach (Johnson and Kasischke 1998), principal component analysis (Li and Kafatos 2000) and the Savitzky–Golay filter approach (Chen et al. 2004). Several spectral-frequency techniques have also been used, including Fourier-based fitting methods for separating the high-frequency components of noise and the low-frequency components of seasonal changes of VI (Sellers et al. 1994; Azzali and Menenti 2000; Lunetta et al. 2006) and, most recently, wavelet decomposition, which has been used to denoise and to build a vegetation signatures through feature extraction of MODIS multi-temporal (Bruce, Mathur, and Byrd 2006), to characterize crop phenology (Sakamoto et al. 2005), to detect temporal changes of annual flooding (Sakamoto et al. 2007), to determine changes in the expansion and intensification of crops (Galford et al. 2008) as well as to investigate vegetation dynamics in terrestrial ecosystem (Martinez and Gilabert 2009).

Image processing

DeFries and Townshend (1994) explained that supervised classification based on temporal change differences faces several issues. First, data sets of a sufficient quality must be available in order to carry out the analysis. Second, in order to carry out a classification each land cover/use type to be included, each category must exhibit a distinct temporal signature based on the temporal pattern change. Third, there must be a systematic procedure for validation of the result. The image processing was designed to addresses these issues and Figure 2 illustrates the process of image data handling and analysis. We used wavelet transforms to reduce the residual noise and reconstruct the temporal pattern of image data. The denoised data were then clustered into several significant temporal patterns, and then the supervised classification was performed based on those temporal patterns.

Figure 2. Image preprocessing and analysis flow. LU, land use; ETM+, Enhanced Thematic Mapper Plus; MODIS, Moderate Resolution Imaging Spectroradiometer; EVI, enhanced vegetation index.

Image transformation using wavelet function

The wavelet transform can decompose signals in time-frequency space in order to identify and reduce the overall noise as well as to maintain useful information in the time-series data. Furthermore, Sakamoto et al. (2005) concluded that the wavelet transform gave the best time profile of EVI for the purpose of detecting phenology stages, and they proposed the Wavelet-based Filter for determining Crop Phenology to remove high-frequency noise in crop phenology data. Galford et al. (2008) has applied the Wavelet-based Filter for determining Crop Phenology in areas influenced by the tropical rainy season in southwestern Brazilian Amazon. They have demonstrated the stability of the wavelet approach to filter out noises in EVI time-series data. In this study, the wavelet transform is applied to a time series of MODIS EVI data in order to generate a land use map based on the annual pattern of land use types on Java Island, Indonesia.

The wavelet transform decomposes a signal into different scales by successively translating and convolving the elements of a high-pass and low-pass scaling filter associated with the mother wavelet. These filters retain the large- and small-scale components of the signals also known as the approximation (A) and detail (D) series, respectively. For most practical applications, the discrete wavelet transform, which analyses signals over a discrete set of scales, follows a hierarchical algorithm for such decomposition, known as the pyramid algorithm (Mallat 1989).

In the first level of the decomposition f(t) = A ₁ + D ₁ the signal has a low-pass filtered component, A ₁, and a high-pass filtered component, D ₁. In a second step, the approximation A ₁ is split as A ₁ = A ₂ + D ₂ and so on. The relationship D_j = A_{j– 1} – A_j gives us information about the portion of the signal that can be attributed to variations between the scales [j–1, j]. Figure 3 shows an example of such decomposition corresponding to the paddy field in Java Island.

Figure 3. The structure of multilevel wavelet decomposition tree for the EVI time series.

In order to analyse the fluctuation of the EVI value of each land use type, we used the Coiflet mother wavelet because this wavelet shape gives remarkably good results in determining phenological stages as pointed out by Sakamoto et al. (2005). The order in the wavelet function is a measure of the wavelet's smoothness, where a larger order produces a smoother wavelet (Burke 1994). In this processing, order 2 of the Coiflet function was used since the trend of that order is similar to the trend of the original data. This analysis was performed using MATLAB through the 1D multisignal wavelet analysis function (MathWorks 2009).

The MODIS EVI data preprocessing was conducted to provide a filtered data set to support multi-temporal (phenological) analysis. Figure 4 shows that the wavelet transform application effectively deals with discontinuous/sharp spike data (resulting from some unspecified error, i.e. residual errors).

Figure 4. Illustration of the filtering result of Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI) by wavelet transform.

Land use classification

Corresponding the previous land cover map

A land cover map provided by the Department of Forestry (DoF) (2008) was used as a reference to select the pixels used to develop the time-series EVI. We used the land cover data set recognizing that the actual types of land use are distinguished based on land cover; therefore, each land cover might be attributed to a number of different land use types.

The DoF land cover map was generated based on the interpretation of Landsat 7 ETM+ acquired in 2005/2006 and performed by visual classification. This classification was performed by digitizing polygons based on image interpretation with minimum mapping unit (MMU) defined by means of the smallest polygon identified (Rocchini, Perry, Salerno, Maccherini, and Chiarucci 2006). The land cover map was produced at a scale of 1:250,000 with a MMU of 40 ha.

The DoF is concerned primarily with natural vegetation categories. Using the thematic baseline forestry and designation of forestry-aquatic area maps for 30 provinces at the scale of 1:250,000 (Department of Forestry 2008), delineated polygons for those classes could be labelled with two codes, representing categories undisturbed by human activities (primary type) and those disturbed (secondary type). Finally, the satellite data were classified into 23 land cover types (Table 1). The 23 categories seem sufficient to represent land cover types in Indonesia. This is near the maximum number of categories that could be discriminated effectively; identification of other categories would likely be difficult to achieve.

Table 1. The 23 categories of land cover type by the Department of Forestry

Number	Land cover type	Description
1	Primary dry land forest	Forest in the lowlands, hills and mountains, where there was no logging activity.
2	Secondary dry land forest	Forest in the lowlands, hills and mountains, where logging had occurred.
3	Primary swamp forest	Forest in swampy areas, including peat swamps not yet showing signs of logging activity.
4	Secondary swamp forest	Forest in swampy areas, where logging had occurred.
5	Primary mangrove forest	Mangrove forests and palm around the coast that have not been logged.
6	Secondary mangrove forest	Mangrove forests and palm that have been logged, as revealed by the pattern of grooves.
7	Bush/slash	Natural plants with multiple stems and canopy cover is greater than agriculture but is not a forested area.
8	Swamp	Vegetation or swamp area that was not forested swamp.
9	Swamp bush	Bush/shrub converted from forest in swampy areas, including bush in swamp areas that will be converted into pond, forest plantation or plantation forests.
10	Savannah	Natural vegetation non-forest, grassland with few trees.
11	Industrial forest plantation	Large-scale of agriculture in which a crop plantation is planted.
12	Plantation	Plantation area covered by a plantation crops. Plantation areas that are bare land or shrubs are classified according to those land covers (bare land or bush).
13	Dry land agriculture	Agricultural activities on dry lands, that is, mixed garden and uplands, include agricultural crops, fruit trees.
14	Mixed dry land agriculture	Agricultural activities on dry land, which is mixed with bushes or shrubs.
15	Rice land	Agricultural area in wetlands, planted with rice and characterized by a ridge pattern.
16	Fishpond	Fish farming area along the shore.
17	Bare land	Area without vegetation (rock outcrops and mountain peaks, volcanic craters, sand).
18	Transmigration	Settlement area including agricultural land, when not classified as agricultural land.
19	Mining	Area without vegetation, open area.
20	Airport area	Airport area identified by a runway.
21	Built-up area/housing	Settlement areas, rural areas, ports, airports and so on industries.
22	Water body	Water features, including sea, rivers, lakes, reservoirs, coral reefs and seagrass.
23	Cloud	Cloud cover over an area.

In addition, although the land cover map provided by the DoF was based on the Landsat ETM+ imagery, for which one pixel of image represents approximately 0.1 hectares of the earth surface (spatial resolution of Landsat about 30 m), the product was generalized based on the smallest polygons that could be identified by visual interpretation. As a result, the DoF classification resolution was effectively courser than a spectral classification using the Landsat data. This classification indicated that approximately 63–83% of 250 m pixel grids contained homogeneous cover types.

In the context of temporal pattern of land use, we need to resolve the differences among the classes in the data set in order to provide a consistent definition for pixel-based analysis. The process of reassigning land use categories based on the knowledge of the phenological pattern of each land use is shown in Table 2.

Table 2. The reassigning of land use categories

	Land cover category			Land use category
Number	Department of Forestry (DoF)		Number	Reassignment of land use category	Number of pixel samples
1	Primary dry land forest
2	Secondary dry land forest	→	1	Natural forest	35,835
3	Primary swamp forest
4	Secondary swamp forest
5	Slash
6	Swamp bush	→	2	Bush	3,240
7	Swamp
8	Savannah
9	Primary mangrove	→	3	Mangrove	25,930
10	Secondary mangrove
11	Industrial forest plantation	→	4	Industrial forest plantation	109,498
12	Plantation		5	Plantation	30,490
13	Rice land	→	6	Paddy rice field	111,729
14	Mixed dry land agriculture	→	7	Mixed garden	41,586
15	Dry land agriculture	→	8	Upland	82,373
16	Built-up area/housing
17	Transmigration			Non-vegetation
18	Airport area		9	Built-up (settlement)	29,855
19	Mining	→	10	Fishpond	9,028
20	Bare land		11	Water body	2,262
21	Fishpond
22	Water body
23	Cloud

A number of pixels corresponding to each of the land use categories were selected randomly from the parent category (11 categories). More than one typical temporal pattern might be identified from each land use category, which represents a specific class of land use. The number of pixel samples to be taken for each category is given in Table 2.

Training sample sets

Previous studies (Tucker, Townshend, and Goff 1985; Townshend, Justice, and Kalb 1987; Knight et al. 2006) have derived land cover maps from NDVI at continental and regional scales based on phenological differences among cover types as reflected in temporal profiles of NDVI. Similar assumptions were used in this study, that is, that locations displaying similar temporal patterns are inferred to have a similar land use type. The land use types are in turn inferred to have relatively homogeneous characteristics such as fractions of vegetation covers and amount of biomass. The use of the long time series (2001–2007) of EVI was based on the hypothesis that the land use can be distinguished based on land cover dynamics because consistent land use has a typical, distinct and repeated temporal pattern of vegetation index inter-annually.

Distinguishing among temporal patterns in order to recognize land use types was accomplished using the k-means clustering method based on Euclidean distance in an EVI-space in which each EVI image provides one dimension of the cluster space, analogous to spectral clustering. Distinct clusters then represent different temporal patterns–temporal signatures as opposed to spectral signatures. The resulting patterns were then used to discriminate among pixels selected from the reassigned land cover map described in Table 2. The clustering method subdivides data sets into k clusters through an iterative process in order to optimize a criterion function (Jain, Murty, and Flynn 1999). Clustering yields a number of significant patterns which correspond to 11 categories (illustrated by Figure 5).

Figure 5. Clustering process to differentiate the variation of the temporal patterns. LU, land use; EVI, enhanced vegetation index.

The number of significant patterns, each representing a specific class of land use, was determined based on a statistical measurement of how similar that pattern is to patterns in its own cluster compared to patterns in other clusters, and then combining clusters as needed. Furthermore, the significant pattern was identified based on its temporal profile, a corresponding ground survey data set and other image data such as Landsat or detailed information which was obtained using the ground survey as a reference. Figure 6 illustrates the use of multiple-date Landsat images in conjunction with a 7-year MODIS EVI pattern. Then, based on the information of those images and supported with local knowledge, the type of land use of that EVI pattern is identified. The same method was applied to other EVI patterns in study area. Once land use types were associated with their EVI patterns, we used those identified pixels as a signature sample (training sets) for supervised classification.

Figure 6. (a) The multiple dates of Landsat Enhanced Thematic Mapper Plus (ETM+) in conjunction with 7-year data sets; (b) temporal pattern of Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI) profile. Note that the EVI profile above is assigned to double a cropping system in an agriculture area. The EVI value of the red rectangle of each image on Figure 6a is located by the red point on Figure 6b.

Supervised classification

Several specific land use types that broadly represent the major land use types on Java Island were distinguished according to their EVI temporal patterns. Then, to explore the feasibility of land use classifications for the whole study area (landscape scale), we ran a maximum likelihood supervised classification. The maximum likelihood procedure also classifies each pixel based on its temporal pattern to assign the land use type. The maximum likelihood algorithm (Richards and Jia 1999) belongs to the class of parametric classification methods. The mean vector and the variance–covariance matrix for each land use type are estimated using a number of pixels from training sets. Then, using the maximum likelihood rule, the multidimensional space is partitioned into subspaces each uniquely associated with one land use type (Keuchel, Naumann, Heiler, and Siegmund 2003).

The whole of the land use type in Java Island was then classified according to the temporal pattern properties of each pixel. Thus, if a pixel falls within the subspace associated with land use type c_i , it is labelled as type c_i . If the pixel falls within the subspace associated with land use type c_j , it is labelled as that land use type, c_j .

Accuracy assessment

The creation of a reference data set in order to assess the accuracy of the classification result is one important step in this land use study. A complicating factor in assessing the accuracy of a classification derived from 250 m image data in the Java Island was the disparity between the image pixel size and the average patch size of the landscape. A previous classification of the land cover by Ministry of Environment (MoE) (2008) using Landsat ETM+ had an effective MMU of 2.4 ha. Overlaying a 250 m grid (corresponding to pixel size of MODIS) over this classification revealed that approximately only 20% of the grid contained homogeneous cover types at that MMU. The same issue was mentioned by Atkinson, Cutler, and Lewis (1997), that one problem of moderate spatial resolution imagery is that the intrinsic scale of spatial variation in many land covers is finer than the scale of sampling imposed by the image pixels. The result is that most image pixels contain a mixture of land cover classes, referred to as sub-pixel mixing (Cracknell 1998).

The creation of a reference data set to assess the accuracy of a MODIS classification by providing sub-pixel land cover information was described by Knight et al. (2006). They performed an accuracy assessment by sub-pixel reference information in which the USGS Digital Orthophoto Quadrangles were selected as the base data source. The same approach was carried out to our classification result.

To capture the variability, the study area was stratified by landscape, from mountains to coastal, and 2110 locations were selected based on random sampling. Each location was overlaid with a 250 m fishnet corresponding to the pixels of the MODIS EVI data, then 10 by 10 cells were selected as our plot size (Figure 7a), and then each cell was overlaid with a 100-point dot (Figure 7b). Each point was then associated with a specific land use based on its location on the Landsat image, a high-resolution image (Google Earth) and the reference data derived from the ground survey points. Whenever possible we provided more detailed land use information (specific classes). In the case of agricultural land, we were supported by a local investigation and other secondary data such as irrigated area, irrigation infrastructure, harvested area, data production and annual rainfall data sets. In situations where detailed information of land use was difficult to ascertain, we provided only a land use category (the parent category) label. This procedure provided the percentage of each land use class within the corresponding MODIS pixels; cells with over 70% of a single land use were taken to be homogeneous and used as a reference for accuracy assessment. The final number of samples in the reference data was 18,505 pixels and 10,000 pixels on the category level and the specific class level, respectively.

Figure 7. (a) A 250 m fishnet (Moderate Resolution Imaging Spectroradiometer (MODIS) pixel size) overlaid with Landsat Enhanced Thematic Mapper Plus (ETM+) and (b) cell on high-resolution image with 100-point dot grid overlay.

In this study, the accuracy of the classification was assessed using several accuracy parameters: total accuracy, user's accuracy and producer's accuracy. These accuracy parameters involve error matrices describing the number of sample units/points correctly classified, misclassified or not classified which can be characterized as commission error (user's accuracy) and omission error (producer's accuracy), respectively (Congalton and Green 1999).

On specific types of land use, for example paddy rice, we performed a statistical comparison with the statistics data at the village level (BPS-Statistics of Indonesia 2003) in order to assess the classification result for these specific classes. This was possible due to the availability of data regarding irrigated area, irrigation infrastructure, harvested area, field rate and data production for each village. The statistics’ data of 24,962 villages (1988 sub-districts) in Java Island were used. These data provided more detail about the accuracy, because they were based on a local census. However, this level of detail was only attempted for the land use type of paddy rice.

Results and discussions

Identification of temporal vegetation dynamics

As mentioned earlier, we identified each land use type from a long time series of land surface data sets. The seasonal pattern can be clearly understood from the temporal pattern of EVI. For example, the paddy field area is submerged during the rainy season and then vegetation starts to emerge.

In the upland category, the number of cultivation cycles (seedling to harvesting) that occur yearly depends on the availability of water. The fluctuation of average EVI in the upland area was differentiated into three significant patterns (Figure 8a). These were identified as mixed garden (LU-1, upland type 1), upland mixed bush (LU-2, upland type 2) and upland with intensive agriculture farming (LU-3, upland type 3), respectively. In fact, intensive farming activity also existed in the upland type 1, but the area was relatively small. Another important feature of upland type 1 was that agriculture made extensive use of mixed cropping and intercropping of annual crops. These practices provide farmers an opportunity for harvesting diverse crops from the same land. On the other hand, the upland type 3 clearly shows the annual vegetation growth cycle, representing intensive cropping with multiple harvests in the upland area.

Figure 8. Identification results of land use type, which can be differentiated from the temporal patterns of enhanced vegetation index (EVI). Group of land use category: (a) upland, (b) paddy field, (c) forest, (d) mangrove, (e) plantation and (f) non-vegetation.

The EVI pattern of a paddy field segmented into seven significant types (Figure 8b); nevertheless, the successfully classified patterns of a paddy rice cropping system were only represented by paddy type 1 (LU-4), paddy type 2 (LU-5), paddy type 3 (LU-6) and paddy type 4 (LU-7), while paddy type 5 (LU-8) corresponded to the upland which was planted at the end of a rainy season. Paddy type 6 (LU-9) was identified as fish/shrimp pond and became land during the long dry season. Paddy type 7 (LU-10) was classified into upland with intensive cropping.

Based on these patterns, the rice paddies could be sorted into several rice cropping systems. The farmers apply single, double or triple cropping patterns depending on the availability of water, either by irrigation or by rain water. In the irrigated areas, paddy rice can grow during dry season, so that rice crops can be grown 1 to 3 times in this land within a year. For that reason, the crops which are planted in irrigated areas display continuous EVI fluctuation, even in the dry season.

The EVI pattern of forest over 7 years was separable into three types (Figure 8c). The examination of each pattern indicated that discrimination among them was related to differences in the forest composition. Forest composed of a homogeneous distribution of trees showed a different pattern than that of a heterogeneous mixture of trees. The forest type 1 (LU-11), type 2 (LU-12) and type 3 (LU-13) of the discriminated patterns were identified as forest mixed with bush, bush–shrub mixed grass and dry–land forest, respectively.

Leaf phenology in natural forest has been studied in many tropical regions, such as the Amazon forest (Nepstad et al. 1994; Jipp, Nepstad, Cassel, and Carvalho 1998; Asner, Townsend, and Braswell 2000), South America (Xiao, Hagen, Zhang, Keller, and Moore 2006) and Borneo forest (Shimada, Takahashi, and Limin 2006). Field measurements in Borneo forests have shown that the canopies respond to a seasonal dry period, with decrease of NDVI of 3–5%. On the other hand, in the tropical forests of South America, the maximum EVI typically occurred between the late dry season and the early rainy season. Moreover, Xiao et al. (2006) pointed out that leaf phenology in the tropical forests is not determined by the seasonality of precipitation, but rather may be driven by the availability of solar radiation. Although the above-referenced studies ascribed different causes for the trends, they successfully demonstrated the leaf phenological phenomena in tropical forest areas.

Similar to the pattern of forest described above, the EVI patterns of mangrove also indicated a stable EVI pattern that did not show any clear periodic fluctuation. Besides the density of vegetation, the variety of their patterns was affected by a number of geophysical factors to which the mangroves have adapted, such as the changing of salinity concentrations, sediment influxes and regular inundation by water (Bird and Ongkosongo 1980). Based on the different patterns, the mangrove could be classified into three classes (Figure 8d). Mangrove type 1 (LU-14) represents a heterogeneous patch of mixed mangrove species, examples of which are distributed around Southern Java. Mangrove type 2 (LU-15) represents mangroves that are significantly influenced by water fluctuation. These areas are typically inundated longer. As with the zonation of mangroves suggested by Watson (1928), mangrove type 2 was dominated by Avicennia spp., Sonneratia spp. and Bruguiera gymnorrhiza, while mangrove type 3 (LU-16) was identified as a mixed garden with bush.

In the case of plantations, different temporal patterns of EVI not only follow the growing season, but also depend on their maintenance, such as weeding and replanting. The plantations are widely distributed in Java: industrial forest plantations, rubber and oil palm plantations. The temporal pattern EVI of plantations was grouped into four types (Figure 8e). Plantation type 1 (LU-17), type 2 (LU-18) and type 3 (LU-19) were identified as rubber plantations, oil palm plantations and industrial forest plantations, respectively. The industrial forest plantation was related primarily to teak plantations which exhibit leaf falls under water stress during a dry season, that is, between June and October. Meanwhile, plantation type 4 (LU-20) did not fit a plantation pattern, but is more characteristic of an intensive cropping system similar to upland type 3 (LU-3). Thus, plantation type 4 was identified as uplands in the plantation area.

The non-vegetation class is considered to represent water bodies, settlement and open areas. The clustering method differentiated the temporal patterns into five different types (Figure 8f). Non-vegetation type 1 represented lands which are inundated during the rainy season, but start to be planted in the last days of the rainy season and turn into dry land during the dry season. Thus, we classified that type into rice-paddy in lowland areas (LU-21). This kind of land use is distributed around Eastern Java. Meanwhile, non-vegetation type 2 was classified as a water body, such as lake or river (LU-22). Non-vegetation type 3 was identified as settlement or mining area (open area) (LU-23). Similar to type 1, the fluctuations of EVI in non-vegetation type 4 (LU-24) and type 5 (LU-25) are characteristic of vegetation growth or cropping activities. Therefore, the LU-24 was classified into paddy rice with a double cropping system per year, and LU-25 was classified as a mixed cropping which becomes barren during a dry season.

From the identification of the temporal patterns of EVI, 25 types of land use were distinguished which broadly represent the major land use types of Java Island. These are listed in Table 3.

Table 3. The list of land use types as a summary of the analysis of EVI pattern

Number of LU type	Type of land use based on the temporal pattern of EVI	Signature samples for supervised classification
LU-1	Mixed garden	72,059
LU-2	Upland mixed bush	28,392
LU-3	Upland with intensive agriculture farming	23,508
LU-4	Paddy rice with double-irrigated cropping system in wetland areas	10,869
LU-5	Paddy rice with double-rainfed cropping system, paddy–secondary crops pattern, where paddy is planted first followed by a second (different) crop after which the land remains fallow for the rest of the year (during the dry season). Most of this type is distributed around eastern Java.	12,160
LU-6	Paddy rice with double-irrigated rice cropping system in upland areas	23,015
LU-7	Paddy rice with triple irrigated cropping system in wetland areas	10,383
LU-8	Upland/dry land	36,174
LU-9	Fish/shrimp pond	15,269
LU-10	Upland with intensive upland crops with triple cropping system	3859
LU-11	Forest mixed with bush	10,235
LU-12	Bush–shrub mixed grass	13,487
LU-13	Dry-land forest (natural forest)	15,353
LU-14	Mangroves mixed (heterogeneous)	5,641
LU-15	Mangrove (in area with significant water fluctuation)	8,168
LU-16	Mixed garden with bush	12,121
LU-17	Rubber plantation	21,278
LU-18	Oil palm plantation	11,333
LU-19	Industrial forest plantation	45,860
LU-20	Upland with intensive agriculture farming	61,517
LU-21	Paddy rice in lowland areas	3,355
LU-22	Water body (lake, clear river)	9,387
LU-23	Settlement (built-up), mining area (open area)	19,479
LU-24	Paddy rice with double-rainfed upland-mixed cropping system in upland areas	3,891
LU-25	Mixed garden in upland area	5,032

Land use mapping and accuracy assessment

In order to understand the distribution of each land use type, the 25 type of patterns were used as signature samples to classify the entire study area (Table 3). The result of that classification is shown in Figure 9.

Figure 9. The distribution of land use type which was classified by the temporal pattern of enhanced vegetation index (EVI).

The use of temporal vegetation patterns allowed land use types to be differentiated into many specific land use classes, such as double and triple cropping system type of paddy rice, upland with intensive cropping and industrial forest plantation. Classification of specific classes was evaluated by 10,000 reference pixels, which revealed the overall accuracy and the overall κ statistics to be 71.07% and 64.94%, respectively (Table 4).

The greatest error was associated with the mixed garden, mixed bush and forest mixed bush classes. Specifically, 50.76% of mixed garden, 52.71% of bush–shrub mixed grass and 60.79% of forest mixed with bush were assigned to an incorrect category. The lower accuracy in these classes indicated the inability of MODIS to define the mixed type of vegetation and also reflected a difficulty with using temporal vegetation patterns to differentiate classes with no discernible periodical pattern. However, there is high confidence in pixels categorized as paddy rice and upland with several different cropping systems.

The classification results have been regrouped in terms of the 11 land use categories as follows: natural forest, mangrove, industrial forest plantation, plantation, bush, paddy rice field, mixed garden, upland, built-up/settlement, fishpond and water body. The overall accuracy and the overall κ statistics for the 11 land use categories are 80.16% and 70.94%, respectively. In addition, the omission error (producer's accuracy) and commission error (user's accuracy) are presented in Table 5.

The bulk of misclassification events contributing to the overall accuracy (19.84%) occurred mostly among the bush, plantation and built-up categories. Interpretation of bush is the most difficult, since this type did not exhibit a consistent temporal pattern, consequently almost a half of bush falls into other categories, such as mix garden and forest. The plantation types such as rubber plantation and oil palm plantation are also difficult to differentiate from forest based on their temporal vegetation profile, so almost 33% of the classified plantation falls into forest category in Table 5. In the case of the built-up category, the effective MMU is much smaller than each MODIS pixel, consequently each pixel contains a mixture of built-up and other classes, especially the agriculture lands, for example paddy rice, mixed garden and upland.

In the plantation and mixed garden cases, 52.11% and 46.80% of the error, respectively, is due to omission error (producer's accuracy), meaning that samples in these categories were assigned to other classes. The plantation type is determined by the area covered by a particular plantation crop (rubber and oil palm) which are widely distributed in Java Island. The poor accuracy of plantation and mixed garden types indicated that there is a difficulty with using MODIS data to resolve those areas.

Table 4. Confusion matrix for 25 specific classes of land use

	Reference data
Classified data	LU-1	LU-2	LU-3	LU-4	LU-5	LU-6	LU-7	LU-8	LU-9	LU-10	LU-11	LU-12	LU-13	LU-14	LU-15	LU-16	LU-17	LU-18	LU-19	LU-20	LU-21	LU-22	LU-23	LU-24	LU-25	Row total	Percent correct	Percent commission
LU-1	286	0	0	0	0	5	0	0	0	0	0	0	81	0	0	64	3	13	0	95	0	0	1	0	21	569	50.33	49.67
LU-2	75	680	0	0	0	0	0	0	0	0	0	0	6	6	1	0	1	0	0	0	0	0	61	0	3	834	81.50	18.50
LU-3	0	0	297	9	0	11	0	0	8	0	0	5	0	0	0	0	0	0	0	0	24	0	15	5	0	375	79.30	20.70
LU-4	0	0	11	278	0	1	8	0	14	0	0	0	0	0	1	0	0	0	0	5	16	0	19	0	9	361	76.92	23.08
LU-5	0	0	8	4	300	1	5	0	0	0	0	0	0	0	0	0	0	0	0	8	7	0	16	1	0	349	85.84	14.16
LU-6	3	0	0	8	0	342	3	0	3	9	61	0	6	1	0	0	1	0	0	0	0	0	13	0	0	449	76.22	23.78
LU-7	0	0	7	0	2	0	361	0	13	0	0	0	0	0	0	0	0	0	0	0	4	4	39	3	0	432	83.55	16.45
LU-8	0	0	0	0	0	2	2	144	0	0	0	0	2	0	0	0	1	2	15	23	0	0	6	0	22	219	65.73	34.27
LU-9	0	0	0	0	0	0	0	0	184	0	0	0	0	0	0	0	0	0	0	0	0	13	4	0	0	201	91.54	8.46
LU-10	0	0	0	0	0	1	2	0	0	332	6	0	15	0	2	0	3	0	0	18	0	0	2	1	0	382	86.93	13.07
LU-11	0	0	0	0	0	0	0	24	0	0	182	91	54	0	4	0	0	0	32	79	0	0	1	0	0	467	39.01	60.99
LU-12	0	31	0	0	0	3	0	14	1	0	9	138	0	0	0	0	8	0	0	0	0	0	5	4	0	213	64.81	35.19
LU-13	0	0	0	0	0	11	0	0	3	0	0	58	498	0	0	0	2	2	0	27	0	0	2	0	0	603	82.55	17.45
LU-14	0	0	0	0	0	0	1	4	1	0	29	0	47	181	0	0	1	0	0	0	0	2	29	0	6	301	60.23	39.77
LU-15	0	0	0	0	0	0	0	0	0	0	23	0	0	0	98	0	0	0	0	0	0	0	0	0	0	121	80.63	19.37
LU-16	0	0	0	0	0	0	0	0	0	0	5	0	6	0	0	266	0	0	0	22	0	0	4	0	0	303	87.68	12.32
LU-17	0	0	0	0	0	0	0	0	8	0	0	0	19	24	2	0	303	3	0	0	0	0	43	0	0	403	75.26	24.74
LU-18	0	3	0	0	0	7	2	0	0	1	87	0	0	0	0	0	12	468	0	19	0	0	6	0	0	604	77.43	22.57
LU-19	0	0	0	0	0	4	0	0	3	0	10	0	8	0	16	0	0	0	140	0	0	0	3	0	0	184	76.23	23.77
LU-20	217	0	0	0	0	0	2	0	0	0	0	0	12	0	0	0	4	11	0	601	0	0	12	0	0	859	69.94	30.06
LU-21	0	0	12	2	4	5	2	0	5	0	0	0	0	0	0	25	0	0	0	0	363	0	0	4	0	422	86.08	13.92
LU-22	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	114	0	0	0	115	99.17	0.83
LU-23	0	11	21	0	11	11	4	9	2	0	0	0	0	0	0	0	0	0	0	9	0	0	561	5	0	645	86.98	13.02
LU-24	0	0	0	1	12	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	2	4	4	299	0	325	92.00	8.00
LU-25	0	0	0	0	0	5	0	0	0	0	52	0	7	0	0	34	4	0	0	0	0	0	15	0	148	265	55.68	44.32
Column total	581	725	356	302	328	410	390	195	248	341	464	292	763	212	124	389	343	499	187	907	416	137	861	321	208
Percent correct	49.24	93.81	83.44	92.17	91.24	83.42	92.52	73.94	74.19	97.22	39.22	47.29	65.22	85.27	78.67	68.35	88.34	93.79	74.87	66.29	87.34	83.21	65.16	93.04	70.86
Percent omission	50.76	6.19	16.56	7.83	8.76	16.58	7.48	26.06	25.81	2.78	60.78	52.71	34.78	14.73	21.33	31.65	11.66	6.21	25.13	33.71	12.66	16.79	34.84	6.96	29.14		n = 10,000
				Overall accuracy = 71.07%
				Overall κ statistics = 64.94%

Table 5. Confusion matrix for 11 categories of land use

	Reference data
Classified data	Forest	Bush	Mangrove	Industrial forest	Plantation	Paddy rice	Mix garden	Upland	Built-up	Fishpond	Water-body	Row total	Percent correct	Percent commission
Forest	1680	15	39	0	19	35	20	0	0	0	0	1808	92.92	7.08
Bush	96	894	27	17	95	206	153	255	3	0	5	1751	51.06	48.94
Mangrove	0	0	124	0	0	0	0	0	0	0	0	124	100.00	0.00
Industrial forest	39	73	0	1365	60	81	74	52	5	0	0	1749	78.04	21.96
Plantation	256	1	0	6	499	3	0	3	0	0	0	768	64.97	35.03
Paddy rice	9	0	0	1	126	5091	97	38	48	0	2	5412	94.07	5.93
Mix garden	85	1	0	0	6	35	807	43	6	0	1	984	82.01	17.99
Upland	44	39	0	0	218	202	111	2554	16	0	4	3188	80.11	19.89
Built-up	0	41	0	0	19	293	255	131	1379	3	27	2148	64.20	35.80
Fishpond	0	0	0	0	0	132	0	0	0	341	0	473	72.09	27.91
Water body	0	0	0	0	0	0	0	0	0	0	100	100	100.00	0.00
Column total	2209	1064	190	1389	1042	6078	1517	3076	1457	344	139
Percent correct	76.05	84.02	65.26	98.27	47.89	83.76	53.20	83.03	94.65	99.13	71.94
Percent omission	23.95	15.98	34.74	1.73	52.11	16.24	46.80	16.97	5.35	0.87	28.06	n = 18,505
	Overall accuracy = 80.16%
	Overall κ statistics = 70.94%

Comparison of classification result of this study and land cover map derived by DoF

The classification result in terms of the 11 land use categories allows visual comparison with the land cover map provided by the DoF. These data sets were then overlaid to determine locations where the different land use type was present in both data sets.

Figure 10 shows that for certain types of land use, namely paddy rice, natural forest and uplands/drylands, both land use classification results appeared similar. In the case of the built-up class, the result of this study is slightly lower than the area extent indicated in the previous land cover map because of the inability of MODIS data to resolve the built-up areas. Meanwhile, the area extent of bush/shrub and mixed garden in the classification result is slightly greater than in the DoF map, because of many bare lands included in both classes. That probably occurred because the classification procedure used in this study was based on available imagery, which was not necessarily coincident with the presence of bare land.

Figure 10. Comparison of the 11 land cover/use categories derived by Department of Forestry and the classification result using of temporal pattern of Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI).

The most significant difference between our results and the previous land cover map is for the plantation and industrial forest plantation categories, as the areas of the categories seem to be reversed. This result indicates a difficulty with distinguishing between plantation and industrial forest plantation based on their temporal vegetation profile. Table 5 shows that almost 35% and 22% of the classified plantation and industrial forest plantation, respectively, fall into other categories.

The comparison of the land use map derived by classifying the temporal patterns and the existing land cover map is shown in the Figure 11. Significantly, the land use classification assigned by long-term fluctuation in vegetative activity results a simpler land use map than the land cover map based on a single-date image due to short-term fluctuations in vegetative activity.

Figure 11. (a) The land cover map by Department of Forestry; (b) the land use map produced using the temporal pattern of Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI).

Assessment of classification result for paddy rice

For the specific class of paddy rice, there is the possibility of comparing the temporal classification results with statistical data collected at the village scale (BPS-Statistics Indonesia 2003) through information about the irrigated paddy area. Paddy rice field can be planted 2 to 3 times per year in irrigated land; however, there are many aspects that impact the actual pattern, such as social capital, farmer welfare, irrigated water, the price of rice and climate.

Figure 12a shows comparison of the irrigated paddy rice area at the province scale between the classification result of MODIS EVI and the tabulations by BPS-Statistics Indonesia (2003). The areas estimated by MODIS EVI classification for western, central and eastern Java were 10% higher, 26% lower and 42% lower, respectively, than the areas determined from the statistical data. Moreover, Figure 12b indicated that the MODIS EVI procedure overestimates the area represented by large paddy areas and underestimates the number of small paddy fields. Analysis of the distribution revealed that most of the overestimates (points distributed below the best fit straight line) are located in western part of Java; most of the underestimates (points above the correlation line) are located in central and eastern Java. Such results suggest different characteristics of paddy rice on western Java as compared with central and eastern Java.

Figure 12. Comparison of the irrigated paddy rice area between the classification result of temporal pattern of Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI) and BPS-Statistics of Indonesia (2003). (a) comparison at the province scale, (b) scatterplot at sub-district scale.

The result above is related to the different size MMU for irrigated paddies in each region. In some areas of central and eastern Java, it is difficult to resolve the minimum unit of paddy area with MODIS EVI; consequently, it is also difficult to characterize the temporal dynamics of paddy rice. Even if the irrigation infrastructure exists locally in these areas, if irrigated water is relatively limited, particularly likely in east Java, double or triple cropping may not be possible in a given year.

Although there is a difficulty caused by the MMU of paddy fields, the characterization of paddy fields in this study is based on available imagery, which was not necessarily coincident with temporarily bare land when employing a single-date image. Temporal vegetation dynamics would provide sufficient information of the use of land; meanwhile the single-date image interpretation is not sufficient to identify the actual land use type because of insufficient documentation of change events.

Conclusion

Monitoring land surfaces continuously in space and time allows characterization of the temporal vegetation dynamics. Considering seasonal vegetation dynamics in a multi-year series leads to a broader view of LULC. In some cases, the land cover may change without any changes in land use type. The mixed use of land cover and land use in maps leads to inconsistencies in land change assessment. This study is based on the hypothesis that land use has a typical, distinct and repeated temporal pattern of the vegetation index inter-annually.

The temporal patterns of land cover, repeated over several years, was sufficient to distinguish a number of significant temporal signatures which were characteristic of different land use types. On Java Island, the temporal pattern of EVI profiles was differentiated into 25 types of specific land use and was then used for land use classification.

Accuracy assessment on 11 types of land use categories (category class level) and 25 types of specific land use classes (specific class level) showed the need to evaluate the effectiveness of the method when applied to land use for which there was no periodic pattern such as bush, mixed garden and forest. However, identification of the intensive agriculture lands (such as paddy rice and upland) in Java Island showed satisfactory results.

Land use classification using the MODIS EVI 250 m data is substantially limited by the moderate spatial data resolution (250 m) (Lobell and Asner 2004); land use types smaller than approximately 6.25 ha have a low probability of being identified. As indicated by the accuracy assessment, the resolution limits classification accuracy; smaller features (e.g. bush and built-up) are mixed with other classes. Although the pixel mixing issue is a significant problem with the MODIS data, this study shows that it still offers great promise for characterizing temporal dynamic vegetation at a large scale, since the methodology presented in this research showed successfully the identification of several specific land use classes, such as paddy rice and upland, based on the characteristics of temporal vegetation patterns.

Characterizing of temporal vegetation dynamics patterns would provide sufficient, significant and useful information regarding the patterns of land use; consequently the set of analyses and techniques in this study should be applicable to a large range of applications such as quantifying dynamics of carbon stock, water resource management as well as land use/agricultural planning in regional scale. The land use map developed in this article will be available for potential users by contacting the authors.

Acknowledgements

The authors would like to thank the NASA LP DAAC for making MODIS data available; Brett Lien and Janice Wilson (NASA LP DAAC) for help when downloading the raw data, the GLOVIS USGS and Global Land Cover Facility, University of Maryland for providing us with Landsat image data as supporting image, also the Center for Environmental Research (PPLH IPB) for providing the image and secondary data in this research. We would like to thank the anonymous reviewers for the helpful comments on the manuscript. The field survey in this research was financially supported by JSPS with the Grant-in-Aid for Scientific Research (B) No. 22405033 and the JST Basic Research Programs, Core Research for Evolutional Science and Technology (CREST) in Japan.