Can a land cover map derived through hyper-temporal NDVI images be improved using NDWI information

(1)

CAN A LAND COVER MAP DERIVED THROUGH HYPER- TEMPORAL NDVI IMAGES BE IMPROVED USING NDWI

INFORMATION?

DAMAYANTI SARODJA March, 2011

SUPERVISORS:

Dr. Ir. C.A.J.M. (Kees) de Bie Dr. Ir. Anton Vrieling

(2)

Thesis submitted to the Faculty of Geo-Information Science and Earth Observation of the University of Twente in partial fulfilment of the

requirements for the degree of Master of Science in Geo-information Science and Earth Observation.

Specialization: Natural Resources Management

SUPERVISORS:

Dr. Ir. C.A.J.M. (Kees) de Bie Dr. Ir. Anton Vrieling

THESIS ASSESSMENT BOARD:

Dr. Y. Hussin (Chair)

Dr. A.J.W. de Wit (External Examiner, Alterra, Wageningen University)

CAN A LAND COVER MAP DERIVED THROUGH HYPER-TEMPORAL NDVI IMAGES BE IMPROVED USING

NDWI INFORMATION?

DAMAYANTI SARODJA

Enschede, The Netherlands, March, 2011

(3)

DISCLAIMER

This document describes work undertaken as part of a programme of study at the Faculty of Geo-Information Science and Earth Observation of the University of Twente. All views and opinions expressed therein remain the sole responsibility of the author, and do not necessarily represent those of the Faculty.

(4)

Land cover can be successfully mapped through-the use of hyper-temporal Normalized Difference Vegetation Index (NDVI) imagery. Normalized Difference Water Index (NDWI) is a different indicator, not correlated to NDVI, which responds to differences in vegetation water content. As different land covers have different temporal behaviour of vegetation water content, theoretically NDWI can provide input to land cover mapping.

This research investigates if additional features of land cover can be mapped after the main variability in land cover is captured through the use of hyper-temporal NDVI imagery. The main data used in this research are the 10-day Maximum Value Composite of SPOT VEGETATION NDVI and NDWI acquired from April 1998 to May 2010. The study area is in Andalucia, Southern Spain.

Three NDVI units were selected for study in this research. They represent one area dominated by agriculture (NDVI unit 19), one area dominated by forest and semi natural vegetation mixed with transitional woodland shrub and herbaceous (NDVI unit 45), and one area dominated by forest and semi natural vegetation (NDVI unit 57). The information in the NDWI images are reduced through an unsupervised classification procedure called Iterative Self-Organizing Data Analysis (ISODATA) in ERDAS Imagine software. The differences in NDWI profiles between classes proved consistent across the years. These differences can only be caused by true differences between the ecosystems they represent.

The process followed by extracting four parameters that characterized the shape of the NDWI-profiles by class. The four derived NDWI parameters are maximum, minimum, length of cycle and amplitude.

The land cover related data collected were percentage cover of tree, shrub, herb and grass, litter, stone and soil. Using orthophotos (2004), field data were converted into one square kilometre pixel information. In turn, that information was correlated with the four NDWI parameters.

NDWI parameter explains 32% variability of tree characteristics in NDVI unit 19, 49% shrubs characteristics in NDVI unit 45, and 25% shrubs characteristics in NDVI unit 57 at 99% confidence level.

It is evidence that among the four parameters, amplitude is the parameter that give a better relationship with cover characteristics.

The results of this research show that linear regression models had overall low goodness-of-fit between NDWI parameters and cover characteristics. It indicates that the chosen cover related data do not represent the anticipated cover related variability that would explain the differences in derived-NDWI profiles.

(5)

ii

ACKNOWLEDGEMENTS

My sincere thank and appreciation to several individuals who have supported me in completing this thesis.

I would like to express my gratitude to:

The European Union, for giving me the opportunity and the fellowship to pursue this MSc.

Dr. (Kees) de Bie, my first supervisor, for the guidance, criticism, especially the patience during the research phases. I would not complete all these phases work without his help.

Dr. Anton Vrieling, my second supervisor, for the valuable comments and suggestions during my thesis writing.

Dr. Y. Hussin (Chairman) and Dr. A.J.W. de Wit (External Examiner), for the willingness to examine my work.

Dr. Michael Weir, NRM program director, for the kind efforts, considerations and supports.

Drs. Tom Loran, for the kind help and support all the way long.

All my Indonesian family in Enschede, Syarif and Astrid, Elham, Andi, Aniq…. all my NRM classmates and ITC friends for the very nice friendship and memorable times.

Twinny, Sally and Sarita for all the help, nice friendships, and much more; Foxtorn for the support, help, and share the ideas; Ha, Amjad, and Riaz for sharing their expertise.

Finally, my heart-felt gratitude and thanks for my beloved family; my husband Bobby Berlianto, my parents, my brothers and sisters, also my nephews and nieces, for their extended prayer, endless love, and encouragement. There is no word can express how grateful I am to be part of you.

(6)

1. INTRODUCTION... 1

1.1. Background ... 1

1.1.1. Retrieving information from remotely sensed data for vegetation monitoring ... 1

1.1.2. Hyper-temporal vegetation indices for land cover classification ... 2

1.1.3. NDWI, a complementary index to NDVI land cover classification ... 3

1.2. Problem Statements (Research gap) and Research Objective ... 4

1.2.1. Problem statements ... 4

1.2.2. Research objective ... 5

1.3. Research questions ... 5

1.4. Research hypothesis ... 5

2. MATERIALS AND METHOD ... 6

2.1. Materials ... 6

2.1.1. Description of the study area ... 6

2.1.2. Description of Data used ... 7

2.2. Method ... 8

