**2. Material and methods**

118 Remote Sensing – Applications

determined by soil properties such as pH, electrical conductivity, salt content and exchangeable sodium percentage (Csillag et al., 1993; Farifteh et al., 2008). In this sense, soil reflectance is derived from the particular spectral behaviour of the heterogeneous combination of minerals, organic matter and soil water (Ben-Dor and Banin, 1994). Salt-

be detected by optical spectrometers since salt minerals have diagnostic spectral features occurring in the visible and near infrared (VNIR) and short-wave infrared (SWIR) spectral regions (Farifteh et al., 2008). Saline soils usually have evaporate minerals, which spectral features that can be explained by vibrational absorption due to water molecules chemically bound as part of the crystal structure (Howari et al., 2000). In this sense, the spectral differences of evaporates of single salt compounds are determinant of the type and

Remote sensing has been extensively employed in soil salinity studies. Data from aerial photography, videography, and optical, thermal, microwave or geophysical sensors has been used in soil salinity mapping (Metternich and Zinck, 2003). Perhaps, the most widely used remote sensing data in recent decades have been provided by multispectral (Landsat, SPOT, IRS, ASTER) or hyperspectral (DAIS, HyMap, AVIRIS, Hyperion) sensors in the spectral range approximately between 400 and 2500 nm. Researchers have frequently employed remote sensing data to map soil salinity with multispectral (Metternich and Zinck, 1997; Dwivedi et al., 2001; Melendez-Pastor et al., 2010a) and hyperspectral images (Dehaan and Taylor, 2002, 2003; Schmid et al., 2009, Ghrefat and Goodell, 2011). Pioneering studies in the 1970s employed air-borne and satellite-borne multispectral scanners to detect soil salinity, indicating the better capability of infrared bands over visible bands to locate saline soils and the low contribution of thermal bands to improve the delineation of saline areas (Richardson et al., 1976; Dalsted et al., 1979). Nowadays, imaging spectroscopy techniques are employed for the automatic detection of soil salinization with airborne or satellite sensor (Dehaan and Taylor, 2002, 2003; Dutkiewicz et al., 2009; Schmid et al., 2009; Weng et al., 2009; Melendez-Pastor et al., 2010a; Ghrefat and Goodell, 2011). Imaging spectroscopy deals with the mapping of ground materials by detecting and analysing reflectance/absorbance features in hyperspectral (or multispectral) images (Clark, 1999). Imaging spectroscopy adds a new dimension of remote sensing by expanding point spectrometry into a spatial domain and under field conditions, which is a very good approach for the study of soil properties

The aim of this chapter is the application of remote sensing for the study of soil salinity of an agricultural area in southeast coast of the Iberian Peninsula. Different digital image processing techniques were applied to satellite multispectral images (Landsat TM). 'Conventional' hard classification techniques were combined with spectral mixture analysis and soil properties to achieve a better understanding of the soil salinization process in the

Multispectral satellite images such as those obtained by the Landsat program provide low or free cost worldwide coverage for four decades. Moreover, salinization problems are concentrated in arid and semi-arid regions, often in developing countries with few economic resources. Although there are more advanced sensors that can provide a more precise quantification of the extent of soil salinity (e.g. hyperspectral), their high cost difficult its

2- and HCO3-) can

