journa l homepage: www.e lsev ier .com/ locate /eco l ind
ombining geostatistical models and remotely sensed data to improve tropicalree richness mapping
. Luis Hernández-Stefanonia,∗, J. Alberto Gallardo-Cruzb, Jorge A. Meaveb, Juan Manuel Dupuya
Centro de Investigación Científica de Yucatán A.C., Unidad de Recursos Naturales, Calle 43 # 130, Colonia Chuburná de Hidalgo, C.P. 97200, Mérida, Yucatán, MexicoDepartamento de Ecología y Recursos Naturales, Facultad de Ciencias, Universidad Nacional Autónoma de México, México 04510, D.F., Mexico
r t i c l e i n f o
rticle history:eceived 18 December 2009eceived in revised form 8 November 2010ccepted 11 November 2010
eywords:o-kriging
mage textureriging with external driftegression krigingree richnessropical forest
a b s t r a c t
Information on the spatial distribution and composition of biological communities is essential in design-ing effective strategies for biodiversity conservation and management. Reliable maps of species richnessacross the landscape can be useful tools for these purposes. Acquiring such information through tradi-tional survey techniques is costly and logistically difficult. The kriging interpolation method has beenwidely used as an alternative to predict spatial distributions of species richness, as long as the data arespatially dependent. However, even when this requirement is met, researchers often have few sampledsites in relation to the area to be mapped. Remote sensing provides an inexpensive means to derive com-plete spatial coverage for large areas and can be extremely useful for estimating biodiversity. The aim ofthis study was to combine remotely sensed data with kriging estimates (hybrid procedures) to evaluatethe possibility of improving the accuracy of tree species richness maps. We did this through the compari-son of the predictive performance of three hybrid geostatistical procedures, based on tree species densityrecorded in 141 sampling quadrats: co-kriging (COK), kriging with external drift (KED), and regressionkriging (RK). Reflectance values of spectral bands, computed NDVI and texture measurements of Landsat7 TM imagery were used as ancillary variables in all methods. The R2 values of the models increased from0.35 for ordinary kriging to 0.41 for COK, and from 0.39 for simple regression estimates to 0.52 and 0.53when using simple KED and RK, respectively. The R2 values of the models also increased from 0.60 for
multiple regression estimates to 0.62 and 0.66 when using multiple KED and RK, respectively. Overall, ourresults demonstrate that these procedures are capable of greatly improving estimation accuracy, withmultivariate RK being clearly superior, because it produces the most accurate predictions, and becauseof its flexibility in modeling multivariate relationships between tree richness and remotely sensed data.We conclude that this is a valuable tool for guiding future efforts aimed at conservation and management
fores
of highly diverse tropical
. Introduction
Comprehensive information on the spatial distribution andomposition of plant species richness is fundamental to designffective strategies for conservation and management of increas-ngly fragmented tropical forests (Rodríguez et al., 2007; Gillespiet al., 2008). Unfortunately, such information cannot be obtainedxclusively from traditional survey techniques, due to logistical dif-culties and the costs involved. An alternative for predicting thepatial distribution of plant species richness is through the usef kriging, a spatial interpolation procedure, as long as the data
re spatially dependent (Isaaks and Srivastava, 1989; Webster andliver, 2001). Even if field samples for estimating species rich-ess are spatially autocorrelated, they are often sparse, which can
lead to considerable uncertainty in the ordinary kriging estimations(Baxter and Oliver, 2005). Several approaches have been proposedto enhance the accuracy of estimations of spatial distribution ofspecies richness. One of the most common approaches is to assessspecies diversity based on the average values of species richnesswithin vegetation or land cover classes. In this method, the mappedclasses are viewed as habitats and diversity within those classesis assessed through field samples (Nagendra and Gadgil, 1999;Wagner et al., 2000). However, a major drawback of this approach,aside from the simplification of having a single mean value predict-ing all non-measured points within each class, is that it assumessample independence, i.e. it does not account for spatial depen-dence or auto-correlation (Dormann et al., 2007).
Biodiversity estimates using remotely sensed data have
achieved some success based on indicators of both environmen-tal heterogeneity and net primary productivity (Gillespie, 2006;Rocchini, 2007). In this sense, hybrid geostatistical procedures, acombination of kriging with ancillary variables, have been sug-
ested to enhance the accuracy of estimation (Odeh et al., 1994;ebster and Oliver, 2001). Researchers have used ancillary vari-
bles that are correlated with the primary variable to improve thestimation process. The ancillary variable is more intensely sam-led than the original variable because it is easier or cheaper toeasure, reducing the cost of sampling and increasing the accu-
acy of prediction of the sparsely sampled target variable. Remoteensing offers an inexpensive means of deriving complete spatialoverage of environmental information for large areas, at regularime intervals, and can therefore be extremely useful in providingncillary variables to estimate species richness (Nagendra, 2001).
There are several methods that use intensive ancillary data toncrease the accuracy of predictions of the variable of interest.mong these procedures are co-kriging (COK), kriging with exter-al drift (KED) and regression kriging (RK) (Goovaerts, 1999). Theseethods differ in the way they utilize the ancillary information to
stimate the target variable at unsampled locations. In COK, forxample, both variables are spatially correlated and have a cross-ariogram that is used to quantify cross-spatial auto-covarianceetween the primary and secondary variable, which in turn is used
n the interpolation procedure (Webster and Oliver, 2001). KED usesvariogram of linear regression residuals of a local trend between
he primary and secondary variables, implicitly estimated throughhe kriging system, to guide the interpolation method (Goovaerts,999). RK combines a simple or a multiple regression model withrdinary kriging of regression residuals to estimate the primaryariable (Odeh et al., 1994; Webster and Oliver, 2001).
Several studies have tried to combine ancillary data with krigingstimates to examine whether it is possible to improve the accu-acy of predictions of different ecological variables. COK has beensed successfully, for instance, for mapping herbaceous biomass
n a tropical savanna (Mutanga and Rugege, 2006) or abundanceistribution of small pelagic species (Georgakarakos and Kitsiou,008) using remotely sensed data as ancillary variables. Similarly,ales et al. (2007) used elevation, vegetation type and soil textureo improve the spatial prediction of forest biomass with KED. Inddition, Hernández-Stefanoni et al. (2009) assessed the spatial dis-ribution of species density and abundance of tropical trees throughhe use of RK and remotely sensed data. However, most studies thatombine geostatistical models with ancillary information to predictcological variables have focused on a single procedure using one orew ancillary variables. Therefore, evaluating the accuracy of differ-nt procedures using single or multivariate secondary informationould improve the spatial mapping of relevant ecological informa-ion, which in turn could provide ecological insights and inform the
anagement and conservation of natural resources.Numerous remotely detectable parameters, such as environ-
ental factors, influence the distribution of species, and thus maye used as ancillary variables. Net primary productivity (NPP) isne factor which has been related with the normalized differenceegetation index (NDVI) and correlated with species richness ategional and local scales (Oindo and Skidmore, 2002; Fairbanks andcGwire, 2004). This is based on the idea that areas with more
esources tend to support a greater number of species (Braakhekkend Hooftman, 1999; Gould, 2000), at least at certain spatial scales.nother factor having a potential predicting power for speciesiversity at local scales is heterogeneity of land cover types ornvironmental heterogeneity measured as heterogeneity of spec-ral indices or spectral variability derived from satellite imageryRocchini, 2007; Oldeland et al., 2010; Rocchini et al., 2010). Heree ascribe to the general hypothesis that greater heterogene-
ty allows a higher number of species to coexist (Huggett, 1995;
ilson, 2000). It is necessary to evaluate which variables derived
rom satellite imagery provide the most accurate estimates of plantpecies richness at unsampled locations using hybrid geostatisticalechniques.
l Indicators 11 (2011) 1046–1056 1047
The goals of this study were twofold: (1) to quantify and com-pare the improvement in the accuracy of tree species richness mapswith different geostatistical methods, resulting from the use of eco-logical indicators derived from spectral data and indices obtainedfrom a remote sensor and (2) to make general suggestions for suit-able mapping methods and the use of remotely sensed data assources of secondary information. We use one or more types ofspectral data and derive indices from a remote sensor, as ancil-lary information in COK, KED and RK, to compare the performanceof these hybrid techniques in predicting the spatial distribution oftree species richness of a tropical forest.
2. Data and methods
2.1. Study area
The study was conducted in a tropical landscape locatedin the Yucatan Peninsula, Mexico, over an area of 64 km2
(18◦53′54′′–18◦58′14′′ N latitude and 88◦10′04′′–88◦14′37′′ W lon-gitude). The terrain is fairly flat and the climate is tropical,sub-humid, with a distinct dry period. Mean annual temperatureis 26 ◦C, and average total annual rainfall ranges between 1000 and1300 mm, concentrated between June and October, with the dryseason extending from December to April (Orellana et al., 2009).
The plant cover mosaic in this landscape comprises three vege-tation types: semi-evergreen forest, savanna and lowland floodedforest. The semi-evergreen forest, dominant in the landscape,grows up to 25 m tall and is a structurally complex community withtwo or three canopy layers consisting mostly of trees, shrubs andvines; for this community four successional stages of different agesof abandonment after traditional slash-and-burn agriculture weredistinguished: Stage 1 (secondary vegetation ≤ 3 yr old); Stage 2(4–10 yr old); Stage 3 (11–19 yr old), and Stage 4 (forests > 20 yr old,including old-growth stands). The other two vegetation types areassociated to seasonally flooded areas where agriculture is rarelypracticed. Savanna is a herb-dominated community, as it lacks acontinuous tree cover, in contrast to the lowland flooded forest,whose tree canopy may be as tall as 10 m (Cabrera et al., 1982).Thus we distinguished in total six vegetation cover types in thestudy.
2.2. Species richness data
Field data were recorded during two plant surveys conductedduring the rainy seasons of 2000 and 2001. The surveys were basedon a stratified random sampling design in which each stratum (veg-etation type) was initially randomly sampled proportionally to itsarea. The number of sites per stratum was subsequently modified toprovide a more adequate representation of species richness; both atthe stratum and the landscape level, based on species accumulationcurves (see Hernández-Stefanoni and Ponce-Hernández (2004) fordetails). The final sample size was 141 quadrats, with sample dis-tribution among strata as follows: 20 sites in successional Stage1 of the semi-evergreen forest, 20 in Stage 2, 25 in Stage 3, 42in Stage 4, 17 in low flooded forest, and 17 in savanna. All sam-pled quadrats were located on the ground with a GPS unit. Dueto a better taxonomic knowledge and reliability of identification,the sampling was restricted to tree species. In determining spe-cific identities of plants, the aid of local inhabitants was essential(Hernández-Stefanoni et al., 2006).
Since most tropical plant species have low local densities, the
assessment of species richness at the community level is difficultand requires sampling large areas in order to capture the rarityof most species in tropical forests. This can be achieved either byusing few very large sample units, or by placing a large number of
mall samples over the area. We used small sample plots (100 m2)ot only because of their practical value, but also because theirse allows a suitable estimation of plant diversity at the landscapecale (Clark et al., 1999; Plotkin et al., 2000), especially in patchynd fragmented landscapes. Moreover, this is the scale at whichnteractions between individual plants occur without substantialariation in geomorphology, hydrology, soil and climate variablesPalmer et al., 2000).
Each sampling site consisted of a 10 m × 10 m quadrat used toample trees taller than 3 m. In each quadrat we computed theumber of tree species per 100 m2 sample site (i.e. species densityensu Gotelli and Colwell, 2001) as a measure of species richness.
.3. Remotely sensed data and imagery processing
A Landsat 7 Thematic Mapper (TM) satellite image (30 m × 30 m)rom April 2000 was geo-referenced to Universal Transverse Mer-ator Projection (WGS 84) and reduced to the study area. Themage was radiometrically and atmospherically corrected to min-mize the effect of atmospheric scattering (Chavez, 1996). TheN (Digital Numbers) values of TM bands were converted to
adiance and transformed to exo-atmospheric reflectance units,sing the equations suggested by the Landsat 7 Handbook (NASA,998). The pixels containing each of the 141 quadrats werehen identified and the reflectance in three wavebands (5—short-ave infrared: 1.55–1.75 �m, 4—near infrared: 0.76–0.90 �m and
—red: 0.63–0.69 �m) were extracted. These bands were selectedecause red, near infrared and middle infrared bands have been rec-mmended for species discrimination (Nagendra, 2001), althoughifferent plant species have differential responses to light in thelectromagnetic spectrum. Moreover, Hernández-Stefanoni andupuy (2007) found in the same study area that the strongest asso-iation between plant species richness and spectral reflectance wasbserved when these bands were used. The normalized differenceegetation index (NDVI) was calculated from spectral reflectanceands near infrared (NIR, band 4) and red (R, band 3) with theormula NDVI = (NIR − R)/(NIR + R), and extracted.
To perform a supervised classification, a false RGB color com-osite image was created from bands 5, 4 and 3 of the Landsat
mage. This composite image was used as a spatial reference frame-ork for selecting suitable training sites for the different land cover
ypes to be identified on the ground. At least two training sites wereelected for each land cover type. Using a GPS receiver, the selectedraining sites were located in the field. In addition to the six veg-tation cover classes, the land cover types also included pasturesnd cropping areas which were excluded from all analyses. Waterodies were also identified in the false color composite image. Theraining areas and the false color image were used to perform atandard supervised classification of the Landsat 7 TM bands, usinghe Maximum Likelihood Algorithm as implemented by ER Mapper.1 (Earth Resource Mapping Ltd, 1998). To preserve the local vari-bility in the classified image, the values of pixels were not changedy median or majority filters, allowing the speckle in the classified
mage. To assess map accuracy, two measures were used: over-ll accuracy and Cohen’s Kappa statistic (Campbell, 1987; Jensen,000; see Hernández-Stefanoni and Ponce-Hernández (2004) foretails).
Heterogeneity of land cover types was also examined as aotential predictor of species diversity, assuming that a greaterabitat heterogeneity allows a higher number of species to coex-
st (Nagendra, 2001; Hernández-Stefanoni and Dupuy, 2008). To
ssess the spatial heterogeneity at the local scale, we used a mov-ng window analysis to calculate two image texture measurementsrom the classified land cover map; window size was 3 × 3 pixels,hich is the smallest window capable of accounting for the local
l Indicators 11 (2011) 1046–1056
variability of the neighborhood area. Texture measurements wererelative richness and dominance indices, which are landscape met-rics that quantify the composition of habitats for a specific group ofpixels (Turner, 1989). The pixels containing each of the 141 sam-pling quadrats were then identified and the values of the two imagetexture measurements were extracted.
Finally, a Principal Component Analysis (PCA) was performedusing reflectance values of bands 5, 4 and 3 of Landsat 7 TM,NDVI values and two texture measures. This was done to developspectral-based indices (new variables) that integrate in a singlevalue different variables that are independent of each other (i.e.orthogonal).
2.4. Geostatistical prediction methods
2.4.1. Co-krigingCo-kriging (COK) is an extension of ordinary kriging (OK) that
takes into account the spatial cross-correlation between two vari-ables, i.e. the primary variable (tree species richness in our case),and another ancillary variable (in this study, derived from remotelysensed data). The primary variable has been sampled at few placeswhile the secondary one has a denser set of data. Both variableshave an auto-variogram that provides information on the structureof the spatial variability which helps define the size and shape of theneighborhood for interpolation. The semi-variogram is computedfrom:
�(h) = 12n
n∑
i=1
(Z(xi) − Z(xi + h))2
where Z(xi) is the primary variable (tree species richness or a vari-able derived from remote sensing data), sampled in quadrat i;Z(xi + h) is the value of variable sampled in other quadrats sepa-rated from xi, by a discrete distance h; n represents the numberof pairs of observations separated by h, and �(h) is the estimated“experimental” semi-variance value for all pairs at a lag distance h.
If the two variables are found to have a cross-variogram, this isthen used to quantify cross-spatial auto-covariance between them(Isaaks and Srivastava, 1989; Webster and Oliver, 2001). The cross-semivariance is computed through the equation:
�12(h) = 12n
n∑
i=1
[Z1(xi + h) − Z1(xi)][Z2(xi + h) − Z2(xi)]
where Z1 is a tree plant species richness value; Z2 is a given ancil-lary variable; �12(h) represents the cross-semivariance and h isthe distance between both locations x + h and x. The three semi-variogram models were fitted to the experimental semi-variogramsusing gstat for R software (Pebesma, 2004). The coefficient of deter-mination (R2) of the least-squares fit was used both as the criterionto select the best model and to fit the three experimental semi-variograms.
The COK predictions of tree richness are obtained from the fol-lowing equation:
Z(x0) =n∑
i=1
�iZ(xi)
where �i are the optimal weights selected to minimize the esti-mation variance (Burrough and McDonnell, 1998); Z(xi) are theobserved values of tree species richness, and Z(x0) is the optimaland unbiased estimate of tree species richness. The estimated maps
were obtained using at least 16 sampling quadrats for each un-estimated location within a maximum radius of 8 km.
To conduct the COK estimation, several experimental cross-variograms were analyzed, using plant species richness in 141
uadrats as primary variable and three correlated ancillary vari-bles derived from remotely sensed data assessed in 850 randomlyelected points. The ancillary variables were the reflectance valuef TM4 band, the NDVI index, and the extracted spectral-basedrincipal component (PC1), because they were the most stronglyorrelated to the main variable (tree species richness).
.4.2. Kriging with external driftKriging with external drift (KED) uses one or more secondary
ariables to assess the local trend m(x) that varies smoothly withinach local neighborhood (Goovaerts, 1999). It assumes a linear rela-ionship between primary and secondary variables and is obtainedrom:
(x) = ˇ0 +n∑
i=i
ˇiyi
here ˇ0 and ˇi′s are the intercept and coefficients associated withhe variables of a linear regression model. The unknown coeffi-ients are regarded as constant within the search neighborhood andmplicitly estimated through the kriging system. The KED estimators:
KED(x0) =n∑
i=1
�KEDi Z(xi)
he KED weights �KEDi
are determined to include the effect of theecondary variables and are obtained through a solution of a systemf equations (Webster and Oliver, 2001), in which the semivari-nces of residuals between the points to be predicted and dataoints are included. The value of the ancillary variables must benown at all primary sampling locations and at all locations beingstimated. Moreover, a regression model is assessed locally in KEDithin the search neighborhood.
In this study we used two models to estimate the spatial distri-ution of tree species richness. The first one (simple KED) used aimple linear regression between species richness and PC1, whilehe second model (multiple KED) used species richness as depen-ent variable and reflectance values of the three spectral bands (3,, and 5) of Landsat, computed NDVI values, and two texture mea-ures (relative richness and dominance) as explanatory variables inmultiple regression model.
.4.3. Regression krigingThis approach combines a simple or multiple regression model
ith OK of regression residuals, i.e. this method addresses both thepatial dependence of observations and the relationship betweenhe dependent variable (tree species richness) and the ancillaryariables (spectral features of a Landsat image) (Odeh et al., 1994;ebster and Oliver, 2001). The regression kriging (RK) estimator of
ree species richness Zrk(x) is defined as the sum of regression esti-ate Zr(x) obtained as a linear function between tree richness and
oth the spectral features of the Landsat 7 image and the krigedstimate of spatially correlated residual values εok(x), using theollowing equation:
rk(x) = Zr(x) + εok(x)
o obtain regression kriging estimates, two linear regression mod-ls were used to predict tree species richness. The first modelas a simple linear regression between plant species richness and
C1. The second model, a multiple linear regression model, used
eflectance values of the three spectral bands (3, 4, and 5) of Landsatimage, computed NDVI values, and two texture measures (relative
ichness and dominance) as explanatory variables. The explana-ory variables were transformed as needed with 1/x, log10 (x), log10
Fig. 1. Land cover map of the study area obtained from a supervised classification.No median or majority filters were applied.
(x + 1) and sqrt (x) to meet the linearity assumptions (Tabachnickand Fidell, 1996).
In total, 15 models considering different combinations ofexplanatory variables were evaluated for the multiple regressionmodel (Table 3). The Akaike Information Criterion (AIC) was usedto select the best model. This procedure is based on parsimony,a trade-off between model fit and the number of parameters inthe model (Anderson and Burnham, 2002). AIC values were calcu-lated from the formula AIC = −2*(log likelihood) + 2K, where K is thenumber of parameters. The models were ranked based on delta AICvalues (�i) (Anderson et al., 2000; Johnson and Omland, 2004). Aset of candidate models were selected using both an approximatecutoff of �i = 4, and those models having wi > 0.1 (Burnham andAnderson, 1998). Finally, we calculated model-averaged parame-ters and unconditional standard errors based on the Akaike weights(Burnham and Anderson, 1998). Residuals were computed by sub-tracting observed values of tree species richness from predictedvalues using each regression model. Estimates of residual valuesat unobserved locations were obtained from the OK. The semi-variogram model and OK residuals were obtained using gstat forR software (Pebesma, 2004).
2.5. Evaluation of geostatistical procedures
The performance of the different interpolation methods wasassessed by leave-one-out cross-validation (Isaaks and Srivastava,1989). In this procedure, one observation is temporally removedfrom the data set, and the remaining sampling values are usedto predict the variable. This cross validation yields a list of esti-mated values of tree species richness paired to those obtained fromthe observed sampling units. The different models were comparedbased on the coefficient of determination (R2) obtained from asimple regression between observed and estimated values of treespecies richness and the root mean square error (RMSE). We alsoevaluated whether patterns of tree species richness among vegeta-tion cover types, produced by the different geostatistical methodsreflected those patterns obtained in the field. This was achievedthrough the use of analysis of variance of species richness estimatesproduced by different methods for the sampled quadrats.
3. Results
3.1. Land cover map
Fig. 1 shows the land cover thematic map of the study area pro-duced by the supervised classification. The total area occupied by
Table 1Total variance explained by two principal components derived from Landsat spectralmeasures and the weights of the variables in each component after rotation.
igures marked in bold typeface are loadings >0.70.
his landscape was 6464.2 ha; the six vegetation cover types jointlyovered 82.6% of this area, whereas pastures, cropping areas andater bodies accounted for the remaining 17.4%. Forest in succes-
ional Stage 4 was the largest land cover type spanning over 36.6%2368.53 ha) of the total area. Younger successional stages repre-ented increasingly smaller proportions of the area: Stage 3, 25.7%1661.40 ha) and Stage 2, 11.3% (728.64 ha). The remaining coverypes represented only 8.0% of the area. Overall accuracy of the
ap was 70.8%, while Cohen’s Kappa was 0.63. The most conspic-ous classification errors were those made between successionaltages, probably due to the fuzzy boundaries between them.
.2. Tree species richness of vegetation classes
A total of 14,417 individuals belonging to 118 plant speciesere recorded in the 141 sampling quadrats. The mean values of
ree species richness by quadrat differed significantly among theix vegetation classes (F[5,140] = 53.81, P < 0.001). Least significantifference tests revealed that the six vegetation classes differedignificantly in the number of plant species (P < 0.001; see columnor observed values in Table 6). The values showed the follow-ng pattern: overall, semi-evergreen forest stages hosted more treepecies than the other two vegetation cover types, and within semi-vergreen forests, species richness increased from the youngest tohe oldest successional stage.
.3. Co-kriging
To conduct the COK estimation, several experimental cross-ariograms were analyzed. In the PCA analysis, the first twoomponents accounted for 58.33% (PC1) and 17.43% (PC2) of theariation after varimax rotation. PC1 was interpreted as a gradientf spectral values of remotely sensed data including the three spec-ral bands (3, 4, and 5) of Landsat 7 image and NDVI. In turn, PC2as associated primarily with texture measurements (Table 1).
The strongest spatial cross-variogram structures were foundor tree species richness and PC1. Table 2 shows the best fitteduto-variograms and cross-variogram models, generated for COK.he least-square fitting technique with minimum reduced sum ofquares (RSS) and maximum R2 values were used to select the bestodels. The structural variance, which determines the variance
ue to spatial dependence explained by the model, ranged from4.4% to 76.7%, suggesting that there is more unexplained vari-bility (nugget variance) of PC1 than of tree species richness andree species richness/PC1 values. It also indicates that PC1 may vary
ore over small scales.
The range of influence was between 422 m and 2839 m. This
uggests that one would reasonably expect tree richness and PC1alues to be still correlated between places separated by such dis-ances. The greater the range of the variogram model, the better
l Indicators 11 (2011) 1046–1056
the continuity of the observed variable, indicating that tree speciesrichness is more continuous (Table 2, for Co-kriging). Fig. 2a showsthe spatial distribution of tree species richness using COK.
3.4. Kriging with external drift
The parameters of the fitted spherical models for the modeledvariograms are given in Table 2, both for simple and multiple KED.The structural variance ranged from 36.7% to 56.3%, for multipleKED and simple KED respectively, suggesting that both models havea relatively large proportion of unexplained variability of the resid-ual values and that this may vary over short distances. The rangeof influence showed values between 862 m and 1054 m, so thatone could reasonably expect residual values in places separatedby such distances to be correlated for simple and multiple KED,respectively.
The model parameters of the variograms from residuals werethen used with PC1 raw data as external drift and the appropriatetrend function to predict tree species richness. Fig. 2b shows themap of simple KED prediction for the study area. Similarly, vari-ograms of residuals together with the same variables for externaldrift were used to predict species richness through multiple KED.The distribution of tree species richness using the multiple KED isshown in Fig. 2c.
3.5. Simple regression kriging
Using simple linear regression to predict tree species richnessfrom PC1, a gradient of different spectral values of remotely senseddata, produced an R2 of 0.39. This means that PC1 contributed with39% of the total variability of tree species richness (Table 5, sim-ple regression). An experimental variogram was computed on theresiduals from this regression model (Table 2, For RK using simpleregression). A spherical model fit well the experimental semi-variogram. The structural variance was 69.8% and corresponds tothe variance due to the spatial dependence explained by the model.This result suggests that the models have a relatively large propor-tion of unexplained variability of the residual values and that thismay vary over small distances. The range of influence showed val-ues of 905 m, which indicates that one would reasonably expect aspatial dependence of the residual values between places separatedby such distance. The regression model parameters were used topredict tree species richness, and the regression kriged estimateswere obtained by adding residual kriging estimates to these val-ues. The distribution of tree species richness using the simple RK isshown in Fig. 2d.
3.6. Multiple regression kriging
Model selection for the linear regressions between the explana-tory variables and tree species richness (Table 3) showed that twomodels were good predictors of tree species richness (Table 4).The best model for predicting the dependent variable includedTM3, TM4, TM5, NDVI and texture measures (relative richness anddominance) as explanatory variables (Akaike weight = 0.78). How-ever, there was an additional model different from the best onethat had Akaike weights > 0.1, thus model averaging was applied tocreate a composite model for tree species richness (Table 5, multi-ple regression). Model-averaging indicated statistically significantrelationships between tree species richness and TM3, TM4, NDVI,relative richness and dominance values (Table 5, multiple regres-sion). The dependent variable was negatively related to TM3, NDVI
and relative richness values, and positively related to the TM4 bandand dominance values.
The spatial variation depicted by the semi-variogram modelrevealed a spatial structure in the residuals of tree species rich-
For KED using simple regressionSpecies richness (residual variogram) Spherical 29.81 68.18 862.6 56.3
For KED using multiple regressionSpecies richness (residual variogram) Spherical 26.11 41.25 1054.6 36.7
For RK using simple regressionSpecies richness (regression residuals) Spherical 30.30 100.28 905.0 69.8
For RK using multiple regressionSpecies richness (regression residuals) Exponential 20.31 62.11 431.0 54.8
Fig. 2. Maps of tree species richness created by co-kriging with PC1 (a), kriging with PC1 as external drift (b), kriging with TM3, TM4, TM5, NDVI, REL-RICH and DOM asexternal drift (c), regression kriging with PC1 (d) and regression kriging with TM3, TM4, TM5, NDVI, REL-RICH and DOM (e).
IC: Akaike Information Criterion; �i: delta AIC values; wi : Akaike weights.a Only models with wi > 0.1 and with �i < 4 are shown (see Section 2).
ess (Table 2, for RK using multiple regression). An exponentialodel fit well the experimental semi-variogram, and explained the
patial autocorrelation present in the model. The structural vari-nce was 54.8%. This result suggests that the model has substantialnexplained variability of the residual values and that this mayary over small distances. The range of influence showed a valuef 431 m. This indicates that one would reasonably expect a spa-ial dependence of the residual values between places separated byuch distance. Fig. 2e shows the distribution of tree richness usinghe multiple RK.
.7. Performance of geostatistical procedures
The performance of three different approaches that combineemotely sensed data and geostatistical models for predicting treepecies richness was compared in terms of the accuracy of esti-ates. Cross-validation predictions show that the three approaches
mproved significantly the accuracy of estimates compared to ordi-
ary kriging (OK), as shown by the increased R2 and decreased rootean square errors (RMSE; Fig. 3). KED and RK models providedbetter performance as well as a visually better resolution of tree
pecies richness maps (Fig. 2) compared to COK. The R2 value in the
able 5arameters of regression models used to estimate local means of tree species richness fro
regression between observed and predicted tree species richnessusing COK (0.41) was higher than that obtained using OK (0.35).COK also decreased in RMSE (8.09) compared to OK (8.52). In turn,the R2 values increased from 0.52 to 0.62 in KED, and from 0.53 to0.66 in RK, respectively, for simple and multiple regression models.Similarly, RMSE decreased from 7.25 to 6.41 in KED, and from 7.23to 6.13 in RK, respectively, for simple and multiple regression mod-els. These results reveal the superiority of multivariate predictionmethods over those methods that used none or only one ancillaryvariable.
Predictions that considered spatial autocorrelation in the modelimproved R2 compared to the analyses that exclusively usedregression. The R2 values increased from 0.39 (Table 5) in simpleregression to 0.52 and 0.53, respectively, for simple KED and simpleRK. Similarly, the R2 values increased from 0.60 in multiple regres-sion to 0.62 and 0.66, respectively, for multiple regression KED andmultiple RK (Fig. 3). Again, multiple RK resulted in the smallestreduction in RMSE and the R2 of the model was consistently higher
than those of the other models.
The comparison of patterns of tree species richness amongvegetation classes produced by different geostatistical procedureswith those observed in the field showed important differences.
m Landsat spectral properties in regression kriging.
Fig. 3. Results of cross validation analyses used to compare the performan
or example, significant differences in tree species richness amonghe six vegetation classes were found for values estimated by COKF[5,140] = 23.88, P < 0.001), simple and multiple KED (F[5,140] = 40.42,< 0.001, F[5,140] = 50.61, P < 0.001 respectively), as well as simple
nd multiple RK (F[5,140] = 38.59, P < 0.001, F[5,140] = 53.98, P < 0.001,espectively). LSD tests showed that the six vegetation cover typesiffered significantly in estimated tree species richness (P < 0.001,able 6), with few exceptions (Table 6). Overall patterns of species
ifferent geostatistical methods for mapping tropical tree species richness.
richness among vegetation classes derived by multiple KED andmultiple RK were the most similar to field observations (Table 6).
4. Discussion
A strong limitation faced by conservation biologists and man-agers of natural resources is the lack of continuous informationconcerning species distribution patterns (Conroy and Noon, 1996;
Table 6Species richness values, both observed and estimated through three hybrid geostatistical models. Figures are means ± 1 S.E. Different letters indicate significant differenceswithin a column according to LSD test (P < 0.05).
Forest Stage 4 27.55 ± 0.63 a 23.29 ± 0.89 a 30.11 ± 0.62 a 23.48 ± 0.49 a 25.27 ± 0.78 a 25.24 ± 0.78 aForest Stage 3 24.64 ± 0.72 b 20.47 ± 1.15 b 26.50 ± 0.81 b 20.30 ± 0.71 b 19.98 ± 0.84 b 20.17 ± 0.85 bForest Stage 2 21.20 ± 0.83 c 17.77 ± 1.43 b 24.15 ± 1.23 c 18.51 ± 1.06 c 18.68 ± 1.30 b 19.04 ± 1.03 bForest Stage 1 10.35 ± 0.73 d 17.56 ± 1.72 b 19.79 ± 1.86 d 14.10 ± 1.53 d 15.64 ± 1.41 c 15.05 ± 1.44 c
8.34 ±8.11 ±
AiphesiIfbiip
giuhsudbtpmicsaeCalR
Rfstctikfi2CamiFtvf
Low flooded forest 4.35 ± 0.32 e 8.07 ± 1.53 c 1Savanna 0.65 ± 0.15 e 7.72 ± 1.49 c
raújo and Guisán, 2006). In addition to providing guidance regard-ng the selection and effectiveness of nature protection areas,recise biodiversity maps produced by accurate modeling can alsoelp to assess species responses to global climate change, and tovaluate the influence of exotic species (Rodríguez et al., 2007). Thisituation has led to multiple efforts to develop different estimat-ng procedures aimed at filling the existing gaps in our knowledge.nstead of being a real solution, the existence of a plethora of dif-erent ways to estimate unavailable information is increasinglyecoming a problem in itself, as it is not obvious which one to use
f one is to make the best possible estimation. In this context, its of utmost interest for the research community to evaluate theerformance of the different procedures that have been proposed.
The results obtained from the comparisons of three hybrideostatistical procedures for estimating plant species richness,ncorporating remotely sensed data as ancillary information, allows to draw several conclusions. First, we must emphasize that theybrid methods performed better than ordinary kriging, demon-trating an improvement in the accuracy of estimations whensing remotely sensed data indicators for species richness pre-iction based on sparse field data. Also relevant is the fact thatoth the simple and multiple regression models performed betterhan ordinary kriging, which implies strong correlations betweenlant species richness and remotely sensed data. However, itust also be stressed that the potential of hybrid methods to
mprove estimations over ordinary kriging is high only if the asso-iations between primary and ancillary variables are robust andignificant (Simbahan et al., 2006). Moreover, hybrid models thatccounted for spatial autocorrelation improved the accuracy ofstimates compared to simple and multiple regression models.ross-validation predictions showed higher R2 for COK, simple KEDnd simple RK than those based on simple regression models. Simi-arly, higher R2 values were observed for multiple KED and multipleK than for the multiple regression model.
Among the compared hybrid estimating procedures, multipleK was the most suitable method according to its superior per-
ormance for accurately mapping the spatial distribution of treepecies richness in the studied landscape. A likely explanation forhis result is that the independent variables (remotely sensed indi-ators) already accounted for 60% of species richness variability inhe multiple regression model. The performance of hybrid estimat-ng procedures has been compared for soil properties, but, to ournowledge, not for species richness. Our results are in line withndings of some studies mapping soil properties (Simbahan et al.,006), while other studies have found a better performances ofOK over RK for estimating soil properties (Rivero et al., 2007),nd a better performance of KED over COK and RK for the esti-ation of soil nitrogen availability (Baxter and Oliver, 2005). The
nconsistencies of these results have two important implications.
irst, they call attention to the need of performing specific evalua-ions of different procedures used to estimate different ecologicalariables. Second, it appears that the observed discrepancies deriverom the strength of the spatial covariance between primary and
1.51 d 13.87 ± 1.17 d 10.37 ± 0.95 d 9.64 ± 0.98 d2.24 e 5.12 ± 1.81 e 3.61 ± 0.75 e 3.55 ± 1.610 e
secondary variables, and from the ability of a regression model toexplain the variability of a target variable (Simbahan et al., 2006).An additional advantage of RK over KED is the flexibility of model-ing multivariate relationships between tree richness and remotelysensed data, since RK combines a multiple regression model withOK of regression residuals.
Tree species richness was shown to be strongly related toremotely sensed data in our landscape. Specifically, regression anal-ysis suggest that tree species richness is strongly related to thereflectance of the red and near infrared bands (TM3 and TM4),the NDVI, and to the two texture measures, namely relative rich-ness and dominance. The negative coefficient of the TM3 bandindicates a low response in the visible red wavelengths due tochlorophyll absorption, whereas the positive coefficient obtainedfor TM4 confirmed its response to green biomass and its consequentpotential for discriminating plant species based on differences intheir reflectance properties (Nagendra, 2001; Vieira et al., 2003;Thenkabail et al., 2004). It is likely that these responses are valid forvegetation during active greening only, however, and thus may notbe universally applicable to vegetation mapping (Muldavin et al.,2001). Consequently, one should exert caution in applying thesemodels for predicting species richness in other (especially drier)tropical landscapes.
As was done in this study, previous reports have highlightedthe significant association between NVDI and plant species rich-ness (Waide et al., 1999; Oindo and Skidmore, 2002; Fairbanks andMcGwire, 2004; Levin et al., 2007). NDVI, a variable commonly asso-ciated to NPP, has been hypothesized to be positively related tospecies richness, in line with species-energy theory (Evans et al.,2005; Gaston, 2000). This theory, however, was formulated for, andgenerally applies to, broad spatial scales (thousands of km2), butnot necessarily to finer scales (such as the one used in this study)where the type and form of the relationship can vary considerably(Waide et al., 1999; Grace et al., 2001; Mittelbach et al., 2001). Inthis study, we found a negative sign of NDVI in the model aver-aged parameters for species richness. Such negative relationshipcould emerge if areas of high foliage density (and thus high NDVI)are dominated by few species with deep crowns capable of limit-ing the establishment of other species (Huston, 1994; Waring et al.,2002).
The two textural variables examined by us behaved contraryto the initial expectations. In theory, relative richness, commonlyinterpreted as an indicator of habitat heterogeneity, should be pos-itively related to species richness. However, this was not the casein our landscape. Equally unexpected was the positive relation ofdominance with plant species richness. These results strongly sug-gest that habitat heterogeneity is a scale dependent phenomenon(Nagendra, 2001). The inclusion of different habitat types at land-scape scales reveals a relationship between the variation in land
cover types and habitat heterogeneity (Hernández-Stefanoni andDupuy, 2008). Similarly, when performing a local-scale analysis,factors unrelated to diversity of land cover classes are correlatedwith species richness. Since the oldest stages of forest succession
Stages 3 and 4) occupy the largest proportion of the landscape,ree species richness may be enhanced by more successful estab-ishment of the highly diverse group of shade tolerant speciesLaurance, 1991; Holl, 1999; Lebrija-Trejos et al., 2010). Of course,n a broad conservation context, the value of the herbaceous flora,
hose richness is not necessarily concentrated in mature or nearature forests should not be disregarded. Future research should
im at mapping this biodiversity component accurately.The lower prediction accuracy observed for species richness in
avannas and seasonally flooded low forest compared to semide-iduous forest may be partly attributed to scale-dependency ofiodiversity drivers (Field et al., 2009). Our finding that envi-onmental heterogeneity as gauged by remotely sensed textureeasures and productivity estimated through NDVI increases the
ccuracy of tree species richness prediction for forest vegetationoncurs with results of previous studies (Gillespie, 2006; Rocchini,007), and suggests that these factors act at relatively broad spa-ial scale allowing the use of medium resolution satellite imageryuch as Landsat TM for biodiversity estimation. Conversely, inavanna and seasonally flooded low forest, both being associatedo relatively homogeneous environment, species richness may benfluenced by other ecological factors, such as soil and micro-limatic conditions, acting at finer spatial scales. Consequently,iodiversity estimations for those areas may require higher res-lution satellite imagery (Nagendra and Rocchini, 2008) capable ofapturing variation at those scales.
Finally, when plant species richness in different habitats (veg-tation classes) exhibits very low within-class variance comparedo total variance across the landscape, as in this case (see Table 6),he final estimate will be improved if the interpolation proceduresre applied separately for each vegetation class (Burrough andcDonnell, 1998; Webster and Oliver, 2001). Using models that
pply hybrid geostatistical procedures, as those used in this study,o each stratum or land cover type may therefore improve the pre-iction accuracy, but would require a considerably large samplingffort within each stratum.
. Conclusion
Multiple RK clearly had a superior performance, not onlyecause of its greater potential to increase map accuracy, but alsoecause of its flexibility in modeling multivariate relationshipsetween tree richness and remotely sensed data. This procedureffers great promise for the production of accurate maps of speciesichness from widely accessible, low-cost, remotely sensed indica-ors. Undoubtedly, these models can be valuable tools for guidinguture efforts aimed at understanding, conserving and managingighly diverse tropical forests.
eferences
nderson, D.R., Burnham, K.P., Thompson, W.L., 2000. Null hypothesis testing:problems, prevalence, and an alternative. Journal of Wildlife Management 64,912–923.
nderson, D.R, Burnham, K.P., 2002. Avoiding pitfalls when using information-theoretic methods. Journal of Wildlife Management 66, 912–918.
raújo, M.B., Guisán, A., 2006. Five (or so) challenges for species distribution mod-elling. Journal of Biogeography 33, 1677–1688.
axter, S.J., Oliver, M.A., 2005. The spatial prediction of soil mineral N and potentiallyavailable N using Elevation. Geoderma 128, 325–339.
raakhekke, W.G., Hooftman, D.A.P., 1999. The source balance hypothesis of plantspecies diversity in grassland. Journal of Vegetation Science 10, 187–200.
urnham, K.P., Anderson, D.R., 1998. Model Selection and Inference: A PracticalInformation-theoretic Approach. Springer-Verlag, New York.
urrough, P.A., McDonnell, R.A., 1998. Principles of Geographical Information Sys-
tems. Spatial Information Systems and Geostatistics. Oxford University Press,Oxford.
abrera, C.E., Souza, M., Téllez, V.O., 1982. Imágenes de la flora Quintanaroense.Centro de Investigaciones de Quintana Roo, Mexico City.
ampbell, J.B., 1987. Introduction to Remote Sensing. The Guilford Press, New York.
l Indicators 11 (2011) 1046–1056 1055
Chavez, P.S., 1996. Image-based atmospheric corrections revisited and improved.Photogrammetric Engineering and Remote Sensing 62, 1025–1036.
Clark, D.B., Palmer, M.W., Clark, D.A., 1999. Edaphic factors and the landscape-scaledistributions of tropical rain forest trees. Ecology 80, 2662–2675.
Conroy, M.J., Noon, B.R., 1996. Mapping of species richness for conservation of bio-logical diversity: conceptual and methodological issues. Ecological Applications6, 763–773.
Dormann, C.F., McPherson, J.M., Araújo, M.B., Bivand, R., Bolliger, J., Carl, G., Davies,R.G., Hirzel, A., Jetz, W., Kissling, W.D., Kühn, I., Ohlemüller, R., Peres-Neto, P.R.,Reineking, B., Schröder, B., Schurr, F.M., Wilson, R., 2007. Methods to account forspatial autocorrelation in the analysis of species distributional data: a review.Ecography 30, 609–628.
Earth Resource Mapping Ltd., 1998. ER Mapper 6.1. User Guide, San Diego.Evans, K.L., Greenwood, J.J.D., Gaston, K.J., 2005. Dissecting the species–energy rela-
tionship. Proceedings of the Royal Society Series B 272, 2155–2163.Field, R., Hawkins, B.A., Cornell, H.V., Currie, D.J., Diniz-Filho, A.F., Gueran, J.F., Kauf-
Fairbanks, H.K., McGwire, K.C., 2004. Patterns of floristic richness in vegetation com-munities of California: regional scale analysis with multi-temporal NDVI. GlobalEcology and Biogeography 13, 221–235.
Gaston, K.J., 2000. Global patterns in biodiversity. Nature 405, 220–227.Georgakarakos, S., Kitsiou, D., 2008. Mapping abundance distribution of small pelagic
species applying hydroacoustics and Co-kriging techniques. Hydrobiologia 612,155–169.
Gould, W., 2000. Remote sensing of vegetation, plant species richness, and regionalbiodiversity hotspots. Ecological Applications 10, 1861–1870.
Gillespie, T.W., 2006. Predicting woody-plant species richness in tropical dryforests: a case study from South Florida, USA. Ecological Applications 15, 27–37.
Gillespie, T.W., Foody, G.M., Rocchini, D., Giorgi, A.P., Saatchi, S., 2008. Measuringand modelling biodiversity from space. Progress in Physical Geography 32, 203–221.
Gotelli, N.J., Colwell, R.K., 2001. Quantifying biodiversity: procedures and pitfallsin the measurement and comparison of species richness. Ecology Letters 4,379–391.
Goovaerts, P., 1999. Geostatistics in soil science: state-of-the-art and perspectives.Geoderma 89, 1–45.
Grace, J., Malhi, Y., Higuchi, N., Meir, P., 2001. Productivity and carbon fluxes of tropi-cal rain forests. In: Mooney, H., Saugier, B. (Eds.), Terrestrial Global Productivity:Past, Present and Future. Academic Press.
Hernández-Stefanoni, J.L., Ponce-Hernández, R., 2004. Mapping the spatial dis-tribution of plant diversity indices using multi-spectral satellite imageclassification and field measurements. Biodiversity and Conservation 13, 2599–2621.
Hernández-Stefanoni, J.L., Bello-Pineda, J., Valdés-Valadez, G., 2006. Comparing theuse of indigenous knowledge with classification and ordination techniques forassessing the species composition and structure of vegetation in a tropical forest.Environmental Management 37, 686–702.
Hernández-Stefanoni, J.L., Dupuy, J.M., 2007. Mapping species density of trees,shrubs and vines in a tropical forest, using field measurements, satellite mul-tispectral imagery and spatial interpolation. Biodiversity and Conservation 16,3817–3833.
Hernández-Stefanoni, J.L., Dupuy, J.M., 2008. Effects of landscape patterns on speciesdensity and abundance of trees in a tropical subdeciduos forest of the YucatanPeninsula. Forest Ecology and Management 255, 3797–3805.
Hernández-Stefanoni, J.L., Dupuy, J.M., Castillo-Santiago, M.A., 2009. Assessingspecies density and abundance of tropical trees from remotely sensed data andgeostatistics. Applied Vegetation Science 12, 398–414.
Holl, K., 1999. Factors limiting tropical rain forest regeneration in abandonedpastures: seed rain, seed germination, microclimate, and soil. Biotropica 31,229–242.
Huggett, R.J., 1995. Geoecology: An Evolutionary Approach. Routledge, New York.Huston, M.A., 1994. Biological Diversity: The Coexistence of Species. Cambridge
University Press, Cambridge.Isaaks, E.H., Srivastava, R.M., 1989. An Introduction to Applied Geostatistics. Oxford
University Press, New York.Jensen, J.R., 2000. Remote Sensing of Environment: An Earth Resource Perspective.
Prentice Hall, New Jersey.Johnson, J.B., Omland, K.S., 2004. Model selection in ecology and evolution. Trends
in Ecology and Evolution 19, 101–108.Laurance, W.F., 1991. Edge effects in tropical forest fragments: application of a model
for the design of nature reserves. Biological Conservation 57, 205–219.Lebrija-Trejos, E., Pérez-García, E.A., Meave, J.A., Bongers, F., Poorter, L., 2010.
Functional traits and environmental filtering drive community assembly in aspecies-rich tropical landscape. Ecology 91, 386–398.
Levin, N., Shimida, A., Levanoni, O., Tamari, H., Kark, S., 2007. Predicting mountainplant richness and rarity from space using satellite-derived vegetation indices.Diversity and Distributions 13, 1–12.
Webster, R., Oliver, M.A., 2001. Geostatistics for Environmental Science. John Wiley
056 J.L. Hernández-Stefanoni et al. / Eco
utanga, O., Rugege, D., 2006. Integrating remote sensing and spatial statisticsto model biomass distribution in a tropical savanna. International Journal ofRemote Sensing 27, 3499–3514.
agendra, H., 2001. Using remote sensing to assess biodiversity. International Jour-nal of Remote Sensing 22, 2377–2400.
agendra, H., Gadgil, M., 1999. Satellite imagery as a tool for monitoring speciesdiversity: an assessment. Journal of Applied Ecology 36, 388–397.
agendra, H., Rocchini, D., 2008. High resolution satellite imagery for tropical bio-diversity studies: the devil is in the detail. Biodiversity and Conservation 17,3431–3442.
ational Aeronautics and Space Administration, 1998. Landsat 7 Science DataUsers Handbook Greenbelt, Maryland. Goddard Space Flight Center, electronicversion available at: http://landsathandbook.gsfc.nasa.gov/handbook.html, lastaccessed 02.07.10.
deh, I.O.A., McBratney, A.B., Chittleborough, D.J., 1994. Spatial prediction of soilproperties from landform attributes derived from a digital elevation model.Geoderma 63, 197–214.
indo, B.O., Skidmore, A.K., 2002. Interannual variability of NDVI and species rich-ness in Kenya. International Journal of Remote Sensing 23, 1195–1198.
ldeland, J., Wesuls, D., Rocchini, D., Schmidt, M., Jürgens, N., 2010. Does usingspecies abundance data improve estimates of species diversity from remotelysensed spectral heterogeneity? Ecological Indicators 10, 390–396.
rellana, R., Espadas, C., Conde, C., Gay, C., 2009. Atlas. Escenarios de CambioClimático en la Peninsula de Yucatán. Centro de Investigación Científica deYucatán, Mérida.
almer, M.W., Clark, D.B., Clark, D.A., 2000. Is the number of tree species in smalltropical forest plots nonrandom? Community Ecology 1, 95–101.
ebesma, E.J., 2004. Multivariable geostatistics in S: the gstat package. Computers &Geosciences 30, 683–691.
lotkin, J.B., Potts, M.D., Yu, D.W., Bunyavejchewin, S., Condit, R., Foster, R., Hubbell,S., LaFrankie, J., Manokaran, N., Lee, H.S., Sukumar, R., Nowak, M.A., Ashton, P.S.,2000. Predicting species diversity in tropical forests. Proceedings of the NationalAcademy of Science USA 97, 10850–10854.
ivero, R.G., Grunwald, S., Bruland, G.L., 2007. Incorporation of spectral data intomultivariate geostatistical models to map soil phosphorus variability in a Floridawetland. Geoderma 140, 428–443.
occhini, D., 2007. Effects of spatial and spectral resolution in estimating ecosystemá-diversity by satellite imagery. Remote Sensing of Environment 111, 423–434.
l Indicators 11 (2011) 1046–1056
Rocchini, D., Balkenhol, N., Carter, G.A., Foody, G.M., Gillespie, T.W., He, K.S., Kark,S., Levin, N., Lucas, K., Luoto, M., Nagendra, H., Oldeland, J., Ricotta, C., South-worth, J., Neteler, M., 2010. Remotely sensed spectral heterogeneity as a proxyof species diversity: recent advances and open challenges. Ecological Informatics5, 318–329.
Rodríguez, J.P., Brotons, L., Bustamante, J., Seoane, J., 2007. The application of pre-dictive modelling of species distribution to biodiversity conservation. Diversityand Distributions 13, 243–251.
Simbahan, G.C., Dobermann, A., Goovearts, P., Ping, J., Haddix, M., 2006. Fine-resolution mapping of soil organic carbon based on multivariate secondary data.Geoderma 132, 471–489.
Sales, M.H., Souza, C.M., Kyriakidis, P.C., Roberts, D.A., Vidal, E., 2007. Improvingspatial distribution estimation of forest biomass with geostatistics: a case studyfor Rondonia, Brazil. Ecological Modeling 205, 221–230.
Thenkabail, P.S., Enclona, E.A., Ashton, M.S., Legg, C., De Dieu, M.J., 2004. Hyperion,IKONOS, ALI and ETM+ sensors in the study of African rainforest. Remote Sensingof Environment 90, 23–43.
Turner, M.G., 1989. Landscape ecology: the effect of pattern on process. AnnualReview of Ecology and Systematics 20, 171–197.
Vieira, I.M.G., De Almeida, A.S., Davidson, E.A., Stone, T.A., Carvalh, C.J.R., Guerrero,J.B., 2003. Classifying successional forest using Landsat spectral properties andecological characteristics in eastern Amazonia. Remote Sensing of Environment87, 470–481.
Waide, R.B., Willig, M.R., Steiner, C.F., Mittelbach, G., Gough, L., Dodson, S.I., Juday,G.P., Parmenter, R., 1999. The relationship between productivity and speciesrichness. Annual Review of Ecology and Systematics 30, 257–300.
Waring, R., Coops, N.C., Ohmann, J.L., Sarr, D.A., 2002. Interpreting woody plantrichness from seasonal ratios of photosynthesis. Ecology 83, 2964–2970.
Wagner, H.H., Wildi, O., Ewald, K.C., 2000. Additive partitioning of plant speciesdiversity in an agricultural mosaic landscape. Landscape Ecology 15, 219–227.
and Sons, Toronto.Wilson, S.D., 2000. Heterogeneity, diversity and scale in plant communities. In:
Hutchings, M.J., John, E.A., Stewart, A.J.A. (Eds.), The Ecological Consequencesof Environmental Heterogeneity. Blackwell Science, Oxford, pp. 53–69.