2.2.1. Data preparation phase ... 10

2.2.2. Fieldwork phase ... 13

2.2.3. Data Processing phase ... 14

2.2.4. Analysis phase ... 16

3. RESULTS AND DISCUSSIONS ... 18

3.1. Data preparation phase ... 18

3.1.1. Generating land cover maps derived through hypertemporal NDVI and NDWI imagery ... 18

3.1.2. Generating derived NDWI parameter maps ... 20

3.2. Fieldwork phase ... 21

3.2.1. Selecting area for fieldwork ... 21

3.2.2. Collecting field data ... 22

3.3. Data processing phase ... 22

3.3.1. Creating image legend and averaging per cent cover ... 22

3.3.2. Upscaling field measurement to a 1 km² pixel size ... 24

3.3.3. Calculating per cent cover type characteristics per pixel ... 24

3.4. Analysis phase ... 24

3.4.1. Regression analysis ... 25

4. CONCLUSIONS AND RECOMMENDATIONS ... 31

4.1. Conclusions ... 31

4.2. Recommendations ... 33

(7)

iv

LIST OF FIGURES

Figure 1.1. The result of twelve-years hyper-temporal derived-NDVI land cover map ... 2

Figure 1.2. Legend of twelve-years hyper-temporal derived-NDVI land cover map ... 3

Figure 1.3. Considered heterogeneous NDWI within a hyper-temporal derived-NDVI land cover map unit (A) Heterogeneous NDWI profiles within a NDVI unit (B) ... 4

Figure 2.1. Research study area ... 6

Figure 2.2. Flowchart of proposed procedures ... 9

Figure 2.3. A profile of NDVI before and after applying a modified Adaptive-Savitzky-Golay filter; taken at a pixel with coordinate (-545664, 4114845), Plate Carree / WGS 84 projection system ... 10

Figure 2.4. A profile of NDWI before and after applying a modified Adaptive-Savitzky-Golay filter, taken at a pixel with coordinate (-545664, 4114845), Plate Carree / WGS 84 projection system ... 10

Figure 2.5. The NDWI derived parameter ... 12

Figure 2.6. Converting NDWI profiles into four parameters ... 12

Figure 2.7. Samples collected within NDVI unit 19 using image objects of orthophotos ... 14

Figure 2.8. An image legend with three groups of image objects ... 15

Figure 2.9. Upscaling approach from field measurement into 1 km² pixel size; the grid frame is using local projection (UTM, WGS 84, ED 50 – 13) ... 16

Figure 3.1. Separability analysis for NDVI; the left Y axis shows the average separability (in red colour) and the right Y axis shows the minimum separability (in blue colour). ... 18

Figure 3.2. Separability analysis for NDWI; the left Y axis shows the average separability (in red colour) and the right Y axis shows the minimum separability (in blue colour). ... 18

Figure 3.3. Land cover map derived hyper-temporal NDVI 60 classes images classification... 19

Figure 3.4. Land cover map derived hyper-temporal NDWI 45 classes images classification ... 19

Figure 3.5. Maps of derived NDWI parameters ... 20

Figure 3.6. Distribution of selected NDVI map unit for fieldwork ... 21

Figure 3.7. Profile of NDVI unit 19 and the included NDWI Profiles ... 24

Figure 3.10. Scatterplot of Litter (%) and Amplitude parameter in NDVI unit 19 (left) and in NDVI unit 57 (right) ... 27

(8)

Table 2.1. Spectral characteristics of SPOT-VGT 1 and 2 ... 7

Table 2.2. Average percent cover per each group per (layer) cover type ... 15

Table 3.1. Standard deviation (Stdev) of the 14 NDVI polygons selected for fieldwork ... 21

Table 3.2. Summary table on field data collected ... 22

Table 3.3. Summary of average cover for NDVI unit 19 per group (%) ... 22

Table 3.6. Count of derived NDWI parameter values and number of pixels digitized per NDWI parameter values within each NDVI unit ... 24

Table 3.7. R-square of individual cover characteristics versus individual NDWI derived parameters ... 26

Table 3.8. R-square of Eucalyptus and Needle leaves characteristics versus individual NDWI derived parameters ... 28

Table 3.9. R-square of cover type characteristics versus NDWI derived parameter from original pixel value ... 29

(9)

vi

ANNEXES

Annex 1. Relevee sheet for field data collection

Annex 2. The legends of landcover maps derived hyper-temporal NDVI (60 classes) and NDWI (45 classes) imagery

Annex 3. Derived hyper-temporal NDWI parameters

Annex 4. The standard deviation of derived NDWI parameters within selected NDVI unit Annex 5. Image legend based on field data collected for NDVI unit 19

Annex 6. Image legend based on field data collected for NDVI unit 45 Annex 7. Image legend based on field data collected for NDVI unit 57 Annex 8. Fraction area of the groups per pixel for NDVI unit 19 Annex 9. Fraction area of the groups per pixel for NDVI unit 45 Annex 10. Fraction area of the groups per pixel for NDVI unit 57 Annex 11. Fraction cover per type per pixel for NDVI unit 19 Annex 12. Fraction cover per type per pixel for NDVI unit 45 Annex 13. Fraction cover per type per pixel for NDVI unit 57

Annex 14. R-square of joint cover characteristics versus individual NDWI derived parameters

Annex 15. R-square of joint group and cover characteristics versus individual NDWI derived parameters

(10)

1. INTRODUCTION

1.1. Background

1.1.1. Retrieving information from remotely sensed data for vegetation monitoring