affected soils cations (Na+, Mg2+, K+, and Ca2+) and anions (Cl-, SO42-, CO3

mineralogy of the soils (Howari et al., 2000).

(Ben-Dor et al., 2009).

study area.

This study will evaluate the applicability of various remote sensing techniques for studying salinization processes in an agricultural coastal area. One of the greatest difficulties in the application of remote sensing techniques to the study area is the fragmentation of the territory by the existence of small plots and buildings that create a dispersed mixture of spectral signals to the scale of a moderate spatial resolution multispectral remote sensing image as those acquired by the Landsat Thematic Mapper sensor. This difficulty motivates the need to evaluate various techniques and methodological approaches to carry out this study as necessary to help monitoring the processes of salinization.

Representative soils of the area were sampled and their properties were characterized at the laboratory by standard methods. Predominant land cover classes at the soil sampling plots and at additional land cover validation points were identified. Land cover is a fundamental variable that impacts on and links many parts of the human and physical environments (Foody, 2002) with a great influence on soil properties (Caravaca et al., 2002; Majaliwa et al., 2010; Biro et al., 2011). Both kinds of information in a GIS database were included. In this sense, the effect of land cover on soil properties was statistically evaluated. Then, multispectral images were employed for a hard land cover mapping with a supervised approach using the k-nearest-neighbour classifier. Accuracy assessment methods highlighted the need to employ a mixed pixel focus to deal with the particularities of the study area. Spectral unmixing techniques allowed the identification of representative spectral endmembers and the obtainment of their corresponding fraction images. Finally, fraction endmembers were employed to characterize land cover classes and to predict soil properties with various statistical methods.

#### **2.1 Description of the study area**

The study area is located in a coastal zone of Southeast Spain, in the province of Alicante. It is located around 38.14°N and 0.73ºW, at the south of the cities of Elche and Alicante. The study area (Figure 1) comprises alluvial plains resulting from the accumulation of sediments from the Segura and Vinalopó rivers. During most of the Holocene (~10,000 years ago to present) the study area was a large lagoon (Blázquez, 2003). In the last centuries, the ancient lagoon was transformed into an irrigated agricultural land draining the wetland. Nowadays, this area is a mixture of small-size cities, coastal urban areas, scattered residential houses, irrigated crops and isolated and scattered wetlands. The perimeter of the study area was delimited according to natural or man-made features in order to enclose a large coastal plain area primarily occupied by irrigated agricultural activities. The study area lies in the north with the natural parks of *El Hondo* and the *Salinas de Santa Pola*. Both natural areas are wetlands included in the RAMSAR list of wetlands of international importance. The east and south boundaries are the *Sierra del Molar* and the *Segura River* respectively. Urban areas and sclerophyllous vegetation mainly occupy the Sierra del Molar, while the Segura River is the most important watercourse in southeast of Iberian Peninsula providing water for irrigation agriculture and to fill the reservoirs that currently comprise

Mapping Soil Salinization of Agricultural Coastal Areas in Southeast Spain 121

Fig. 2. Daily values of precipitation, mean temperature and evapotranspiration for the hydrological year 2010-11 at Catral station. Source data from the Spanish Ministry of

Predominant soil classes are Entisols according to the Soil Taxonomy (Soil Survey Staff, 2006) but affected by agriculture practices along years. They are characterized by a massive presence of carbonates and soluble salt content. In the studied area, irrigation is essential to support agriculture. The water deficit during several months requires irrigation while low quality water is used in the poorly drained soils of these coastal plains, being soil salinization an environmental problem. Thus, the study area soils are subjected to severe risk of physical, chemical and biological degradation (De Paz et al., 2006) that endanger

Field survey was done in the late spring and summer months of the hydrological year 2010- 2011 to collect soil samples and identify land cover classes. An extensive soil sampling was done, and 116 samples were collected and geographically referenced. Samples were obtained from the upper 5 cm as solar radiation in VNIR spectral range has limited penetration capabilities. Soil samples were dried at room temperature and a 2 mm sieve was used to separate the fine fraction to be analysed. Analysed soil characteristics included in the study were electrical conductivity (EC) (1:5 w/v water extraction), pH and organic carbon (OC) by wet chemical oxidation (Walkley and Black, 1934) with potassium

A land cover validation campaign was also conducted along with the soil survey in order to allow accuracy assessment of generated land cover maps. Land cover validation points were randomly generated in a GIS and a database with the land cover category generated. A total of 205 land cover validation points were identified, combining field observation and recent aerial orto-photography (0.5 m of spatial resolution). Land cover classes identified in the

Environment and Rural and Marine Affairs (MARM).

dichromate oxidation (Nelson and Sommers, 1982).

agriculture sustainability.

**2.2 Field survey** 

Fig. 1. Study area with the Landsat scene (false-colour composite RGB:742) and superimposed cartographic information (soil samples, urban areas, natural parks and roads).

the wetland of *El Hondo*. The western boundary of the study area is a motorway that cuts north to south the floodplain.

This coastal region has a semiarid Mediterranean climate, with a mean annual rainfall of less than 300 mm and a mean annual temperature of 17 ºC and defined by the Köppen climate classification system as *Bsk* class (dry climate with a dry season in summer and a mean annual temperature about 18 ºC). The climate is arid or semiarid according to the aridity index of Martonne (De Martonne, 1926) and the aridity index of UNEP (1997) respectively. Figure 2 shows the daily climatic diagram of mean temperature, precipitation and evapotranspiration (by the Penman-Monteith method) for the hydrological year 2010-2011 (from October to September) at Catral meteorological station. Mean daily temperature (blue line) varies from approximately 9ºC in winter to more than 25ºC in summer. Rain events (red bars) mainly occurred from December 2010 to May 2011 with total accumulated precipitation of 182 mm. This very scarce precipitation joint with an accumulated evapotranspiration of 1115 mm implied that the hydrological year was very dry.

Fig. 2. Daily values of precipitation, mean temperature and evapotranspiration for the hydrological year 2010-11 at Catral station. Source data from the Spanish Ministry of Environment and Rural and Marine Affairs (MARM).

Predominant soil classes are Entisols according to the Soil Taxonomy (Soil Survey Staff, 2006) but affected by agriculture practices along years. They are characterized by a massive presence of carbonates and soluble salt content. In the studied area, irrigation is essential to support agriculture. The water deficit during several months requires irrigation while low quality water is used in the poorly drained soils of these coastal plains, being soil salinization an environmental problem. Thus, the study area soils are subjected to severe risk of physical, chemical and biological degradation (De Paz et al., 2006) that endanger agriculture sustainability.

#### **2.2 Field survey**

120 Remote Sensing – Applications

Fig. 1. Study area with the Landsat scene (false-colour composite RGB:742) and superimposed cartographic information (soil samples, urban areas, natural parks and

evapotranspiration of 1115 mm implied that the hydrological year was very dry.

the wetland of *El Hondo*. The western boundary of the study area is a motorway that cuts

This coastal region has a semiarid Mediterranean climate, with a mean annual rainfall of less than 300 mm and a mean annual temperature of 17 ºC and defined by the Köppen climate classification system as *Bsk* class (dry climate with a dry season in summer and a mean annual temperature about 18 ºC). The climate is arid or semiarid according to the aridity index of Martonne (De Martonne, 1926) and the aridity index of UNEP (1997) respectively. Figure 2 shows the daily climatic diagram of mean temperature, precipitation and evapotranspiration (by the Penman-Monteith method) for the hydrological year 2010-2011 (from October to September) at Catral meteorological station. Mean daily temperature (blue line) varies from approximately 9ºC in winter to more than 25ºC in summer. Rain events (red bars) mainly occurred from December 2010 to May 2011 with total accumulated precipitation of 182 mm. This very scarce precipitation joint with an accumulated

roads).

north to south the floodplain.

Field survey was done in the late spring and summer months of the hydrological year 2010- 2011 to collect soil samples and identify land cover classes. An extensive soil sampling was done, and 116 samples were collected and geographically referenced. Samples were obtained from the upper 5 cm as solar radiation in VNIR spectral range has limited penetration capabilities. Soil samples were dried at room temperature and a 2 mm sieve was used to separate the fine fraction to be analysed. Analysed soil characteristics included in the study were electrical conductivity (EC) (1:5 w/v water extraction), pH and organic carbon (OC) by wet chemical oxidation (Walkley and Black, 1934) with potassium dichromate oxidation (Nelson and Sommers, 1982).

A land cover validation campaign was also conducted along with the soil survey in order to allow accuracy assessment of generated land cover maps. Land cover validation points were randomly generated in a GIS and a database with the land cover category generated. A total of 205 land cover validation points were identified, combining field observation and recent aerial orto-photography (0.5 m of spatial resolution). Land cover classes identified in the

Mapping Soil Salinization of Agricultural Coastal Areas in Southeast Spain 123

used for the correction. The nearest neighbour resampling method was selected because it ensures that the original (raw) pixel values are retained in the resulting output image, which is an important requirement in any change detection analysis (Mather 2004). The maximum allowable root mean square error (RMSE) of the geometric correction was less than half a

Atmospheric correction involves the estimation of the atmospheric optical characteristics at the time of image acquisition before applying the correction to the data (Kaufman, 1989). This type of correction is a pre-requisite in many remote sensing applications such as in classification and change detection procedures (Song et al., 2001). Radiometric calibration was applied prior to the atmospheric correction. The conversion of raw digital numbers (DNraw) of Landsat level 1 (L1) image products to at-satellite radiance values (Lsat) required the application of current re-scaling values (Chander et al., 2010) by applying the following

pixel, a reference value frequently cited (Townsend and Walsh 2001; Jensen, 2005).

MAX MIN sat (DN) MIN <sup>255</sup>

 

Where *L*sat is at-satellite radiance [W/(m2 sr m)]; *L*MIN is the spectral radiance that is scaled to *Q*calmin [W/(m2 sr m)] (*Q*calmin is the minimum quantized calibrated pixel value, i.e. DN=0, corresponding to *L*MIN ); *L*MAX is the spectral radiance that is scaled to *Q*calmax [W/(m2 sr m)] (*Q*calmax is the maximum quantized calibrated pixel value, i.e. DN=255, corresponding to *L*MAX ); and DN are digital numbers of the L1 image product. Surface

1996). Path radiance (*L*p) values were computed by using the equation reported in Song et al. (2001) that assumes 1% surface reflectance for dark objects (Chavez Jr, 1989, 1996; Moran et

Image classification procedures aim to automatically categorize all pixels in an image into land cover classes or themes (Lillesand et al., 2003). Thematic mapping from remotely sensed data can be defined as grouping together cases (pixels) by their relative spectral similarity (unsupervised component) with the aim of allocating cases based on their similarity to a set of predefined classes that have been characterized spectrally (supervised component) (Foody 2002). Multispectral images (like Landsat TM scenes) are frequently used to perform the classification based on spectral pattern recognition methods that exploits the pixel-by-pixel spectral information as the basis for automated classification (Lillesand et al., 2003). In this study, a supervised land cover classification of the Landsat TM image was performed with a *k*nearest-neighbour clustering algorithm to obtain a discrete or 'hard' categorical land cover map for the study area. *K*-nearest-neighbour (KNN) classifier searches away from the pixel to be classified in all directions of the spectral space until it encounters *k* user-specified training pixels and then assigns the pixel to the class with the majority of pixels encountered (Jensen, 2005). KNN algorithm has been successfully applied for land cover classification with remote sensing data (Franco-Lopez et al., 2001; Haapanen et al., 2004; Blanzieri and Melgani, 2008). Land cover classes assigned to the soil plots were employed in the training stage of the

 

) were computed by using the image based COST method (Chavez Jr,

(1)

r) was estimated according to the

