Vol49_1_2006def 219 ANNALS OF GEOPHYSICS, VOL. 49, N. 1, February 2006 Key words decomposition rates – hyper-spectral image – NDVI – co-kriging 1. Introduction Environmental assessment at landscape and regional level claims new and efficient technolo- gies and procedures in order to take full account of both structural and functional properties of en- vironmental systems on a broad spatial extent. Nowadays, remote sensing is an essential tool for ecosystem and landscape status monitoring. Multi- and hyper-spectral remotely-sensed data are more and more frequently used for ecological data analysis, landscape unit mapping and spatial analysis of ecosystem structural organization. In the past few years, spatial data analysis has been frequently applied to remotely sensed imageries. Only recently has attention been given to the spa- tial analysis of processes at landscape level (Dugan et al., 1994; Curran, 2001). Yet coupling structures which are remotely detected with functioning still remains one major task of the environmental monitoring activities. Litter de- composition is a key process in terrestrial eco- system functioning because it controls nutrient availability in the soil, thus regulating primary production, improves soil structure and reduces soil erosion. Decomposition is a complex process Mapping litter decomposition by remote-detected indicators Letizia Sabetta (1), Nicola Zaccarelli (2), Giorgio Mancinelli (1), Stefania Mandrone (3), Rosamaria Salvatori (3), Maria Letizia Costantini (4), Giovanni Zurlini (2) and Loreto Rossi (4) (1) Laboratorio di Ecologia, Dipartimento di Scienze e Tecnologie Biologiche ed Ambientali, Università degli Studi di Leccce, Ecotekne, Lecce, Italy (2) Laboratorio di Ecologia del Paesaggio, Dipartimento di Scienze e Tecnologie Biologiche ed Ambientali, Università degli Studi di Leccce, Ecotekne, Lecce, Italy (3) Istituto sull’Inquinamento Atmosferico (IIA), CNR, Roma, Italy (4) Dipartimento di Genetica e Biologia Molecolare - Area Ecologica, Università degli Studi di Roma «La Sapienza», Roma, Italy Abstract Leaf litter decomposition is a key process for the functioning of natural ecosystems. An important limiting fac- tor for this process is detritus availability, which we have estimated by remote sensed indices of canopy green biomass (NDVI). Here, we describe the use of multivariate geostatistical analysis to couple in situ measures with hyper-spectral and multi-spectral remote-sensed data for producing maps of litter decomposition. A direct rela- tionship between the decomposition rates in four different CORINE habitats and NDVI, calculated at different scales from Landsat ETM+ multi-spectral data and MIVIS hyper-spectral data was found. Variogram analysis was used to evaluate the spatial properties of each single variable and their common interaction. Co-variogram and co-kriging analysis of the two variables turned out to be an effective approach for decomposition mapping from remote-sensed spatial explicit data. Mailing address: Dr. Letizia Sabetta, Laboratorio di Eco- logia, Dipartimento di Scienze e Tecnologie Biologiche ed Ambientali, Università degli Studi di Leccce, Ecotekne, SP Lecce-Monteroni, 73100 Lecce, Italy; e-mail: Letizia.sabet- ta@unile.it 220 Letizia Sabetta et al. that occurs through different phases and depends on many environmental factors. The rate of litter decomposition is regulated by abiotic and biotic factors, such as climate, litter quality and decom- poser activity (Swift et al., 1979; Gallardo and Merino, 1993; Aerts, 1997; Perez-Harguindeguy et al., 2000). A variety of litter quality indices based on the initial chemical composition (Aber and Melillo, 1980; Aber et al., 1990) and leaf structure parameters (Gallardo and Merino, 1993) have been used as decay predictors; more- over, works on the use of Near Infrared Re- flectance Spectroscopy (NIRS) have shown that it is a powerful and rapid method for predicting the biochemical composition of forest foliage (Card et al., 1988; Martin and Aber, 1994) and also litter decomposability (Gillon et al., 1999). According to Saunders (1976), decomposition rate is a function of both the substrate and de- composer concentrations as bimolecular second- order reaction. In forested habitats, the substrate concentration strictly depends on litter input from canopies which can be remotely sensed and estimated through the Normalized Difference Vegetation Index (NDVI, Rouse et al., 1974). In this study we present results from multi- variate geostatistical analyses used to couple in situ measures of decomposition rates with hy- per-spectral and multi-spectral remotely sensed data in terrestrial ecosystems. Relationships be- tween decomposition rates, estimated for four CORINE habitats, and NDVI values estimated at various scales from multi-spectral Landsat ETM+ and hyper-spectral MIVIS data, were in- vestigated. Variogram and co-kriging analyses were carried out to evaluate the spatial proper- ties of each ecosystem variable, their common interaction and for mapping them from remote- sensed spatially explicit data. 2. Materials and methods Study area – The study was carried out in the catchment-basin of Lake Vico, a volcanic area located in Central Italy about 50 km north of Rome (42°19lN, 12°10lE). Since 1982 this area (about 3200 ha) has been part of a Regional Re- serve, which includes the lake (510 m a.s.l.) and Mount Venere (851 m a.s.l.). Plant communities follow a typical altitudinal gradient from reed thickets (Phragmites australis) located along the lake shores and hazelnut cultivated areas to mixed forests. In this study we have considered four plant communities referring to four CORINE habitats (CooRdination de l’INforma- tion sur l’Environnament; European Union/Di- rectorate General XI, 1991): hazelnut areas (dominated by Corylus avellana, CORINE code 83.1), chestnut forests (dominated by Castanea sativa, CORINE code 41.9), turkey oak forests (dominated by Quercus cerris, CORINE code 41.74) and beech forests (dominated by Fagus sylvatica, CORINE code 41.181). Decomposition studies – Plant litter break- down was studied using the litterbag technique in 20 sampling sites (5 sampling sites for each CORINE habitat). Sampling sites were georefer- enced through a GPS system. Fifty litterbags (mesh-bags of 0.5 cm mesh size containing 3 dry-g of the dominant species leaves) were ran- domly placed at each sampling site to simulate the natural heaps of detritus. Ten litterbags were retrieved monthly from each site from December to October. Samples were carefully and separate- ly placed in polythene boxes and then brought to laboratory, where their ash-free dry masses were measured after leaf cleaning, oven-drying (at 60°C for 72 h) and ignition (at 800°C for 3 h). Remote sensing data – The remote sensed imageries came from Landsat ETM+ and MIVIS (Daedalus AA5000 Multispectral Infrared & Visible Imaging Spectrometer, Italian National Research Council – LARA Project) and they were acquired by the end of the vegetative peri- od, before leaves falling, in order to evaluate the detritus input. MIVIS scanner is a hyper-spectral airborne imaging system with 102 spectral bands covering visible and near infrared (0.43-0.83 nm, 1st sensor), middle infrared (1.15-1.55 nm, 2nd sensor; 1.98-2.50 nm, 3rd sensor) and ther- mal infrared (8.21-12.70 nm, 4th sensor) regions of the electromagnetic spectrum. MIVIS images had different spatial resolutions (4, 8, 10 m) cor- responding to different fly heights (2000, 4000 and 5000 m a.s.l.). The resolution interval con- sidered in this study thus ranged from 4 to 30 m (from Landsat images). MIVIS data were ac- quired on 22 September 2000, at 12.00 am local time. All MIVIS images were geocoded by the 221 Mapping litter decomposition by remote-detected indicators Delaunay triangulation procedure using a com- mon set of ground referenced points and or- tophotos of the area (Richards, 1994). Maximum root mean square error (RMS) was 0.3. Landsat data were acquired on 16 August 2000, which was the only available scene with no cloud cov- er, and it was geocoded by an affine transforma- tion of first order with an RMS error of 0.2. Landsat data radiometric calibration and conver- sion to radiance were performed by applying stored parameters derived by the calibration pa- rameter file supplied with the original data and using Chander and Markham formulas (USGS, 2004). MIVIS data were acquired as radiance values by the provider. Basic atmospheric cor- rections were applied to all images: Landsat da- ta were corrected using a Dark Object Subtrac- tion (DOS) approach (Chavez, 1988) with Lake Vico as a target; MIVIS data were corrected by an «empirical line calibration» procedures based on in situ collected spectra for the vegetation cover of the four CORINE habitats and bare soil (Kruse et al., 1990). Data analyses – Leaf mass losses were fit- ted to the simple exponential model of decom- position, Mt = M0 e−Kt where M0 is the initial leaf AFDM (Ash Free Dry Mass), Mt is remaining leaf AFDM at time t (in days) and K is a break- down coefficient expressed in days−1 (Olson, 1963). The half-life, i.e. the time necessary to reduce detritus mass to 50% of its initial AFDM, was expressed as ln(2)/K. NDVI values were calculated using band number 3 (0.63-0.69 nm) and 13 (0.67-0.69 nm) as red channels and band number 4 (0.75- 0.90 nm) and 19 (0.79-0.81 nm) as infrared channels from Landsat and MIVIS data, respec- tively. In plots of 60 × 60 m2, centered on each sampling site, we determined the mean NDVI values at each spatial resolution. Plot dimen- sions were based on the position accuracy of Landsat data, which was at least of 30 m, i.e. one pixel, if the geocoding process was per- formed with RMS error lower than 1. Mean plot NDVI values were calculated for each pixel of each image of the study area. The relation between the litter half-life time and NDVI was investigated applying both clas- sical and spatial statistics: classical linear re- gression analysis and ANCOVA were per- formed to evaluate the scale information decay while a geostatistical approach based on vari- ogram evaluation was carried out to model spa- tial dependency within scales. Spatial continuity among observations of a given set of variables may be characterised by a variogram or a co-variogram, which reveals the random and the structured aspects of the spatial dispersion. The variogram and the co-variogram have been widely used to describe the spatial structure of ecological variables (Legendre and Fortin, 1989; Rossi et al., 1992). The traditional estimator of the co-variogram is defined as (Isaaks and Srivastava, 1989) where Z(xi), Z(xi + h) and T(xi), T(xi + h) are measurements of two variables at locations xi and xi + h separated by the vector of directional distance h, and N(h) is the number of pairs of samples considered in the given distance class. This calculation is repeated for different values of h and provides the empirical co-variogram, which is a plot of the values of c(h) as a func- tion of discrete distance h describing the joint spatial cross-covariance of two variables. When only one variable is analysed (i.e. Z and T are the same) c(h) is called empirical variogram. To provide a smooth, continuous description of the covariance and cross-covariance spatial structure of the parameters, variogram models were derived by applying admissible functions (Wackernagel, 1995) to empirical estimates cal- culated from the data sets. Basic variogram models were added together in a nested struc- ture and fitted to experimental estimates by an ordinary least square procedure (Pebesma and Wesseling, 1998). Modeled variograms and co- variograms were used in ordinary block kriging and cokriging techniques for spatial interpola- tion to produce decomposition maps (Isaaks and Srivastava, 1989; Rossi et al., 1992). Kriging is a weighted, moving-average interpolation method where the set of weights to estimate val- ues at unmeasured points is computed as a func- tion of the variogram model and locations of the ( ) ( ) ( ( ) ( )) ( ( ) ( )) N Z Z T Tx x x x h h h h 2 1 i i i i i $ $ ) = - + - + c / Fig. 1. NDVI values for MIVIS and Landsat images at different spatial resolutions (lines) and leaf half-time life (bars) in the four vegetation types. Data are mean values from five sampling sites per vegetation type. 222 Letizia Sabetta et al. samples. Kriging is a BLUE estimator (Best Linear Unbiased Estimator, Isaaks and Srivasta- va, 1989) attempting to minimize the estimation variance of predicted values. Estimation can be done for single points or blocks (i.e. block krig- ing) and exploiting the spatial properties of just one variable or the joint variation of two (i.e. co- kriging). For a detailed description of the proce- dure refer to Isaaks and Srivastava (1989) or Wackernagel (1995). 3. Results Litter mass loss showed a significant fit with the negative exponential model (Olson, 1963) in all study sites and for each plant litter type. Senescing leaves of hazelnut (Corylus avellana) had shorter average half-life (199±22 days) showing faster decomposability than, in order, chestnut (Castanea sativa, 356 ± 28 days), oak (Quercus cerris, 237±17 days) and beech (Fagus sylvatica, 461±37 days; table I; fig. 1); differ- ences in the litter decomposition rate among veg- etation types were significant (nested ANCOVA, table II(a)). We also observed significant differ- ences of the decomposition rate among study sites within each vegetation type (nested ANCO- VA, table II(a)). NDVI value differed both among vegetation types and among spatial scales; hazelnut area had lower average NDVI value at each considered spatial scale, showing lower vegetation cover, than chestnut, oak and beech areas (table I; fig. 1). On an average NDVI was lower from Landsat than from MIVIS im- ages, but the differences observed among plant types were scale-invariant (two-way ANOVA, table II(b)). In the studied area litter decomposi- tion rate (half-life) and vegetation cover (NDVI) were positively related to all considered spatial resolutions (table III). Slopes of this observed re- lationships did not differ from each other (AN- COVA, variable: litter half-life time, covariable: NDVI); on the contrary, intercepts of the linear relationships differed significantly between Landsat and MIVIS (table IV), the former being lower than the latter (903 versus 1178 ± 85). Significant spatial autocorrelations of the half-life and NDVI value were found within the study area (table V). A significant amount of the observed variation in half-life and NDVI was spatially structured; semivariograms showed that Table I. Litter half-life time and NDVI for MIVIS (three flight heights) and Landsat images in the four vege- tation types. Mean values from 5 sampling points and standard errors (in brackets) are reported. Type Half-life (d) MIVIS NDVI Landsat NDVI 5000 m 4000 m 2000 m Corylus avellana 199 (22) 0.567 (0.013) 0.593 (0.021) 0.615 (0.015) 0.479 (0.062) Castanea sativa 356 (28) 0.683 (0.042) 0.697 (0.049) 0.707 (0.051) 0.619 (0.093) Quercus cerris 237 (17) 0.673 (0.012) 0.701 (0.013) 0.725 (0.009) 0.586 (0.013) Fagus sylvatica 461 (37) 0.711 (0.014) 0.730 (0.014) 0.754 (0.011) 0.594 (0.030) 223 Mapping litter decomposition by remote-detected indicators Table II. (a) Nested ANCOVA analysis of remaining litter mass as variable and time as covariable (sampling sites nested in vegetation type as main factor); (b) two-way ANOVA analysis of NDVI as variate (factor A = spa- tial resolution, factor B = vegetation type) followed by Tukey HSD test (main factor: vegetation type on the left and spatial resolution on the right). Source DF MS effect F p-level (a) Nested ANCOVA: litter mass (covariate: time) Between vegetation types 3 0.34 8.699 1.22∗10−5 Within vegetation types 16 0.076 1.943 0.015 (b) Two-way ANOVA: NDVI Spatial resolution (A) 3 0.066 22.28 5.4∗10−10 Vegetation type (B) 3 0.073 24.37 1.2∗10−10 Interaction A×B 9 0.001 0.3 0.97 Tukey HSD test Corylus avellana Quercus cerris Fagus sylvatica Landsat 2000 m 4000 m Quercus cerris 1.5∗10−4 2000 m 1.5∗10−4 Fagus sylvatica 1.5∗10−4 0.44 4000 m 1.5∗10−4 0.64 Castanea sativa 1.5∗10−4 0.99 0.63 5000 m 1.7∗10−4 1.5∗10−4 0.59 Table III. Regression and cross-variogram statistics of half-life time (y) and NDVI (x) from Landsat ETM+ and MIVIS images. Regression Cross-variogram Function r p Model r p Landsat y = 903x−180 0.53 < 0.05 Linear 0.98 < 0.001 5000 m y = 1256x−492 0.72 < 0.001 Linear 0.95 < 0.001 MIVIS 4000 m y = 1087x−407 0.65 < 0.01 - - - 2000 m y = 1190x−500 0.71 < 0.001 - - - Table IV. Summary statistics for estimated marginal means (a) and p-values for pairwise comparisons follow- ing Tukey HSD test (b) from ANCOVA analysis of half-life time as variate and NDVI as covariate after test for parallelism. (a) Spatial scale Mean (1) Standard error 95% Confidence interval Lower bound Upper bound Landsat 397.93 26.87 344.34 451.51 2000 m 270.11 24.39 221.47 318.75 4000 m 289.86 23.48 243.04 336.68 5000 m 311.09 23.02 265.18 357.01 (b) Landsat 2000 m 4000 m 2000 m 0.002 4000 m 0.005 0.548 5000 m 0.018 0.222 0.519 (1) Adjusted means for NDVI equal to 0.646. Fig. 2a,b. Spatial variation of half-life time in 60 × 60 m-support block size. a) Ordinary block kriging map. b) Ordinary block co-kriging map based on NDVI from MIVIS data (5000 m). Open circles locate the sampling sites (Corylus avellana, no; Quercus cerris, ce; Castanea sativa, ca; Fagus sylvatica, fa). 224 Letizia Sabetta et al. the proportion of sample variance (C0 + C) ac- counted for by spatially structured variance (C) was on average 81 ± 12%. The semivariance of NDVI increased regu- larly with the separation distance, up to 1541 m for Landsat data and 2364 m for MIVIS data, following spherical model, and showed negligi- ble nugget variance (table V). The decomposi- tion rate had an increase of semivariance up to 1627 m and a high nugget variance. Table V. Variogram model parameters of NDVI from the Landsat ETM+ and MIVIS (5000 m) images and the litter half-life time. Variable Model Semivariogram parameters Proportion r 2 RSS Nugget (C0) Sill (C0+C) Range (m) C/(C0+C) Landsat NDVI Spherical 0.006 0.02 1541 0.700 0.993 9.2∗10−7 MIVIS NDVI Spherical 0.002 0.01 2364 0.800 0.986 1.3∗10−6 Half-life time Linear with sill 1230 17870 1627 0.931 0.895 2.05∗10+7 a b 225 Mapping litter decomposition by remote-detected indicators Litter half-life and NDVI were also spatial- ly co-structured; the cross-variogram analysis showed inter-relationships among the two auto- correlated variates with positive co-regionaliza- tion (table III). The half-life increased north- wards (fig. 2) from the hazelnut area to the beech forests. Cokriging maps from Landsat and MIVIS data did not differ. 4. Discussion and conclusions Some major points can be drawn from the results: – Decomposition rate is both plant species-dependent and spatially dependent. In fact it varies with the vegetation type and with- in each forest type; it is significantly heteroge- neous because of the heterogeneous decompos- er (microorganisms and detritivores) distribu- tion. NDVI, a remote-sensed index reflecting structural (i.e. cover and density of the leaf lay- er) and functional (i.e. physiological status of the leaf layer) properties of vegetation, is pecu- liar to the four CORINE habitats under study and satisfactorily relates with the decomposi- tion rates, measured in situ. – The information content of NDVI varies at each resolution scale. Differences between MIVIS and Landsat data are related mainly to spatial and spectral resolution of the sensors, even if a residual atmospheric interference not accounted for DOS approach has to be suggest- ed. Further investigation is needed to evaluate their relative importance. However, the infor- mation decay from MIVIS to Landsat images does not modify the relationship existing with the litter decomposition rate. – Using NDVI as co-variable, the co-krig- ing map of litter decomposition is better shaped than the kriging map. Co-kriging offers addi- tional advantages over kriging as it involves co- variate that is cross-correlated with the variable of interest and that can be more easily sampled (Isaaks and Srivastava, 1989). NDVI, as spatial- ly explicit co-variable, can thus sustains the mapping algorithm where decomposition infor- mation is scarce or missing. In conclusion, the results show that remote sensing techniques can provide outstanding help in developing the cartography of ecosystem functions, also when availability of the ground- truth data is low. Landsat data are as useful as MIVIS data for the remote sensing analysis and multiple scales mapping of decomposition. This is a clear advantage since the Landsat data can be more cheaply and regularly acquired than the MIVIS data. Acknowledgements This research was founded by the joint pro- gram «Progetto n. 2: Conoscenza, conservazio- ne e gestione della biodiversità» of the Italian Ministry of the Environment and the Italian Na- tional Research Council. We would like to ac- knowledge the contributions made by the two anonymous reviewers who helped with valu- able comments. REFERENCES ABER, J.D. and J.M. MELILLO (1980): Litter decomposition: measuring relative contribution of organic matter and nitrogen to forest soil, Can. J. Bot., 58, 416-421. ABER, J.D., J.M. MELILLO and C. MCCLAUGHERTY (1990): Predicting long-term patterns of mass loss, nitrogen dynamics, and soil organic matter formation from ini- tial fine litter chemistry in temperate forest ecosystems. Can. J. Bot., 68, 2201-260. AERTS, R. (1997): Climate, leaf litter chemistry and leaf lit- ter decomposition in terrestrial ecosystems: a triangu- lar relationship, Oikos, 79, 439-449. CARD, D.H., D.J. PETERSON, P.A. MATSON and J.D. ABER (1988): Prediction of leaf chemistry by the use of visi- ble and near-infrared reflectance spectroscopy, Remote Sensing Environ., 26, 123-147. CHAVEZ, P.S. JR. (1988): An improved dark-object subtraction technique for atmospheric scattering correction of multi- spectral data, Remote Sensing Environ., 24, 459-479. CURRAN, P.J. (2001): Remote sensing: Using the spatial do- main, Environ. Ecol. Stat., 8, 331-344. DUGAN, J.P., D.L. PETERSON and P.J. CURRAN (1994): Alter- native approaches for mapping vegetation quantities using ground and image data, in Environmental Infor- mation Management and Analysis: Ecosystem to Glob- al Scales, edited by W. MICHENER, J. BRUNT and S. STAFFORD (Taylor and Francis, London), 237-261. EUROPEAN UNION/DIRECTORATE GENERAL XI (1991): CO- RINE Biotopes Manual, Habitats of the European Community. A Method to Identify and Describe Consis- tently Sites of Major Importance for Nature Conserva- tion (Bruxelles), EUR 12587/3. GALLARDO, A. and J. MERINO (1993): Leaf decomposition 226 Letizia Sabetta et al. in two Mediterranean ecosystems of Southwest Spain – Influence of substrate quality, Ecology, 74 (1), 152- 161. GILLON, D., R. JOFFRE and A. IBRAHIMA (1999): Can litter decomposability be predicted by near infrared re- flectance spectroscopy?, Ecology, 80 (1), 175-186. ISAAKS, E.H. and R.M. SRIVASTAVA (1989): Applied Geosta- tistics (Oxford University Press, Oxford). KRUSE, F.A., K.S. KIEREN-YOUNG and J.W BOARDMAN (1990): Mineral Mapping at Cuprite, Nevada with a 63 channel imaging spectrometer, Photogramm. Eng. Re- mote Sensing Environ., 56 (1), 83-92. LEGENDRE, P. and M.-J. FORTIN (1989): Spatial pattern and ecological analysis, Vegetatio, 80, 107-138. MARTIN, M.E. and J.D. ABER (1994): Analysis of forest fo- liage, III. Determining nitrogen, lignin and cellulose in fresh leaves using near infrared reflectance data, J. Near Infrared Spectrosc., 2, 25-32. OLSON, J.S. (1963): Energy storage and the balance of pro- ducers and decomposers in ecological systems, Ecolo- gy, 44, 322-331. PEBESMA, E.J. and C.G. WESSELING (1998): GSTAT: a pro- gram for geostatistical modelling, prediction and simu- lation, Comput. Geosci., 24 (1), 17-31. PEREZ-HARGUINDEGUY, N., S. DIAZ, J.H.C. CORNELISSEN, F. VERDRAMINI, M. CABIDO and A. CASTELLANOS (2000): Chemistry and toughness predict leaf litter decomposi- tion rates over a wide spectrum of functional types in Central Argentina, Plant Soil, 218, 212-30. RICHARDS, J.A. (1994): Remote Sensing Digital Image Ana- lysis (Springer-Verlag, Berlin), p. 340. ROSSI, R.E., D.J. MULLA, A.G. JOURNEL and E.H. FRANZ (1992): Geostatistical tools for modeling and interpret- ing ecological spatial dependence, Ecol. Monogr., 62, 277-314. ROUSE, J.W., R.H. HAAS, J.A. SHELL, D.W. DEERING and J.C. HARLAN (1974): Monitoring the vernal advance- ment of retrogradation of natural vegetation, NASA/ GSFC, Greenbelt, MD, Final Rep. SAUNDERS, G.W. (1976): Decomposition in freshwater, in The Role of Terrestrial and Aquatic Organism in De- composition Processes, edited by J.M. ANDERSON and A. MACFADYEN (Blackwell Scientific), 341-374. SWIFT, M.J., O.W. HEAL and J.M. ANDERSON (1979): De- composition in Terrestrial Ecosystems (Univ. Calif. Press, Berkeley), pp. 509. USGS (2004): Revised Landsat 5 TM Radiometric Calibra- tion Procedures and Post-Calibration Dynamic Ranges, edited by J. CHANDER and B. MARKHAM (USGS, Landsat Project; on line http://landsat7.usgs.gov/resource.html). WACKERNAGEL, H. (1995): Multivariate Geostatistics. An In- troduction with Applications (Springer-Verlag, Berlin).