Remote sensing technology are widely applied to understand the physical and biological characteristics of the land surface (Maselli, et al., 1998). Remotely sensed data are used to retrieve information on vegetation parameter through the spectral indices (Ceccato, et al., 2001; Ceccato, et al., 2002; Delbart, et al., 2005). Spectral indices are “mathematical transformations designed to assess the spectral contribution of green plants to multispectral observations”(Maselli, et al., 1998). This approach is advantageous to apply in a large-scale area because of its simplicity and efficiency computation (Yi, et al., 2007).

NDVI is a vegetation index developed by Rouse et al. (1973) and widely used for remote sensing of vegetation (Gao, 1996; Reed, et al., 1994). This index is derived from reflectance values of the red and the NIR band using the formula (Tucker, 1979) ;

NDVI NIR — Red NIR Red

where,

Red: wavelength located in the strong chlorophyll absorption region,

NIR: wavelength located in the region of high reflectance of vegetation canopy.

The algorithm normalizes the difference between the spectral reflectance in NIR and red bands. Derived- NDVI values are related to ground cover percentage, total green biomass, leaf area index (LAI), and photosynthetic activity of the vegetation (Justice et al., 1985). NDVI values range between -1 to 1. Higher NDVI values indicate higher level of greenness of the vegetation.

However, NDVI is known to have some limitations. NDVI values tend to get saturated, or become no longer responsive to crop growth, in a multilayer closed canopy with a large leaf area index (LAI) values (Gao, 1996). This phenomenon is most likely due to the strong chlorophyll absorption by the red wavelength (Ceccato et al., 2002) also the fact that red wavelength could only detect one leaf layer or less (Gao, 1996). NDVI is also sensitive to atmospheric conditions and soil background (Gao, 1996; Huete, et al., 2002). Other spectral regions have been examined to answer this limitation.

Tucker (1980) first proposed that Short Wave Infra Red wavelength (SWIR) is sensitive to vegetation canopy water content. Region between ~1.3 – 2.5 μm spectral interval shows the quantity of liquid water present in the leaf biomass of the plant canopy (Tucker, 1980). SWIR absorption increases as leaf liquid water content or soil moisture increases, resulting in the decrease of SWIR reflectance.

Gao (1996) introduced a new vegetation index named Normalized Difference Water Index (NDWI) to assess the vegetation liquid water content. This index is using the same algorithm as NDVI; however NDWI combines the NIR and SWIR spectral bands. The formula is,

(11)

CAN A LAND COVER MAP DERIVED THROUGH HYPER-TEMPORAL NDVI IMAGES BE IMPROVED USING NDWI INFORMATION?

2

NDWI NIR — SWIR

NIR SWIR where

SWIR wavelength located in the region of vegetation canopy reflectance with water absorption

The SWIR band is sensitive to soil moisture, and leaf moisture content which can improve the discrimination of vegetation and other land covers (Xiao, et al., 2002). Since it is influenced by the desiccation and wilting in the vegetation canopy, NDWI can be used to detect and monitor the moisture condition of vegetation canopies (Xiao et al., 2002; Jackson et al. 2004) and thus give a better prediction for the start of the growing season (Boles, et al., 2004). Unlike NDVI, NDWI becomes saturated later than NDVI in a closed canopy cover and high LAI, and is less sensitive to atmospheric conditions (Gao, 1996; Huete, et al., 2002).

1.1.2. Hyper-temporal vegetation indices for land cover classification

Information on land cover variability is very useful in natural resources management (Boles, et al., 2004).

However, land cover information is time sensitive and can change dramatically over seasons (Agrawal, 2003). Thus, the repeated vegetation capture can improve the land cover information, especially in an area where the land cover or land use rapidly changes (de Bie, 2008).

Figure 1. shows a land cover map derived from twelve-years hyper-temporal SPOT-VGT NDVI data.

The map has a legend of NDVI-profiles (Figure 1.2) which correspond to the seasonal characteristics of vegetation in the area over twelve years. The figure shows the difference between NDVI-profiles caused by a combination of land cover type and its ecosystem such as weather, soil, and terrain (de Bie, 2008).

Figure 1.1. The result of twelve-years hyper-temporal derived-NDVI land cover map

(12)

Figure 1.2. Legend of twelve-years hyper-temporal derived-NDVI land cover map

The interpretation of land cover derived from hyper-temporal NDVI classification is based on the variability of the profile patterns (Figure 1.2). Reed, et al. (1994) defined 12 seasonal parameters to interpret the phenological variability, adapted from Malingreau’s (1986) and Lloyd’s (1990) research. The 12 parameters are divided into 3 types: temporal (based on the timing of an event), NDVI-based (the value of NDVI when the events occurred), and derived metrics (based on time-series characteristics).

The temporal NDVI metric is related to the time of onset, end, duration, and maximum of greenness, which is related to the time of measurable photosynthetic activity. The derived metrics are related to the time-integrated NDVI (NPP), rate of greenup and senescence (acceleration and deceleration of photosynthesis), and to modality (periodicity of photosynthetic activity).

Time series analysis of NDVI has been widely used for land cover and vegetation classification as well as vegetation change monitoring (Ferreira & Huete, 2004). It reflects the phenological patterns of vegetation. Using time series NDVI data is an advantage (Defries & Townshend, 1994) and provides information on the vegetation types (Knight, 2006). Another advantage of using a time series vegetation index is that atmospheric effects such as the effects of clouds are reduced (Doraiswamy, et al., 2007).

However, the use of high temporal dimension compensates for the low spectral and spatial dimension (Udelhoven, et al., 2009).

Multi temporal, time series or hyper-temporal data have the same meaning in the sense that they deal with data from many different observation (capture) times. The term “time series” means that the repeated time of observation is periodical. The term time series and hyper-temporal essentially have the same meaning, however the use of hyper-temporal indicates that the number of repeated time of observation used is huge. “The term ‘hyper’ in remote sensing refers to a condition where the question of interest is fully determined in spectral resolution (hyperspectral), in the spatial dimension (hyperspatial), or in time (hypertemporal)” (Chambers et al., 2007).