*L L L L* 

expression (Chander and Markham, 2003; Chander et al., 2010):

reflectance values (

equation given in Kaufman (1989).

**2.4 Land cover classification** 

al., 1992). The optical thickness for Rayleigh scattering (

study area were: water bodies, seasonal or permanent crops, saltmarshes and misused agricultural field that tends to be saltmarshes, palm groves, marshes with almost permanent inundation, and anthropic areas (Table 1).


Table 1. Descriptions of land cover classes identified in the study area.

Land cover categories at soil sampling points were also identified and included along with soil properties in a GIS database for the land cover classification training stage and for further spatial analyses. Note that land cover (i.e. biophysical materials found on the land) and land use (i.e. how the land is being used by human beings) (Jensen, 2007) are different terms but often used together or interchangeably. In this chapter, we adopt the term land cover because we are interested in knowing about the biophysical characteristics of the study area, but the knowledge of both land use and land cover are important for land planning and land management activities (Lillesand et al., 2003).

#### **2.3 Satellite imagery preprocessing**

Remote sensing data were acquired by the Thematic Mapper (TM) sensor on-board the Landsat 5 satellite. Meteorological conditions and the satellite pass over the study area conditioned the date of image acquisition. A scene acquired on 28th June 2011 (path 199 row 33) was employed for analyses. A vertical black line on Figure 2 indicates the time of acquisition of the scene. No rain events happened 16 days prior to the scene acquisition date. Typically summer meteorological conditions without cloud coverage and high temperature were registered on the date of image acquisition, and thus the image quality was optimal.

Satellite image preprocessing included geometric and atmospheric corrections with the aim to ensure the spatial comparability with other data sources and to obtain at-ground reflectance pixel spectra, respectively. Various georeferenced data types were used for the geometric correction: aerial orthophotos (0.5 m of pixel resolution) and digital cartography (scale = 1:10000). The Landsat 5 TM scene was geometrically corrected using Ground Control Points (GCP) identified on the orthophotos and cartographic maps. A quadratic mapping function of polynomial fit and the nearest neighbour resampling method were

study area were: water bodies, seasonal or permanent crops, saltmarshes and misused agricultural field that tends to be saltmarshes, palm groves, marshes with almost permanent

Arable land 2 Herbaceous (e.g. alfalfa, barley) and horticultural (e.g. melon,

Saltmarsh 5 Halophyte vegetation (e.g. *Salicornia sp*., *Suaeda sp*., *Limonium sp*.,

Land cover categories at soil sampling points were also identified and included along with soil properties in a GIS database for the land cover classification training stage and for further spatial analyses. Note that land cover (i.e. biophysical materials found on the land) and land use (i.e. how the land is being used by human beings) (Jensen, 2007) are different terms but often used together or interchangeably. In this chapter, we adopt the term land cover because we are interested in knowing about the biophysical characteristics of the study area, but the knowledge of both land use and land cover are important for land

Remote sensing data were acquired by the Thematic Mapper (TM) sensor on-board the Landsat 5 satellite. Meteorological conditions and the satellite pass over the study area conditioned the date of image acquisition. A scene acquired on 28th June 2011 (path 199 row 33) was employed for analyses. A vertical black line on Figure 2 indicates the time of acquisition of the scene. No rain events happened 16 days prior to the scene acquisition date. Typically summer meteorological conditions without cloud coverage and high temperature were registered on the date of image acquisition, and thus the image quality

Satellite image preprocessing included geometric and atmospheric corrections with the aim to ensure the spatial comparability with other data sources and to obtain at-ground reflectance pixel spectra, respectively. Various georeferenced data types were used for the geometric correction: aerial orthophotos (0.5 m of pixel resolution) and digital cartography (scale = 1:10000). The Landsat 5 TM scene was geometrically corrected using Ground Control Points (GCP) identified on the orthophotos and cartographic maps. A quadratic mapping function of polynomial fit and the nearest neighbour resampling method were

Palm groves 6 Palm trees plantations and nurseries, mainly from *Phoenix* 

Marsh 7 *Phragmites australis* dominated wetland vegetation Man-made/urban 8 Urban areas, roads, farms or industrial areas Table 1. Descriptions of land cover classes identified in the study area.

**Land cover ID Features**  Water 1 Wetlands water tables and irrigation ponds

broccoli) crops Permanent crops 3 Fruit trees (e.g. orange, lemon, pomegranate) Fallow/abandoned 4 Fallow or recently abandoned agricultural land.

*Halocnemum sp*.)

*dactilifera*

planning and land management activities (Lillesand et al., 2003).

**2.3 Satellite imagery preprocessing** 

was optimal.

inundation, and anthropic areas (Table 1).

used for the correction. The nearest neighbour resampling method was selected because it ensures that the original (raw) pixel values are retained in the resulting output image, which is an important requirement in any change detection analysis (Mather 2004). The maximum allowable root mean square error (RMSE) of the geometric correction was less than half a pixel, a reference value frequently cited (Townsend and Walsh 2001; Jensen, 2005).

Atmospheric correction involves the estimation of the atmospheric optical characteristics at the time of image acquisition before applying the correction to the data (Kaufman, 1989). This type of correction is a pre-requisite in many remote sensing applications such as in classification and change detection procedures (Song et al., 2001). Radiometric calibration was applied prior to the atmospheric correction. The conversion of raw digital numbers (DNraw) of Landsat level 1 (L1) image products to at-satellite radiance values (Lsat) required the application of current re-scaling values (Chander et al., 2010) by applying the following expression (Chander and Markham, 2003; Chander et al., 2010):

$$L\_{\rm sat} = \left(\frac{L\_{\rm MAX\lambda} - L\_{\rm MIN\lambda}}{255}\right) \text{(DN)} + L\_{\rm MIN\lambda} \tag{1}$$

Where *L*sat is at-satellite radiance [W/(m2 sr m)]; *L*MIN is the spectral radiance that is scaled to *Q*calmin [W/(m2 sr m)] (*Q*calmin is the minimum quantized calibrated pixel value, i.e. DN=0, corresponding to *L*MIN ); *L*MAX is the spectral radiance that is scaled to *Q*calmax [W/(m2 sr m)] (*Q*calmax is the maximum quantized calibrated pixel value, i.e. DN=255, corresponding to *L*MAX ); and DN are digital numbers of the L1 image product. Surface reflectance values () were computed by using the image based COST method (Chavez Jr, 1996). Path radiance (*L*p) values were computed by using the equation reported in Song et al. (2001) that assumes 1% surface reflectance for dark objects (Chavez Jr, 1989, 1996; Moran et al., 1992). The optical thickness for Rayleigh scattering (r) was estimated according to the equation given in Kaufman (1989).

#### **2.4 Land cover classification**

Image classification procedures aim to automatically categorize all pixels in an image into land cover classes or themes (Lillesand et al., 2003). Thematic mapping from remotely sensed data can be defined as grouping together cases (pixels) by their relative spectral similarity (unsupervised component) with the aim of allocating cases based on their similarity to a set of predefined classes that have been characterized spectrally (supervised component) (Foody 2002). Multispectral images (like Landsat TM scenes) are frequently used to perform the classification based on spectral pattern recognition methods that exploits the pixel-by-pixel spectral information as the basis for automated classification (Lillesand et al., 2003). In this study, a supervised land cover classification of the Landsat TM image was performed with a *k*nearest-neighbour clustering algorithm to obtain a discrete or 'hard' categorical land cover map for the study area. *K*-nearest-neighbour (KNN) classifier searches away from the pixel to be classified in all directions of the spectral space until it encounters *k* user-specified training pixels and then assigns the pixel to the class with the majority of pixels encountered (Jensen, 2005). KNN algorithm has been successfully applied for land cover classification with remote sensing data (Franco-Lopez et al., 2001; Haapanen et al., 2004; Blanzieri and Melgani, 2008). Land cover classes assigned to the soil plots were employed in the training stage of the

Mapping Soil Salinization of Agricultural Coastal Areas in Southeast Spain 125

error associated to the estimation of proportions for each pixel *i*, *j*. The Least Square Mixing Model proposed by Shimabukuro and Smith (1991) is commonly used to resolve linear spectral mixture models. The method proposed by Shimabukuro and Smith (1991) assumes two initial restrictions for the computation of the proportions of spectrally pure endmembers. The first one implies that pure endmember proportions must range between 0 and 1. This means that the proportions of the components are normalized to a common

The second restriction is that the sum of the fractions for every component is equal to the total pixel surface. In this way, it is quite simple to express the individual contribution or

The choice of a LSU model must consider both the landscape of the test site and the ability of the model to depict the structure, shape and distribution of the basic landscape components (Ferreira et al. 2007). Well-chosen endmembers not only represent materials found in the scene, but provide an intuitive basis for understanding and describing the information in the image (Adams and Gillespie 2006). Endmembers were obtained after applying a spatial and spectral remote sensing data dimensionality reduction with the minimum noise fraction (MNF) and pixel purity index (PPI) techniques, respectively. The MNF is used to detect the inherent dimensionality of image data, segregating noise from the signal in the data and reducing computational requirements for subsequent processing tasks (Boardman and Kruse, 1994). The MNF as modified from Green et al. (1988) consists in two steps: 1) applying a transformation, based on an estimated noise covariance matrix to decorrelate and rescale the noise in the data (noise has unit variance and no band-to-band correlations); and 2) performing a standard principal component transformation of the noise-whitened data. A final dataset of coherent and almost noise-free bands are selected from the MNF output and can be used for subsequent processing steps. Pixel Purity Index (PPI) is a procedure for finding the most spectrally pure (extreme) pixels that typically correspond to mixing endmembers in multispectral and hyperspectral images (ITT VIS, 2008). PPI is computed by repeatedly projecting *n*-dimensional scatterplots onto a random unit vector; the extreme pixels in each projection (those pixels that fall onto the ends of the unit vector) are recorded and the total number of times each pixel is marked as extreme is noted. The selection of extreme pixels corresponding to analogous surface features is complex due to the great number of pixels typically found in remote sensing image data. The *n*-dimensional visualizer implemented in ENVI software (ITT Visual Information Solutions) is a tool to locate, identify, and cluster the purest pixels and most extreme spectral responses in a data set. The distribution of these points in *n*-space can be used to estimate the number of spectral endmembers and their pure spectral signatures (Boardman, 1993). Three endmembers were used in the LSU model of the study area, namely green vegetation (GV), non-photosynthetic vegetation (NPV) and shade (S). The GV endmember represents the signature of green dense vegetation, the NPV endmember is the signature of bare soil or sparse non-photosynthetic vegetation, and the shade endmember represents the signature of

, , 0 1 *Fijm* (3)

, , 1, <sup>1</sup> *m p Fijm* (4)

range of potential values. The following expression summarizes this first restriction:

fraction of an endmember in relation to the total reflectance of the pixel.

dark pixels and water bodies.

algorithm. Major urban areas were digitized with a GIS and masked-out of the supervised classification procedure since urban areas induce a great spectral confusion. Water areas training points were also included in the training dataset.

The land cover validation database was employed to evaluate the performance of the classification. Land cover map accuracy assessment was quantified with statistical methods such as the error matrix and the kappa statistic. The error matrix is a square array of numbers organized in rows and columns that express the number of sample units (i.e. pixels) assigned to a particular category relative to the actual category as indicated by the reference data (Congalton, 2004). Reference data are in the columns while the rows indicate the map categories to be assessed. This form of expressing accuracy as an error matrix is an effective way to evaluate both errors of inclusion (commission errors) and errors of exclusion (omission errors) present in the classification as well as the overall accuracy (Congalton et al., 1983). In addition to the error matrix, the Kappa coefficient developed by Cohen (1960) was employed to quantify the accuracy of the land cover map. Cohen's Kappa (or KHAT) is a measure of agreement for nominal scales based on the difference between the actual agreement of the classification (i.e., agreement between computer classification and reference data as indicated by the diagonal elements) and the chance agreement, which is indicated by the product of the row and column marginal (Congalton et al., 1993).