1.1.3. NDWI, a complementary index to NDVI land cover classification

Recently, NDWI has been used together with NDVI as input to efforts for land cover mapping ( Xiao et al., 2002; Bartalev et al., 2003; Lu et al., 2003). The combination of NDVI and NDWI can be used to

0 50 100 150 200 250

1 37 73 109 145 181 217 253 289 325 361 397 433

class46 class47 class48 class49

Time (decade)

NDVI (0‐255)

(13)

4

generate a relevant land use thematic map and also can obviously provide improved description of the growth status of vegetation canopies (deAlwis et al., 2007).

Xiao, et al. (2002) revealed that VGT-derived NDVI time series can be used to separate dominated evergreen tree species forest from the one dominated by deciduous species in Northeastern China. The inclusion of VGT-derived NDWI time series added unique information on forest type and dynamics of evergreen needle leaves, deciduous needle leaves and deciduous broadleaves. Furthermore, Boles, et al.

(2004) stated that the use of NDWI time series has the potential to define the vegetation growing season.

Jackson, et al. (2004) proved that NDWI is superior to NDVI in mapping the daily vegetation water content for corn and soybean. Yi, et al. (2007) also tested the use of NDWI for winter wheat, and concluded that NDWI can be used as an indicator of crop water deficiency.

However, NDVI and NDWI are independent, and each contains different vegetation information. Gao (1996) stated that one cannot replace the other. NDVI shows the chlorophyll content of vegetation while NDWI shows the vegetation water content.

1.2. Problem Statements (Research gap) and Research Objective 1.2.1. Problem statements

The accuracy of land cover maps relies on the input information available (Lu et al., 2003). Hyper- temporal NDVI images classification developed by ITC through unsupervised classification technique of ISODATA and separability analysis using ERDAS Imagine software is a sufficiently proven method for land cover mapping capturing dynamic land cover/vegetation aspects (de Bie et al., 2008; Khan et al., 2010), although the method may be further improved.

There is a need to assess the potential additional use of other indices to produce an improved-land cover classification method. Previous research proved that NDVI and NDWI can be combined to give additional information (Gao, 1996; Xiao et al., 2002). However, study on how the relationship between NDWI and cover characteristics using field measurements and statistical analysis is still limited.

A preliminary investigation was made to define the objective of this research. The result showed that after SPOT VGT-derived NDVI and NDWI hypertemporal classification, variability of NDWI classes exists within one NDVI-derived land cover unit as shown in Figure 1.3.

Figure 1.3. Considered heterogeneous NDWI within a hyper-temporal derived-NDVI land cover map unit (A) Heterogeneous NDWI profiles within a NDVI unit (B)

A B

(14)

1.2.2. Research objective

The objective of this research is to investigate whether NDWI can be used to explain the remaining land cover variability after a NDVI hyper-temporal images classification process is carried out.

1.3. Research questions

1. Can the remaining variability within a land cover map unit as derived from NDVI hyper- temporal classification be explained by a NDWI derived parameter?

2. Which of the NDWI derived parameters can better explain the variability of land cover characteristic?

3. Does combination of four NDWI derived parameters give a significant difference in explaining the remaining variability within a land cover map unit, compared to a single parameter as derived fromNDVI hyper-temporal classification

1.4. Research hypothesis

Related to research question 1 and 3, the following research hypothesis will be tested,

Ho: within a NDVI-derived land cover map unit, a NDWI derived parameter cannot explain for at least 60% of the remaining variability of a specific land cover characteristic at a 95% confidence level

Ho: within a NDVI-derived land cover map unit, all four NDWI derived parameters cannot explain for at least 60% of the remaining variability of a specific land cover characteristic at a 95% confidence level

Since previous research have not mentioned per cent cover characteristics explained by the NDWI, author selected 60% to be tested in this research. 60% is selected otherwise the effort provides too little contribution to any improvement made to the land cover map (author’s opinion).

(15)

6

2. MATERIALS AND METHOD

2.1. Materials

2.1.1. Description of the study area

The study area of this research covers almost all the administrative area of Andalucia, Southern Spain.

Andalucia is an autonomous community of Spain and is divided into eight provinces. It is located between latitudes 36º00' and 38º72'N and longitude 7º51' and 1º63'W, covering an area of 87,597 square kilometres. The northern part of Andalucia, is not part of the study area due to the unavailability of free NDWI data from the SPOT-VGT sensor (Figure 2.1)

Figure 2.1. Research study area

The climate of Andalucia is characterised by mild, with rainy winters and hot dry summers (Khan, et al., 2010). Two municipalities have the highest average temperature in Spain. These are Almeria with 18.6 ºC and Seville with 18.7ºC. The hottest months are July and August, while January is the coldest. Annual precipitation in the study area is very heterogeneous following a decreasing gradient from west to east, ranging from a maximum of 2000 mm to a minimum of 170 mm (Font, 2000). Autumn and winter are the wettest seasons, with heavy rainfall typically occurring during September and October. The gradual change of climatic conditions, from west to east, was one of the factors taken into consideration when selecting the study area.

(16)

2.1.2. Description of Data used

SPOT Vegetation NDVI and NDWI