#### **2.5 Spectral unmixing**

A mixed pixel results when a sensor's Instantaneous Field of View (IFOV) includes more than one land cover type on the ground (Lillesand et al., 2003). The spectrum of a single pixel is a complex measurement that integrates the radiant flux from all the spatially unresolved materials in the IFOV, regardless of whether or not we know their identities (Adams and Gillespie, 2006). Spectral mixture analysis (SMA) has been developed as a method to transform the reflectance in the bands of multispectral images to fractions of reference endmembers, which are reflectance spectra of well-characterized materials that mix to produce spectra equivalent to those of pixels of interest in the image (Adams et al., 1995). As part of SMA techniques, linear spectral unmixing (LSU) models tread the radiation recorded by a sensor as the result of a linear mixture of spectrally pure endmember radiances (Small and Lu 2006). This method is based on the assumptions that: 1) the recorded radiation by the sensor for each pixel is limited to the sensor's IFOV, and assumes no influences by reflected radiation from neighbouring pixels (Settle and Drake 1993), 2) the overall global radiance is proportional to the surface occupied by each land cover type, and 3) the spectrally pure endmembers are valid for the whole study area (Quarmby et al. 1992). LSU models describe radiation reflected by an individual pixel (*i,j*) of a band *k* as the result of the product of reflectance for each land cover type by their respective mixture fraction plus an additional associated error for each pixel. The general expression of the model is presented in the following equation:

$$
\rho\_{i,j,k} = \sum\_{m=1,p} F\_{i,j,m} \, \rho\_{m,k} + e\_{i,j} \tag{2}
$$

Where i,j,k is the observed reflectance of a pixel for row *i*, column *j*, and band *k* ; *F*i,j,m is the proportion of component *m* of a pixel for row *i*, column *j*, for each one of the pure components; m,k is the characteristic reflectance for component *m* in band *k* ; and *e*i,j is the

algorithm. Major urban areas were digitized with a GIS and masked-out of the supervised classification procedure since urban areas induce a great spectral confusion. Water areas

The land cover validation database was employed to evaluate the performance of the classification. Land cover map accuracy assessment was quantified with statistical methods such as the error matrix and the kappa statistic. The error matrix is a square array of numbers organized in rows and columns that express the number of sample units (i.e. pixels) assigned to a particular category relative to the actual category as indicated by the reference data (Congalton, 2004). Reference data are in the columns while the rows indicate the map categories to be assessed. This form of expressing accuracy as an error matrix is an effective way to evaluate both errors of inclusion (commission errors) and errors of exclusion (omission errors) present in the classification as well as the overall accuracy (Congalton et al., 1983). In addition to the error matrix, the Kappa coefficient developed by Cohen (1960) was employed to quantify the accuracy of the land cover map. Cohen's Kappa (or KHAT) is a measure of agreement for nominal scales based on the difference between the actual agreement of the classification (i.e., agreement between computer classification and reference data as indicated by the diagonal elements) and the chance agreement, which is indicated by the product of the row and column marginal