Several factors were taken into consideration when deciding to use the SPOT VEGETATION (VGT) sensor for the NDVI and NDWI data. The data provided was geometrically, atmospherically and radiometrically corrected, and also projected in a standard projection. In addition, the NDVI data are freely available to download as the S10 SPOT VGT data product through the website www.vgt.vito.be, and NDWI data are available at www.vgt4africa.org. The S10 is a 10-day VEGETATION synthesis product. It is a composite of daily images (of NDVI and NDWI) over a 10-day period generated using the Maximum Value Composites (MVC) algorithm. The 10-day period referred are the 1st-10th, 11th-20th, and 21st to the last day of the month. The algorithm selects the pixels with best ground reflectance values during the 10-day period (Vegetation programme webpage).

The data are available continuously starting from April 1998 to recently; the latest freely available data are three months old. Although there are other datasets existing (e.g. MODIS products), NDVI and NDWI of SPOT-VGT with 1 kilometre spatial resolution were considered sufficient to investigate the given objective.

The NDVI and NDWI datasets used in this study span from 1^st April 1998 to 31^st May 2010 (altogether 438 images). These are standard products derived from the SPOT VGT 1 instrument on board the SPOT 4 satellite and the SPOT VGT 2 instrument on board SPOT 5. The SPOT is a joint programme by France, the European Commission, Belgium, Italy and Sweden (Vegetation programme webpage). The spatial resolution of the data is 1 kilometre and the temporal resolution is 1 day. The spectral characteristics of the two sensors are shown in table 2.1.

Table 2.1. Spectral characteristics of SPOT-VGT 1 and 2 Spectral

bands

Specified (μm)

VGT‐1 (μm) (actual values)

VGT‐2 (μm) (actual values)

Surface

reflectance range Blue (B0) 0.430 ‐ 0.470 0.437 – 0.480 0.438 – 0.475 0.0 – 0.5