A mixed pixel results when a sensor's Instantaneous Field of View (IFOV) includes more than one land cover type on the ground (Lillesand et al., 2003). The spectrum of a single pixel is a complex measurement that integrates the radiant flux from all the spatially unresolved materials in the IFOV, regardless of whether or not we know their identities (Adams and Gillespie, 2006). Spectral mixture analysis (SMA) has been developed as a method to transform the reflectance in the bands of multispectral images to fractions of reference endmembers, which are reflectance spectra of well-characterized materials that mix to produce spectra equivalent to those of pixels of interest in the image (Adams et al., 1995). As part of SMA techniques, linear spectral unmixing (LSU) models tread the radiation recorded by a sensor as the result of a linear mixture of spectrally pure endmember radiances (Small and Lu 2006). This method is based on the assumptions that: 1) the recorded radiation by the sensor for each pixel is limited to the sensor's IFOV, and assumes no influences by reflected radiation from neighbouring pixels (Settle and Drake 1993), 2) the overall global radiance is proportional to the surface occupied by each land cover type, and 3) the spectrally pure endmembers are valid for the whole study area (Quarmby et al. 1992). LSU models describe radiation reflected by an individual pixel (*i,j*) of a band *k* as the result of the product of reflectance for each land cover type by their respective mixture fraction plus an additional associated error for each pixel. The general expression of the model is

, , ,, , , 1,

proportion of component *m* of a pixel for row *i*, column *j*, for each one of the pure

 

i,j,k is the observed reflectance of a pixel for row *i*, column *j*, and band *k* ; *F*i,j,m is the

m,k is the characteristic reflectance for component *m* in band *k* ; and *e*i,j is the

*<sup>i</sup> <sup>j</sup> k i m p F e <sup>j</sup> m mk i <sup>j</sup>* (2)

training points were also included in the training dataset.

(Congalton et al., 1993).

**2.5 Spectral unmixing** 

presented in the following equation:

Where 

components;

error associated to the estimation of proportions for each pixel *i*, *j*. The Least Square Mixing Model proposed by Shimabukuro and Smith (1991) is commonly used to resolve linear spectral mixture models. The method proposed by Shimabukuro and Smith (1991) assumes two initial restrictions for the computation of the proportions of spectrally pure endmembers. The first one implies that pure endmember proportions must range between 0 and 1. This means that the proportions of the components are normalized to a common range of potential values. The following expression summarizes this first restriction:

$$0 \le F\_{i,j,m} \le 1\tag{3}$$

The second restriction is that the sum of the fractions for every component is equal to the total pixel surface. In this way, it is quite simple to express the individual contribution or fraction of an endmember in relation to the total reflectance of the pixel.

$$\sum\_{m=1,p} F\_{i,j,m} = \mathbf{1} \tag{4}$$

The choice of a LSU model must consider both the landscape of the test site and the ability of the model to depict the structure, shape and distribution of the basic landscape components (Ferreira et al. 2007). Well-chosen endmembers not only represent materials found in the scene, but provide an intuitive basis for understanding and describing the information in the image (Adams and Gillespie 2006). Endmembers were obtained after applying a spatial and spectral remote sensing data dimensionality reduction with the minimum noise fraction (MNF) and pixel purity index (PPI) techniques, respectively. The MNF is used to detect the inherent dimensionality of image data, segregating noise from the signal in the data and reducing computational requirements for subsequent processing tasks (Boardman and Kruse, 1994). The MNF as modified from Green et al. (1988) consists in two steps: 1) applying a transformation, based on an estimated noise covariance matrix to decorrelate and rescale the noise in the data (noise has unit variance and no band-to-band correlations); and 2) performing a standard principal component transformation of the noise-whitened data. A final dataset of coherent and almost noise-free bands are selected from the MNF output and can be used for subsequent processing steps. Pixel Purity Index (PPI) is a procedure for finding the most spectrally pure (extreme) pixels that typically correspond to mixing endmembers in multispectral and hyperspectral images (ITT VIS, 2008). PPI is computed by repeatedly projecting *n*-dimensional scatterplots onto a random unit vector; the extreme pixels in each projection (those pixels that fall onto the ends of the unit vector) are recorded and the total number of times each pixel is marked as extreme is noted. The selection of extreme pixels corresponding to analogous surface features is complex due to the great number of pixels typically found in remote sensing image data. The *n*-dimensional visualizer implemented in ENVI software (ITT Visual Information Solutions) is a tool to locate, identify, and cluster the purest pixels and most extreme spectral responses in a data set. The distribution of these points in *n*-space can be used to estimate the number of spectral endmembers and their pure spectral signatures (Boardman, 1993). Three endmembers were used in the LSU model of the study area, namely green vegetation (GV), non-photosynthetic vegetation (NPV) and shade (S). The GV endmember represents the signature of green dense vegetation, the NPV endmember is the signature of bare soil or sparse non-photosynthetic vegetation, and the shade endmember represents the signature of dark pixels and water bodies.

Mapping Soil Salinization of Agricultural Coastal Areas in Southeast Spain 127

2011). Cultivated areas (i.e. arable land and permanent crops, homogenous subgroup *b*) have lower electrical conductivity values than natural or semi-natural vegetation. The construction of drainage systems at agricultural areas to encourage the leaching for salinity control has been a traditional amelioration strategy (Qadir et al., 2000). This fact explains that marshes and saltmarshes soils have higher EC values than the other land cover classes, since they are areas with poor drainage and temporally flooded. Salinity increases when farming finishes (i.e. fallow/abandoned) because irrigation water is not available to

The pH values were slightly alkaline but significantly different for marshes. Wetland soils are characterized by the permanent or seasonal inundation of the land, promoting anaerobic conditions and thus reduced redox conditions (high concentration of H+ which implies low pH)(Reddy et al., 2000). The organic carbon content was also different depending on the type of land cover. Arable land and permanent crops soils have organic carbon content ranging from 1.46 to 2.11% that is not very high (Pérez-Sirvent et al., 2003). Opposite, wetland soils (i.e. stable saltmarshes and saltmarshes) exhibited the highest organic carbon contents. Compared to upland areas, most wetland soils show an accumulation of organic matter by the higher rates of photosynthesis in wetlands than other ecosystems and the lower rates of decomposition due to anaerobic conditions

All soil properties are significantly correlated (P<0.01) according to the Pearson bivariate correlation test applied to the full dataset. Figure 3 shows two scatterplots of the land cover classes average values (error bars represent the standard deviation) of pH and organics carbon versus electrical conductivity (EC). EC is negatively correlated with pH (R=-0.61) and positively correlated with organic carbon (R=0.34), while pH and organic carbon are negatively correlated (R=-0.32). Two sets of distinct land uses mainly dependent on the EC values are distinguished: 1) active cultures: with low EC and OC values, and 2) natural vegetation and crops of low requirements (palm groves): with high EC and OC values, increasing as the land cover is more similar to the wetland. Palm groves are the most

A land cover map was obtained with the *k*-nearest-neighbours algorithm (Figure 4). Optimum results were obtained with *k*=4. The area occupied by the land cover classes

The study area is mainly agricultural but largely occupied by urban areas. Urban/manmade areas represent 14.4% of the study area. There is a clear distinction between the northeast portion (area between the two natural parks and the *Sierra del Molar*) with large fields and less presence of buildings, and the rest of the study area, with numerous buildings scattered, villages, small-size towns and smaller parcels. Dominant land cover classes are also different at these two sectors. Close to the natural parks, there are many saltmarshes (7.09% of the study area), marshes (1.95% of the study area), palm groves (2.21% of the study area) and arable land (45.76% of the area and mainly forage, barley and melons). The other sector has a massive presence of permanent crops (16.3% of the study

(hectares and percentage of the total area) was quantified (Table 3).

area and mainly citrus trees such as orange and lemon trees).

promote salts leaching.

(Reddy et al., 2000).

halotolerant crop and require little tillage.

**3.2 Land cover classification** 

### **2.6 Statistical methods**

The possible existence of differences in soil properties based on the land cover classes was determined by the use of the analysis of variance (ANOVA). ANOVA is used to evaluate significant differences between means of independent variables. The observed variance of independent variables is partitioned into components by several explanatory variables (factors). Land cover class was the factor employed in the analysis. Post-hoc analysis was performed using Tukey method.

Relationships between fraction endmembers and soil properties were studied by the principal component analysis (PCA). PCA is a technique of data dimensionality reduction that performs an orthogonal transformation to convert potentially correlated input variables into uncorrelated variables or principal components. The first components accumulate most of the variance and therefore, the most useful information about the variables.

Regression analyses between soil properties and fraction endmembers were performed for quantitative estimation of soil properties. A linear regression analysis applying a stepwise method for variable entry and removal was the selected statistical technique. Model selection was based on the lower typical error of the estimation and minimum collinearity.