Red (B2) 0.610 – 0.680 0.615 – 0.700 0.615 – 0.690 0.0 – 0.5 NIR (B3) 0.780 – 0.890 0.772 – 0.892 0.782 – 0.890 0.0 – 0.7 SWIR (MIR) 1.580 – 1.750 1.600 – 1.692 1.582 – 1.685 0.0 – 0.6 (Source: http://www.spot-vegetation.com/vegetationprogramme/index.htm)

Although the MVC algorithm has been applied, the original dataset still shows some remaining noise such as clouds, ozone, dust or other aerosol (Goward, et al., 1991; Hird & McDermid, 2009) which lead to the decrease of NIR radiance. Further cleaning to remove remaining cloud effects and to improve the quality of the NDVI and NDWI hyper-temporal data is often required.

The NDVI and NDWI data of the SPOT VGT were already converted from its original values (-1 to 1) into unsigned-8-bit integers. The conversion of NDVI from its original values into 0-255 digital number (DN) is based on the formula (VG4Africa user manual, 2006):

DN _. ^. (equation 1)

The same formula is applied for the conversion of NDWI values; thus the physical NDVI and NDWI values are between -0.1 and 0.92 (VGT4Africa user manual, 2006).

(17)

8

Topographic maps of Andalucia

A topographic map of Andalucia at 1:100.000 scale in shapefile format was also used. The map was obtained from the Junta de Andalucia. The map was mainly used to support field data collection and map layout.

Orthophotos

Orthophotos of year 2004 were used since they were the latest data available when the research started in September. The scale of the orthophotos is 1:10.000, and available to downloading freely through www.juntadeandalucia.es.

2.2. Method

The research work is divided into 4 phases (Figure 2.2): data preparation, fieldwork, data processing, and data analysis phase. There are several main steps in each phase;

1. The data preparation phase concerns all the work by:

 Generating land cover maps derived through hyper-temporal NDVI & NDWI imagery (DP-1)

 Generating derived NDWI parameter maps (DP-2) 2. The fieldwork phase concerns to all the work by:

 Selecting area for fieldwork (FW-1)

 Collecting field data (FW-2)

3. The data processing phase concerns to all the work by:

 Creating image legend and averaging per cent cover (PR-1)

 Up scaling field measurement into 1 km² pixel size (PR-2)

 Calculating per cent cover type characteristics per pixel (PR-3) 4. Data analysis phase is mainly the regression analysis (AN)

Leica geosystems - ERDAS Imagine was used for image processing and ESRI - ArcGIS software for GIS spatial analysis; while R statistics software were chosen for the regression analysis.

(18)

Figure 2.2. Flowchart of proposed procedures

(19)

10

2.2.1. Data preparation phase

All steps in this phase relate to hyper-temporal NDVI and NDWI image classification. The NDVI and NDWI dataset were already cleaned from noise and classified into 90 predefined classes. The output consists of maps of classified hyper-temporal NDVI and NDWI, and their annual mean profiles.

DP-1. Generating land cover maps through hyper-temporal NDVI & NDWI imagery

 The NDVI and NDWI datasets were cleaned from noise using a modified Adaptive-Savitzky-Golay filter. This filter is a modification from the Adaptive-Savitzky-Golay that is originally based on Savitzky-Golay (Beltran, 2009). SAVGOL (Savitzky & Golay, 1964) is a filter used to approximate the series using a windowed polynomial least-square fit which preserves the higher order moments of the series (Jönsson & Eklundh, 2004). Figure 2.3 shows a hyper-temporal profile of NDVI from a pixel before and after the Savitzky-Golay filtering process, while Figure 2.4 shows a profile for NDWI.

Figure 2.3. A profile of NDVI before and after applying a modified Adaptive-Savitzky-Golay filter; taken at a pixel with coordinate (-545664, 4114845), Plate Carree / WGS 84 projection system

Figure 2.4. A profile of NDWI before and after applying a modified Adaptive-Savitzky-Golay filter, taken at a pixel with coordinate (-545664, 4114845), Plate Carree / WGS 84 projection system

The profiles represent a hyper-temporal NDVI and NDWI data taken from a single pixel at the same location. It shows that the application of the filter resulted in a smoother NDVI and NDWI profile.

It is known that NDVI and NDWI data are susceptible to cloud cover which lead to underestimated value. Furthermore, the SWIR band of SPOT VGT bring in an interpolated or missing lines or pixels due to the blind detectors (VG4Africa user manual, 2006). This explained the zero values on NDWI profiles as shows in Figure 2.4.

0 50 100 150 200 250

1 37 73 109 145 181 217 253 289 325 361 397 433

NDVI Before NDVI After

time (decade)

NDVI (0‐255)

0 50 100 150 200 250

1 37 73 109 145 181 217 253 289 325 361 397 433

NDWI Before NDWI After

NDWI (0‐255)

time (decade)

(20)

 An unsupervised classification method available in ERDAS-imagine software, called the Iterative Self-Organizing Data Analysis (ISODATA) technique, used to classify the 438 stacked images of NDVI and NDWI separately. The classification run to generate pre-defined classified maps from 10 to 100 classes. Each classification was set to 50 iterations and 1.0 convergence threshold. The ISODATA procedure is iterative; it repeatedly performs an entire classification and recalculates statistics (Khan, et al., 2010). It is self-organizing regarding the way in which it locates the clusters that are inherent in the data and minimizes the Euclidian distances to form clusters (Leica Geosystem, 2005; De Bie et al., 2008).

 To define the best number of NDVI and NDWI classes (between 10 and 100), the cluster separability results were examined. This separability was measured by the divergence between each pair of classes of the unsupervised classification result using the signature evaluator function in ERDAS (Leica Geosystem, 2005). The best number of classes to be chosen was based on the maximum values of average and minimum divergence statistics with a clear and distinct peak in the separability values of divergence. The minimum separability is considered less important, and may be ignored when it is lower than the value of 26, especially when gradients are part of the dataset (de Bie, 2010; personal communication). Furthermore, the lesser number of classes is chosen when more than one divergence combination is suitable (de Bie, 2010; personal communication).

Assumption: The chosen number of classes, by using distinct peak of average and minimum divergence values, is the best for unsupervised classification of NDVI and NDWI. The land cover map derived from 12 years hyper-temporal NDVI or NDWI imagery represents sufficiently the variability in land cover of the study area.

DP-2. Generating derived NDWI parameter maps

 The first step is to generate the mean annual profile of NDWI to get a yearly representative profile by class. Mean annual profiles of NDWI are defined by calculating the average values of each decade (1^st to 36^th decades) over 12 years for each class of NDWI.

Assumption: A class profile represents all variability in the respective class; mean annual profiles represent sufficiently the main variability of the 12 years profiles

 The second step is to convert NDWI profiles into parameters. The average annual NDWI profile by class was, converted into four parameter values. The logic behind this conversion method is to get some measured parameters to use as explanatory variables for the anticipated regression analysis. The interpretation was carried out at an NDWI class level, and not at the pixel level. This generalisation process was chosen in order to simplify the complex variability provided by the data.

The parameters are adapted from Reed et al. (1994), originally used for NDVI phenological interpretation. Reed et al. (1994) divided seasonal NDVI metrics into 3 categorical metrics, which are temporal, value and derived metrics. Parameters selected for this study are categorised as value and temporal metrics. They are considered to be a proper representation of the NDWI profiles when explaining the relationship with cover characteristics using regression analysis.

(21)

(22)

2.2.2. Fieldwork phase

Output of this phase is all data collected in the field. The fieldwork took place between 13th September and 2nd October 2010. IPAQ with ArcPad software was used as an instrument to collect geo-located data in the field.

The data were collected jointly with two other ITC students (1 PhD and 1 MSc) working in the same area, using the same NDVI dataset, but for different research topics. Since time was limited, several adjustments had to be taken when preparing the sampling scheme in order to satisfy all three topics.

The sampling scheme used for field data collection was clustered random (representative) sampling. The sampling points visited in the field were image objects as taken from 2004 orthophotos. All image objects within a pre-selected NDVI map unit were sampled.

FW-1. Selecting area for fieldwork

 Only NDVI map units with a polygon area of at least 50 kilometres square or larger were selected as potential units of observation. This was done to avoid selection of areas with only a few pixels in the unit, and to ensure that the investigation of NDWI variability was carried out in a sufficient large area.

Assumption: An area of at least 50 kilometres square or larger was considered sufficient to study the variability of NDWI parameter values

 Another criterion considered when selecting a sampling area is the variability of the derived NDWI parameters within the potential NDVI map units. Hence, the standard deviation by parameter within the potential NDVI map unit was calculated.

 NDVI map units with high standard deviations are potential sampling areas. Therefore, the sum of all standard deviation values from the derived NDWI parameters is calculated. To ensure that all parameter have an equal weight, the standard deviation values of each parameter were normalized to a 0 to 1 scale.

FW-2. Collecting field data

Orthophotos of year 2004 with 1:10.000 scale were used for field data collection in the selected NDVI units. A sample taken in the field represent an image object shown by the orthophotos. An image object has a unique characteristic in the orthophoto and represents an object in the field (e.g. olive field, pine forest). Land cover related data such as percentage cover of tree, shrub, herb and grass, stone, litter, and bare soil were collected. The percentage were estimated as a vertical cover, thus the total percentage for all cover in one sample point is 100. The datasheet used for field data collection is presented in Annex 1.

The nomenclature used for the fieldwork data is as follows:

Tree layer : All tree with height = > 2 meter

Shrub layer : All vegetation with height between = 50 centimeter to 2 meter Herb and grass layer : All vegetation with height < 50 centimeter

Stone layer : All rock with diameter = > 2 centimeter Litter layer : All dead matter on the ground

Bare soil : Soil in the ground

Assumption: There is no major change in land cover from 2004 to 2010

(23)

14

2.2.3. Data Processing phase

Output of this phase is a summary table of upscaled land cover type characteristics per one kilometre square pixel.

PR-1. Creating image legend and averaging per cent cover

Similar image objects were put into groups to create an image legend of an NDVI unit. The image legend is used to interpret the orthophoto.

 Using the orthophotos, snapshots of all image objects were prepared. A snapshot of an image object represents a sampled point of estimated cover characteristics in the field (Figure 2.7).

Figure 2.7. Samples collected within NDVI unit 19 using image objects of orthophotos

Image objects were put into groups based on their similarities in pattern (regular grid, dots, stripes etc), field (olive field, fallow field, etc) and texture (smooth, rough, irregular, etc). Groups of image objects (e.g. group A, B, C, etc) formed an image legend; this is created for a specific NDVI unit (Figure 2.7). Samples located at the edge between two image objects were excluded from the list.

Samples with descriptions that did not match the image object were also excluded (change did happen).

Figure 2.8 is an example of the image objects grouping. Group A is a field characterized by the regular dot patterns; group C is characterized by regular stripes and rough patterns; while group D is characterized by the irregular smooth to rough green patterns. Differences in soil colour were not considered.

(24)

Figure 2.8. An image legend with three groups of image objects

 An average of all cover characteristics estimated during fieldwork is calculated for each image object group.

Table 2.2. Average percent cover per each group per (layer) cover type

PR-2. Up scaling field measurement to 1 km² pixel size

The upscaling approach in this study was done through a digitizing process at a pixel level.

 The NDWI raster image is used to create grid pixels.

 Digitizing the image object was done in ArcGIS through visual interpretation at orthophotos pixel resolution. The image legend was used to assign labels to the created polygons. The digitizing was extended to pixels that were not visited in the field but of which all objects were present in the image legend (Figure 2.9).

(25)

16

Figure 2.9. Upscaling approach from field measurement into 1 km² pixel size; the grid frame is using local projection (UTM, WGS 84, ED 50 – 13)

Assumption: Within one NDVI unit, similar image objects shown on the orthophotos represent similar land cover characteristics as described by the image legend.

 After the digitizing process, the area of each polygon is calculated. The results of this step is a table showing the area of each group (e.g. group A, B, C, etc) for a specific pixel grid.

PR-3. Calculating land cover characteristics per pixel

Land cover characteristics per one kilometre square pixel size were calculated by multiplying the area fraction each group (A, B, C, etc), obtained through the digitizing process, with the averages of land cover characteristics from table 2.2.

2.2.4. Analysis phase

In this phase, named AN, linear regression analysis between derived land cover characteristics per pixel and variability of derived NDWI parameter between pixels was carried out using several approaches;

 First, the regression between a single cover type characteristic versus a single derived NDWI parameter.

 Second, the regression between a joint cover type characteristic (tree and shrub layer together, stone, litter, and soil together) versus a single derived NDWI parameter.

 Third, the regression analysis also applied for a land cover type characteristics. Instead of tree cover, the first layer defined as eucalyptus forest, needle leaves forest.

 Fourth, multiple regression between a single cover type characteristic versus the combination of all four derived NDWI parameters

 Fifth, as step 1-3 but this time using the original NDWI values (without classification) All regression was executed using R statistical software.

Adjusted R-square is selected to be use to report the amount of explained variability of each regression.

This adjusted R-square was reported in a table and checked for the 95% confidence level.

(26)

Adjusted R-squared is calculated by using following equation :

1 1

Where,

R² : coefficient determination n : number of samples

p : number of regressor

Thus, it is possible to have a negative value of the coefficient adjusted R-square. In the result table, this negative value of coefficient adjusted R-square is reported as “ 0 “.

(27)

18

3. RESULTS AND DISCUSSIONS

3.1. Data preparation phase

3.1.1. Generating land cover maps derived through hypertemporal NDVI and NDWI imagery

Figure 3.1 and 3.2 show the separability analysis for NDVI and NDWI hypertemporal analysis using ISODATA clustering in ERDAS for a pre-defined 90 classes (from class 10 to class 100). The figures shows the average separability in the left Y axis (left side with red colour), while the minimum separability shows in the right Y axis (right side with blue colour). The analysis determined that the best number of classes is 60 for NDVI and 45 for NDWI (black circle), because they are located at the distinct peak of average and minimum seperability. In order to keep lesser number of classes, Class 45 is selected for NDWI though the minimum seperability is not at its maximum.

Figure 3.1. Separability analysis for NDVI; the left Y axis shows the average separability (in red colour) and the right Y axis shows the minimum separability (in blue colour).

Figure 3.2. Separability analysis for NDWI; the left Y axis shows the average separability (in red colour) and the right Y axis shows the minimum separability (in blue colour).

0 30 60 90 120 150 180 210

1000 2000 3000 4000 5000 6000 7000 8000

10 20 30 40 50 60 70 80 90 100

Average Minimum

average separability minimum separability

Pre‐defined number of classes

0 60 120 180 240 300 360

0 1000 2000 3000 4000 5000 6000

10 20 30 40 50 60 70 80 90 100

Average Minimium

average separability minimum separability

Pre‐defined number of classes

(28)

Figure 3.3 illustrates the land cover map derived from hyper-temporal NDVI classification and figure 3.4 shows the NDWI classification. NDVI map units appear to be more scattered than the NDWI units.

NDVI presents more variability on land cover than NDWI. Overlying the NDVI 60 ISODATA classification and NDWI 45 ISODATA classification resulted in 1100 combinations, indicating their poor correlation. The legend of the maps are given in Annex 2.

Some polygons show similarities in shape in both maps, e.g. the one identified using a black circle. This means that the information given by both indices is similar. Based on CORINE (Coordination of information on the environment of the European Environmental Agency (EEA)) land cover map of the year 2000 (CLC 2000), these polygons are classified as rice and irrigated agriculture. This confirmed the results of (Xiao et al., 2005; Xiao et al., 2006 & Sari et al. 2010); NDWI has potential to identify paddy rice fields due to the specific phenology of the irrigated crops.

Figure 3.3. Land cover map derived hyper-temporal NDVI 60 classes images classification

Figure 3.4. Land cover map derived hyper-temporal NDWI 45 classes images classification

(29)

20

3.1.2. Generating derived NDWI parameter maps

Figure 3.5 shows the derived NDWI parameter maps of maximum, LoC, minimum, and amplitude.

Annex 3 shows the results of derived NDWI parameters. The conversion approach generalized some NDWI classes. After the conversion, the NDWI map contained 30 different values for maximum, 14 values for LoC, 27 values for minimum, and 28 values for amplitude (see Figure 3.5).

The maximum parameter map shows that the lowest maximum value appears in the eastern part and it gradually increases to the west (Figure 3.5). The situation is most likely corresponding to precipitation amounts. Some part in the eastern area have a combination of a low maximum and a low minimum which means this part is drier than the others.

Figure 3.5. Maps of derived NDWI parameters

(30)

Figure 3.5 shows that the study area has little variability in LoC; overall, the area seems to have a high value. In the eastern part, this high value of length of cycle has a small fluctuation in terms of moisture, as shown by the amplitude map.

3.2. Fieldwork phase

3.2.1. Selecting area for fieldwork

Selecting out NDVI with areas ≥ 50 square kilometres gives 195 potential polygons of NDVI. The standard deviations of NDWI parameters of these polygons ranged from a minimum of 0 to a maximum of 3.14 (Annex 4). After several adjustments made during the fieldwork, 14 polygons were selected as sampling areas representing 12 NDVI units (Table 3.1). These polygons were considered as having a high standard deviation in NDWI parameters. The highlighted rows show the three selected NDVI units for further analysis. Distribution of the polygons is shown in Figure 3.6.

According to CLC 2000, two polygons out of 14 represent agriculture (with yellow background in Figure 3.6), while the rest are representing non agriculture.

Table 3.1. Standard deviation (Stdev) of the 14 NDVI polygons selected for fieldwork NDVI unit Stdev of

Max Param Stdev of

LoC Param Stdev of

Min Param Stdev of

Amp Param Sum of Stdev of all Param

4 0.19 0.25 0.13 0.21 0.78

5 0.24 0.27 0.22 0.11 0.85

7 0.16 0.31 0.17 0.07 0.7

8 0.15 0.23 0.11 0.14 0.63

13 0.28 0.26 0.19 0.24 0.98

19 0.16 0.14 0.12 0.36 0.78

45 0.16 0.33 0.15 0.31 0.96

48 0.09 0.21 0.15 0.2 0.64

50 0.13 0.68 0.21 0.24 1.26

53 0.18 0.17 0.23 0.22 0.8

53 0.14 0.22 0.17 0.26 0.79

57 0.18 0.81 0.21 0.42 1.61

57 0.13 0.83 0.18 0.1 1.24

60 0.12 1 0.2 0.16 1.48

Figure 3.6. Distribution of selected NDVI map unit for fieldwork

(31)

22

3.2.2. Collecting field data

The total amount of data collected in the field is 590 samples spread over 12 NDVI units (Table 3.2). The samples comprise field estimations of cover characteristics. The distribution of these sample points is shown in Figure 3.6 (red dots).

Table 3.2. Summary table on field data collected NDVI unit Area (km2) Measured samples

4 245 46

5 490 121

7 676 132

8 76 8

13 62 19

19 279 42

45 91 54

48 306 23

50 337 37

53 63 45

53 65 ‐

57 171 49

57 274 10

60 122 4

T o t a l 590

3.3. Data processing phase

3.3.1. Creating image legend and averaging per cent cover

Three image legends were created; one for NDVI unit 19 consisting of 4 groups (Annex 5), one for NDVI unit 45 with 8 groups (Annex 6), and one for NDVI unit 57 with 8 groups (Annex 7).

Regarding the image legends, table 3.3, 3.4, and 3.5 show the summary of average percentage of cover characteristics for NDVI unit 19, unit 45, and unit 57. The summary is presented for each land cover characteristics (tree, shrub, etc). TS (tree and shrub) and SLS (soil, litter, and stone) represent the total per cent cover of the respective layers.

Table 3.3. Summary of average cover for NDVI unit 19 per group (%)

Group Tree Shrub Herb&Grass Stone Litter Soil TS SLS

A 21.88 0.25 2.50 5.06 2.44 67.88 22.13 75.38

B 12.11 0.00 15.56 9.89 4.78 57.67 12.11 72.33

C 0.00 2.44 0.38 3.38 11.69 82.13 2.44 97.19

D 5.67 3.00 0.83 3.83 2.50 85.50 8.67 91.83

A+B 17.75 0.13 10.00 8.09 3.91 66.38 17.88 78.38

Based on field data collected, Table 3.3 shows that soil is the dominant cover in this unit. In fact, this unit is an active agriculture area for olive and other annual crops. The agriculture field were left fallow in September when the fieldwork was carried out. This situation contributed to high percentage of soil.

The last row of the table represents an average cover when group A and B are integrating. These two groups represent image objects with similar pattern but different density level (Annex 5). The idea behind