Assessment of geostatistical and interpolation methods for mapping forest dieback intensity in Zagros forests

Document Type: Research Paper

Authors

1 Sari University

2 Gorgan University

3 University of Wuerzburg

Abstract

During recent years, oak decline has been widely spread across Brant’s oak (Quercus Brantii Lindl.) stands in the Zagros Mountains, Western Iran, which caused large-area forest dieback in several sites. Mapping the intensity and spatial distribution of forest dieback is essential for developing management and control strategies. This study evaluated a range of geostatistical and interpolation methods to explore the spatial structure and provide area-based maps of the intensity of forest dieback across a representative test site - Ilam Province - that was severely affected by Oak decline. The geostatistical analysis provided in-depth measures of the spatial structure amongst the selective sampling units (120 quadratic sample plots of 1200 m2), which eventually resulted in an area-based maps of dieback intensity. The accuracy of the applied methods was assessed by mean error percentage (%ME), root mean squared error percentage (%RMSE) and coefficient of determination (R2). Results showed moderate spatial structure within the sampling units. Moreover, cokriging (associated with soil humidity and aspect as independent variables) approach resulted in the highest accuracy, followed by two other methods of kriging and Radial Basis Function. Results suggested that cokriging can accurately estimate the intensity of dieback and its spatial distribution in the study area. According to this, an average dieback intensity of 18.12 % was estimated within the study area.

Keywords


Assessment of geostatistical and interpolation methods for mapping forest dieback intensity in Zagros forests

 

Karami O.1, Fallah A.*1, Shataei S.H.2, Latifi H.3

 

1. Department of Forest Sciences, Faculty of Natural Resources, Sari University of Agriculture and Natural Resources, Sari, Iran.

2. Department of Forestry, Faculty of Forest Sciences, Gorgan University of Agriculture and Natural Resources, Gorgan, Iran.

3. Institute of Geography and Geology, University of  Wuerzburg, Oswald-Külpe-Weg 86, D-97074, Wuerzburg,Germany.

 

* Corresponding author’s E-mail: fallaha2007@yahoo.com

(Received: Oct. 07. 2017 Accepted: Feb. 14. 2018)

ABSTRACT

During recent years, oak decline has been widely spread across Brant’s oak (Quercus Brantii Lindl.) stands in the Zagros Mountains, Western Iran, which caused large-area forest dieback in several sites. Mapping the intensity and spatial distribution of forest dieback is essential for developing management and control strategies. This study evaluated a range of geostatistical and interpolation methods to explore the spatial structure and provide area-based maps of the intensity of forest dieback across a representative test site - Ilam Province - that was severely affected by Oak decline. The geostatistical analysis provided in-depth measures of the spatial structure amongst the selective sampling units (120 quadratic sample plots of 1200 m2), which eventually resulted in an area-based maps of dieback intensity. The accuracy of the applied methods was assessed by mean error percentage (%ME), root mean squared error percentage (%RMSE) and coefficient of determination (R2). Results showed moderate spatial structure within the sampling units. Moreover, cokriging (associated with soil humidity and aspect as independent variables) approach resulted in the highest accuracy, followed by two other methods of kriging and Radial Basis Function. Results suggested that cokriging can accurately estimate the intensity of dieback and its spatial distribution in the study area. According to this, an average dieback intensity of 18.12 % was estimated within the study area.

 

Key words: Oak decline, Spatial structure, Interpolation, Geostatistics, Zagros.


INTRODUCTION

Zagros forests in western Iran play an important role in balancing regional water resources, as well as in soil conservation, weather modification and socio-economic stability (Heydari 2006).

In recent years, these forests have been exposed to decline caused by natural and anthropological agents, which resulted in large diebacks in oak stands.

The oak decline was initially reported in 2008 by local forestry offices in Ilam and Lorestan provinces. Ever since numerous reports have been published on dieback of oak trees in multiple spatial scales in Western Iran, mostly indicating a rapid spread throughout the Zagros forests over the time (Mir-Abolfathi 2013).

The phenomenon occurs as a result of  a complex set of biotic and abiotic factors such as biological disturbance agents (pests, diseases), physiographic conditions, soil, physiologic tree mortality, periodic droughts, dust and socio-economic issues (Fan et al. 2006; Baguska et al. 2014). Regardless of the causal factors, the consequences of this phenomenon initially include severe structural changes which eventually degrade variety of forest ecosystem functions (Hosseini et al. 2014).

Access to up-to-date spatiotemporal information is a key factor to develop and maintain monitoring strategies for forest health. That is, one of the crucial issues in predicting the forest reaction to periodic decline is quantification of dieback intensity of trees, which can further help by facilitating the optimization of management strategies (Sánchez et al. 2002).

Therefore, knowledge on spatial distribution  of forest calamities (e.g. stand decline) as mapped on area-based level is considered as a basic prerequisite for forest health monitoring (Ristainio & Gumpertz 2000), for which a range of various methods exists (Holdenrieder et al. 2004).

Whereas common logistic and financial constrains in conventional field-based inventories make timely description of spatial variation of forest disturbances infeasible, these  can be alternatively estimated by means of interpolating previously sampled attributes on sparse locations to those of unknown attributes (Burrough, 2001).

In recent years, use of different interpolation methods from geostatistical field for mapping the spatial distribution of various environmental entities has been substantially increased. Geostatistical methods help by detecting, estimating and mapping the spatial patterns of regional response variables by focusing on modeling and interpretation of semivariograms (Lopez-Granados et al. 2002). While in these methods the location of the samples play a major role in analyzing the unknown values on known locations, the spatial location is largely neglected in classical statistics (Pourreza et al. 2012).

To our knowledge, geostatistical and interpolation methods have been only rarely applied for mapping the spatial distribution of forest decline, whereas the existing works suggest that these methods can provide good assessments of forest decline (Hlasni et al. 2009). The early works by Hohn et al. (1993) and Liebhold et al. (1993) presented the use of kriging techniques for estimating the spatial spread of Gypsy Moth (Lymantria dispar dispar L.). Koch and Smith (2008) used kriging approach for mapping the spatial pattern of dead trees and concluded kriging can accurately estimate the spatial distribution of dead trees.

In addition, geostatistical methods were evaluated for mapping forest defoliation, in which cokriging was shown to return superior results compared to kriging (Cocco et al. 2010). Foster et al. (2013) studied the state of defoliation caused by insect attack, and concluded that modeling of the defoliation along with other correlated factors results in more accurate estimations. Despite a number of studies on geostatistical estimation of forest and soil attributes in Iran (Habashi et al. 2007; Akhavan et al. 2010; Hosseini et al. 2012; Jafarnia & Akbarinia, 2014), there is a lack of reports on geostatistical applications for estimating forest health condition, in particular across semi-Mediterranean Zagros forests in Western Iran.

This study aims to investigate the relationship between oak decline and a range of potentially-influential physiographic and soil factors, followed by evaluating geostatistical and interpolation methods for estimating and spatially mapping oak decline in a representative portion of forests affected by oak decline.

 

MATERIALS AND METHODS

Study area

The study site encompasses an area of ca. 1339 ha in northeast of Ilam City in Ilam Province (longitude: 46˚ 25' 39" - 46˚ 29' 45" E and latitude: 33˚ 40' 4" - 33˚ 37'18" N) (Fig. 1). The slope average is 37%, with elevation ranging between 1440 and 2324 m a.s.l. The average annual rainfall is 562.7 mm, and average annual temperature is 16.8° C.

The study area is a part of the Zagros mountain folds that were formed in the late Triassic period (Mirzaei 2005).

 

 

Fig. 1. The location of the study site in Ilam Province (bottom right), country (top right), overlaid by the sample plots.

 

METHODS

Field sampling

A selective sampling (McCoy 2005) was used as a result of multiple site visits and according to the spatial distribution of dead trees in the study area, in which 120 quadratic sample plots of 1200 m2 each  in  a forest area of ca. 1339 ha in northeast of  Ilam in Ilam Province were inventoried. Within each sample plot dieback intensity, tree crown diameters as well as number of trees were measured. The plot center coordinates were recorded by differential GPS surveys. The sampled data were partitioned into training and test datasets by assigning 70% of samples for training and keeping 30% for testing.  The recorded sample plots were exported into a vectorized shapefile. In addition, a digital elevation model (DEM) was derived from the local topographic map and was used to extract plot-based average values of topographic attributes including slope, elevation and aspect (Cocco et al. 2010).

In each plot, a mixed soil sample from plot center and four corners was taken from 0-20 cm soil depth (Robertson et al. 1999), which was further analyzed in the soil lab (Maranon et al. 1999). The soil samples were air-dried and sieved through a 2-mm mesh sieve. Soil pH was measured with a 1/2.5 ratio using pH meter. Electric conductivity (EC) was measured with a

 

 

1/5 ratio, total nitrogen was measured by Kjeldahl method (Westerman et al. 1990), while phosphorus was measured following Olsen and Sommers (1990). Available potassium was extracted using ammonium acetate and calcium carbonate was measured by neutralizing with hydrochloric acid. In addition, the amount of organic matter was measured by its oxidation (Walkley & Black 1934), bulk density was measured on core soil samples by drying in a low-temperature oven (105 °C) for at least 24 h, the soil moisture contents was measured by oven-drying method and soil texture was determined by the percentage of clay, sand and silt using hydrometric method (Duffera et al. 2007).

 

Statistical preprocessing

Normal distribution of the underlying data is a commonly considered as an assumption for using geostatistical variograms. Therefore, we initially conducted a normality by checking Kolmogorov-Smirnov test (Vannini et al. 2010). We additionally correlated between the dieback intensity and the physiography -ic/edaphic parameters to provide initial impressions prior to the use of cokriging approach (Erre et al. 2009: Cocco et al. 2010).

Analysis of spatial structure

Semivariogram (also referred to as variorum) is a basic tool in geostatistics which is applied to reveal the spatial correlation existing within a given sample set. It is a function which describes the degree of spatial dependence within the data, and is defined as the expected squared increment of the recorded values between two locations (Wackernagel 2003). The spatial correlation amongst the measurement points can be quantified by the following semi-variance function (equation 1) and co-semivariogram (equation 2) (Vannini et al. 2010):

 

 

                                                (1)

                (2)


 

Where γ(h) is the semivariance for the interval distance class h, N(h) is the number of pairs of the lag interval, Z(xi)  is the measured sample value at point i, Z(xi+h) is the measured sample value at position i+h and Z2(xi) and Z2(xi+h) are related to second variable.

If variogram depends on a direction, it is known as anisotropic and otherwise isotropic. A variogram consists of three major characteristics including sill (C + C0), range (R) and nugget semivariance (C0). Nugget is the variance at zero distance, sill which defines the asymptotic value of semi-variance with respect to the lag distance, and range is the distance at which values of one variable reaches spatial independence. To define different classes of spatial dependence for soil variables, the ratio between the nugget semivariance and the total semivariance or sill can be used (Cambardella & Karlen 1994), in which the variable is subject to strong spatial dependence if the ratio is ≤ 25%. The variable is considered to be moderately distributed in patches if the ratio is between 25 and 75%, and it is weakly spatially dependent if the ratio is > 75%.

The variable have not spatially structure if the ratio = 100%. Finally the variable is considered to be spatially uncorrelated (pure nugget), if the semivariogram slope is close to zero (Lopez-Granados et al. 2002).

 

Mapping dieback intensity

We evaluated two geostatistical (kriging and cokriging) and a range of interpolation

 

approaches (Inverse Distance Weighted, Radial Basis Functions, General Polynomial and Local Polynomial) to map dieback intensity in the study area. A summary of the applied methods is as follows:

Kriging: kriging is a technique for unbiased estimations of variables at unsampled spatial locations. It is known as an exact interpolator which is applicable when no assumption is made on the trend in the data. Moreover, it is applied when the variogram function depends on the separation vector and not on location (Nalder & Wein 1998).

It is a method of local estimation in which each estimate is a weighted average of the observed values.  The spatial prediction of the values of a variable Z at an unsampled point x0 is given by the following formula (Mardiks et al. 2005):

 

                               (3)

Where xdenotes a set of spatial coordinates (X, Y), λi is the weights associated with the sampling point xi and the ith observation point.

Cokriging: Collocated kriging is defined as a multivariate version of kriging which uses additional covariates that are ideally sampled at the same location as the estimated variable. When spatial correlation between a covariate and the variable of interest is high, cokriging commonly returns better estimation results. To apply cokriging (equation 4) for modelling the relationship between the response and predictors variables, cross semivariogram is used (Leu et al. 2008).

              (4)   

Where is the estimated value for xi, λi is the weights of the variable x, λj defines the weights of the variable y, is the observed value of the response variable and  is the observed value of the predictors.

Inverse Distance weighted (IDW):IDW combines the idea of proximity espoused by Thiessen polygons (Thiessen 1911) with the gradual change of a trend surface. Therefore, the measured values which are spatially closest to the prediction location, will have higher influence on the predicted value.

IDW assumes that each measured point has a local influence that diminishes with distance (Leu et al. 2008).

Radial Basis Function (RBF): RBF embraces a diverse range of interpolation methods. In general, RBF is a real-valued function in which the value depends only on the distance from the origin. Sums of RBFs are typically used to approximate given functions.

This approximation process can also be interpreted as a simple neural network. In the RBF, values are estimates over the maximum of observed values and are in some cases less than the minimum of observed (Mir-Mousavi et al. 2010).

General Polynomial Interpolation (GPI): Under GPI mapped data are approximated by a polynomial expansion of the geographic coordinates of a set of recorded control points. The coefficients of the polynomial function are found by the least squares method, insuring that the sum of the squared deviations is minimized from the trend surface.

Each original observation is considered to be the sum of a deterministic polynomial function of the geographic coordinates plus a random error term (Leu et al. 2008).

Local Polynomial Interpolation (LPI): LPI is mostly similar to GPI. It fits multiple polynomials using subsets of the measured points, whereas GPI fits a polynomial model to the entire surface based on all measured points. LPI is a moderately quick deterministic interpolator that is smooth. It is more flexible than GPI, with more parameter decisions. When the input dataset exhibits short-range variation, LPI can be an appropriate method to capture finer detail (Akima 1970).

 

Accuracy assessment

To assess the accuracy of estimations across the applied method, 30% of the data was used as independent test dataset. The observed and estimated data were compared by deriving a set of accuracy measures including correlation coefficient (R), mean error percentage (%ME) (equation 5) and the root mean square error percentage (%RMSE) (equation 6) (Machiwal & Jha 2015).

 

                           

                                                                             (5)

 

(6)

Where  is the observed value at the point i and  is the estimated value at the point i.

 

RESULTS

Table 1 shows the results of the statistical analysis across the measured field samples. Average dieback intensity in the measured samples was 23.9%, the skewness was < 1 and Kolmogorov-Smirnov test revealed a normal distribution of the data. This was also supported by plotting the histogram of data (Fig. 2).

 

 

 

 

Table 1. Statistical characteristics of dieback intensity samples.

variable

min

mean

max

Standard deviation

skewness

kurtosis

Kolmogorov-Smirnov test

P-value

Dieback intensity

0.9

23.9

55.46

13.33

0.34

2.13

0.33

 

 

Fig. 2. Histogram of dieback intensity samples.

 

 

The correlation analysis between the dieback intensity, physiographic and soil parameters showed the correlation of dieback intensity with number of trees in plot, slope, aspect, soil bulk density and soil moisture (Table 2).

The analysis of spatial structure of the samples by means of the variogram and covariogram of dieback intensity with other correlated parameters using models such as spherical, circular, Gaussian revealed that Gaussian model fit on the spatial structure analysis of dieback intensity, which returned the lowest rates of %ME and the %RMSE (Fig. 3). “Geostatistical Analyst” module was used and different models were evaluated on ArcGIS 10.3 software. Table 3 lists variogram and covariogram parameters along with the fitted models, which indicates C0/(C+C0) of dieback intensity = 0.48. This suggests a moderate spatial structure of samples (Table 3). Furthermore, the variogram of measured samples of dieback intensity was shown to be isotropic (Fig. 3).

 

 

Table 2. The Pearson correlation between dieback intensity and other parameters.

Variable

Trees number

slope

aspect

Soil bulk density

Soil humidity

Dieback intensity

Pearson coefficient

0.197*

 

0.202*

0.344**

0.278**

0.257**

P-value

0.034

0.03

0.001

0.002

0.005

 

Table 3. variogram and covariogram parameters and fit models.

C0/(C+C0)

Nugget (C0)

Sill (C+C0)

Range (A)

model

variable

0.48

90.9

186.3

611.8

Gaussian

Dieback intensity

0.48

90.9

186.3

611.8

Gaussian

Dieback intensity with Aspect

0.51

96.3

186.3

624.9

Gaussian

Dieback intensity with bulk density

0.51

94.1

184.5

582.2

Gaussian

Dieback intensity with humidity

0.52

91.9

174.4

642.8

J-Bessel

Dieback intensity with trees number

 

 

The model fits for the applied geostatistical and interpolation approaches were compared with results indicating that cokriging, kriging and RBF provide more accurate results compared to other tested methods. cokriging (with soil humidity and aspect) using Gaussian model turned out to be the superior estimator, with average ME percentage = 0.41, average RMSE percentage = 33.30, R2 = 0.52. In contrast, the GPI (average ME percentage = 0.22, average RMSE percentage = 43.41, R2 = 0.17) was associated with the lowest accuracy. The results showed that cokriging associated with two independent variables provide more accurate results compared to cokriging associated with one independent variable. So that the results of cokriging (associated with soil humidity and aspect) had more accuracy than cokriging (associated with aspect), cokriging (associated with soil bulk density), cokriging (associated with soil humidity) and cokriging (associated with trees number) (Table 4).  Fig. 4 shows the difference between measured and predicted dieback intensity in cokriging method (associated with soil humidity and aspect as independent variables). This Fig. indicated that the predicted dieback intensity was biased compared to the measured data due to the mismatch of its line graph with measured values of dieback intensity.

 

 

Fig. 3. a: variogram of dieback intensity, b: covariogram between dieback intensity and aspect, c: covariogram between dieback intensity and soil bulk density, d: covariogram between dieback intensity and soil humidity, e: covariogram between dieback intensity and number of trees, f: variogram surface that was shown to be isotropic.

 

Table 3. variogram and covariograms parameters and fit models.

C0/(C+C0)

Nugget (C0)

Sill (C+C0)

Range (A)

model

variable

0.48

90.9

186.3

611.8

Gaussian

Dieback intensity

0.48

90.9

186.3

611.8

Gaussian

Dieback intensity with Aspect

0.51

96.3

186.3

624.9

Gaussian

Dieback intensity with bulk density

0.51

94.1

184.5

582.2

Gaussian

Dieback intensity with humidity

0.52

91.9

174.4

642.8

J-Bessel

Dieback intensity with trees number

 

Fig. 5 shows the mapped results of forest dieback intensity across the study area. A visual comparison of different mapping results obviously indicates the difference amongst the applied methods.

Using cokriging for estimation of dieback intensity associated with soil humidity and aspect as independent variables returned the highest accuracy (Fig. 5, Table 5).

 

 

Fig. 4. Correlation between measured and predicted dieback intensity in cokriging associated with soil humidity and aspect as predictors.

 

 

 

Fig. 5. The maps of dieback intensity prepared using different methods (a: GPI, b: LPI, c: IDW, d: RBF, e: kriging, f: cokriging with aspect, g: cokriging with soil bulk density, h: cokriging with soil humidity, i: cokriging with trees number, j: cokriging with humidity and aspect together).

Table 4. Assessment of the applied methods in the mapping of dieback intensity.

R2

%RMSE

%ME

Model

Method

0.17

43.41406

0.2215

Order 3

GPI

0.29

39.68547

0.479917

Order 2

LPI

0.34

39.02097

6.054341

Power 1.59

IDW

0.38

37.50738

2.916421

Regularized spline

RBF

0.40

36.80597

1.144418

Gaussian

Kriging

0.44

35.69846

1.070585

Gaussian

Cokriging (with aspect)

0.43

35.9938

0.295334

Gaussian

Cokriging (with soil bulk density)

0.47

35.88305

0.812168

Gaussian

Cokriging (with soil humidity)

0.42

36.28913

0.369167

J-Bessel

Cokriging (with trees number)

0.52

33.29888

0.406084

Gaussian

Cokriging (with soil humidity and aspect)

                                       %ME: mean error percentage

                                       %RMSE: root mean square error percentage

                                       R2: coefficient of determination

 

Table 5. The area of dieback intensity by cokriging (with soil humidity and aspect).

Classes

Dieback intensity (%)

Area (ha)

Area (%)

C1

0 -10

140.47

10.49

C2

10 - 20

781.41

58.38

C3

20 -30

309.92

23.15

C4

30 -40

103.3

7.72

C5

40 <

3.42

0.25

 

DISCUSSION

This study was designed for mapping of dieback intensity in the Zagros forest, west of Iran by means of geostatistical and interpolation methods. Statistical analysis of measured samples indicated the normal distribution of data, with a low degree of skewness. This (and similar conditions) can be attributed to intrinsic properties of the data, sampling method and number of samples (Jafarnia & Akbarinia 2014). In this study, distribution patterns and intrinsic characteristics of the oak decline supported by several site visits led to the decision using a selective sampling, which was eventually supported to be representative by the low level of skewness in the distribution of the field samples. The results did not show any dieback intensity > 42% and < 3% for the cokriging and kriging approaches, while the maximum and minimum dieback intensity of 55.46 and 0.9 % were observed in the recorded field samples, respectively. The changed maximum and minimum of data in result of geostatistic was also mentioned by previous studies by Akhavan & Kleinn (2009) and Pourreza et al. (2012). The reason can be attributed to the intrinsic smoothness of geostatistical methods (Yamamoto 2005; Pourreza et al. 2012). In contrast, the maximum estimated value of 55.74% was observed when using RBF method, which was larger than maximum value measured in the field samples. This is one of the attributes of the RBF in which the estimates can occasionally fall beyond the range of the observed values (Maroufi et al. 2009). This, in turn, caused errors which were eventually influential on the minimum and maximum margins of the area-based mapped values (Yamamoto 2010; Pourreza et al. 2012). The results of this study showed a moderate spatial structure of the forest dieback intensity in the study area. A range of different results has been reported by the available literature, which turn out to be mostly site - response variable and sampling scheme - dependent. For instance, weak spatial structure was observed amongst samples in a study carried out on quantitative characteristics across northern forests of Iran (Akhavan et al. 2010). The studies of Gunnarson et al. (1998) and Touminen et al. (2003) in Nordic Europe both reported weak spatial structure amongst the samples, whereas Kint et al. (2003) (temperate plantations), Montes et al. (2005) (natural oak stands in Spain), Akhavan & Kleinn (2009) (even-aged stands in Northern Iran) and Pourreza et al. (2012) (Zagros forests of Iran) reported moderate spatial structures. We argue that the homogeneity existing within the Zagros forests (in particular in tree species composition which is mostly dominated by Q. brantii) contribute to improving the spatial structure of several forest attributes, including those estimated here. Calculating variogram parameters resulted in the effect range of 611.8 m for forest dieback intensity. That is, the samples have spatial dependence up to 611.8 m which afterwards turns into independence. Effect range also represents the degree of homogeneity of a given variable. That is, higher values of the range lead to decrease in the number of homogeneous samples. Here, deriving the effect range was further applied to design the sampling grid, i.e.  the sampling distance was equal to the derived range (Akhavan & Kleinn 2009). In practice, however, we suggest the 2:3 of the effect range to be taken as sampling distance to account for the uncertainties raised by this (Hassani-Pak 2006; Akhavan & Kleinn 2009). In addition, according to the isotropy that was observed in the variogram, 408 m is recommended as appropriate dimensions of sampling grid for analysis of dieback intensity across this and structurally-similar study sites. The results showed that the highest accuracy in estimating the oak decline intensity was returned by cokriging approach. Cokriging uses correlation and spatial relationships amongst the response and predictor variables, and provide superior results compared to other estimated methods. This result is consistent with results from Cocco et al. (2010), Foster et al. (2012), Misaghi & Mohammadi (2005), Nourzadeh-Haddad et al. (2013), Sarmadian & Taghizadeh (2008) and Wu et al. (2006). Also the results showed that geostatistical methods were generally associated with higher accuracies than the applied interpolation methods, which supports previously published results by Mir-Mousavi et al. (2010), Kazemi-Poshtmasari et al. (2012) and Omidvar et al. (2014). The geostatistical methods did not contain systematic errors which are typically present in interpolation results, and also are estimated with minimum amount of bias and variance. Concerning this, our result is in contrast to those published by Maroufi et al. (2009) and Mohamadi et al. (2012) who reported the IDW to be of superior estimation performance. It turns out that the choice of an appropriate method for estimating a variable depends mainly on the variable type as well as on the regional influential factors. Therefore, an all-round functional estimation method cannot be prescribed due to area - specific character -istics. Also the results showed that cokriging using two predictor variables returned better performances compared to using a single predictor (Table 4).

The IDW method for estimating the spatial structure has been previously suggested by the instructions provided by Iranian Forest, Rangeland and Watershed Management Organization (FRWO) to prevent and control oak decline (Anonymous 2012). However, our results showed that kriging and cokriging methods provide more accurate results compared to IDW. Therefore we suggest to use those methods instead of IDW in the areas that are structurally similar to ours. Information on spatial extension of declined oak classes (Table 5) indicated that oak decline has been expanded and caused extensive diebacks across the study area. According to the results, the average dieback intensity in the study area was estimated to be 18.12%, which can be considered as an alarm. Together with other possible anthropological factors such as charcoal and animal husbandry, oak decline can also be considered as an emerging natural disturbance agent across Zagros forests in Iran. Thus it intensifies the call for well-founded planning to manage the Zagros forest landscapes, in particular those prone to oak decline. A future prospect of this study includes using remote sensing techniques for large - area sampling, which can help by solving problems caused by time and labor constrains and enable multi - temporal investigations. The results will be further used for estimating and mapping forest attributes by means of geostatistical methods.

 

ACKNOWLEDGEMENTS

Authors gratefully acknowledge Aliakbar Jafarzadeh, Soroush Zabiholahi, Vahid Mirzaiizade, Ahmad Amuzad, Reza Badrkhani and Reza Abedinpour for their valuable contributions in sampling and inventorying.

Akhavan, R & Kleinn, C 2009, On the potential of kriging for estimation and mapping of forest plantation stock (Case study: Beneshki plantation). Iranian Journal of Forest and Poplar Research, 17: 303-318.

Akhavan, R, Zahedi Amiri, Gh & Zobeiri, M 2010, Spatial variability of forest growing stock using geostatistics in the Caspian region of Iran. Caspian Journal of Environment Sciences, 8: 43-53.

Akima, H 1970, A new method of interpolation and smooth curve fitting based on local procedures. Journal of Association for Computing Machinery, 17: 589–602.

Anonymous, 2012, Instructions of forest sustainable management in Zagros forest ecosystems for prevention and control of oak dieback. Organization of Forests, Rangelands and Watershed Management, Tehran, Iran,  61 p.

Baguskas, SA, Peterson, SH, Bookhagen, B & Still, CJ 2014, Evaluating spatial patterns of drought - induced tree mortality in a coastal California pine forest. Forest Ecology and Management, 315: 43–53.

Borrough, PA 2001, GIS and geostatistics: Essential partners for spatial analysis, Environmental and Ecological Statistics, 8: 361 –377.

Cambardella, CA & Karlen, DK 1999, Spatial analysis of soil fertility parameters. Precision Agriculture, 1: 5–14.

Cocco, A, Cossu, AQ, Erre, P, Nieddu, G & Luciano, P 2012, Spatial analysis of gypsy moth populations in Sardinia using geostatistical and climate models. Agricultural and Forest Entomology, 12: 417–426.

Duffera, M, Jeffrey, GW & Weisz, R 2007, Spatial variability of Southeastern U.S. Coastal Plain soil physical properties: Implications for site - specific manag -ement. Geoderma, 137: 327–339.

Erre, P, Chessa, I, Nieddu, G & Jones, PG 2009, Diversity and spatial distribution of Opuntia spp. in the Mediterranean Basin.  Journal of Arid Environments, 73: 1058–1066.

Fan, Z, Fan, X, Spetich, MA, Shifley, SR, Moser, WK, Jensen, RG & Kabrick, JM 2011, Developing a stand hazard index for oak decline in upland oak forests of the Ozark Highlands, Missouri. Northern Journal of Applied Forestry, 28: 19-26.

Foster, JR, Townsend, PA & Mladenoff, DJ 2013, Spatial dynamics of a gypsy moth defoliation outbreak and dependence on habitat characteristics. Landscape Ecology, 28: 1307-1320.

Gunnarsson, F, Holm, S, Holmgren, P & Thuresson, T 1998, On the potential of kriging for forest management planning. Scandinavian Journal of Forest Research, 13: 237- 245.

Habashi, H, Hosseini, SM, Mohammadi, J & Rahmani, R 2007, Application of geostatistics techniques in the study of soils in forest areas. Iranian Journal of Agricultural Sciences and Natural Resources, 14: 18-27 (In Persian).

Hassani-Pak, A 2006, Geostatistics. Tehran University Press, Tehran, Iran, 2th edition, 314 p.

Heidari, R 2006, Investigation of different methods of distance inventory in Zagros forests. PhD Dissertation, Department of Natural Resources, Tehran University, 116 p.

Hlásny, T, Vizi, L, Turčáni, M, Koreň, M, Kulla, L & Sitková, Z 2009, Geostatistical simulation of bark beetle infestation for forest protection purposes. Journal of Forest Science, 55: 518–525.

Hohn, ME, Liebhold, AM & Gribko, LS 1993, A geostatistical model for forecasting the spatial dynamics of defoliation caused by the gypsy moth, Lymantria dispar (Lepidoptera: Lymantriidae). Environme -ntal Entomology, 22: 1066–1077.

Holdenrieder, O, Pautasso, M, Weisberg, P & Lonsdale, D 2004, Tree diseases and landscape processes: The challenge of landscape pathology. Trends in Ecology and Evolution, 19: 446− 452.

Hosseini, A, Hosseini, SM, Rahmani, A & Azadfar, D 2014, Comparison between two oak stands (healthy and affected by oak decline) in respect to characteristics of competitive environments at Ilam province. Iranian Journal of Forest and Poplar Research, 21: 565-577 (In Persian).

Hosseini, V, Akhavan, R & Tahmasebi, M 2012, Effect of Pistachio (Pistacia atlantica) canopy on the spatial distribution of soil chemical characteristics (Case study: Sarvabad, Kurdistan Province, Iran). Iranian Journal of Forest, 4: 13-24 (In Persian).

Jafarnia, Sh & Akbarinia, M 2014, Investigation of spatial distribution of soil and water properties by use of geostatistical in Mangrove forest of Qeshm Island. Iranian Journal of Forest and Poplar Research, 22: 673-686 (In Persian).

Kazemi-Poshtmasari, H, Tahmasebi-Sarvestani, Z, Kamkar, B, Shataei, Sh & Sadeghi, S 2012, Evaluation of geostatistical methods for estimating and zoning macronutrients in agricultural lands of Golestan Province, Iran. Iranian Journal of Water and Soil Science, 22: 201-218 (In Persian).

Kint, V, Meirvenne, MV, Nachtergale, L, Geuden, G & Lust, N 2003, Spatial methods for quantifying forest stand structure development: comparison between nearest neighbor indices and variogram analysis. Forest science, 49: 36-49.

Koch, FH & Smith, WD 2008, Spatio-temporal analysis of Xyleborus glabratus (Coleoptera: Circulionidae: Scolytinae) invasion in Eastern US forests. Environmental Entomology, 37: 442-452.

Liebhold, AM, Rossie, RE & Kemp, WP 1993, Geostatistics and geographic information systems in applied insect ecology. Annual Review of Entomology, 38: 303–327.

López-Granados, F, Jurado-Expósito, M, Atenciano, S, García-Ferrer, A, Orden, MS & García-Torres, L 2002, Spatial variability of agricultural soil parameters in southern Spain. Plant and Soil, 246: 97-105.

Luo, W, Taylor, MC & Parker, SR 2008, A comparison of spatial interpolation methods to estimate continuous wind speed surfaces using irregularly distributed data from England and Wales. International Journal of Climatology, 28: 947-959.

Machiwal, D & Jha, MK 2015, Identifying sources of groundwater contamination in a hard-rock aquifer system using multivariate statistical analyses and GIS-based geostatistical modeling techniques. Journal of Hydrology: Regional Studies, Article in Press: 31 p.

Maranon, T, Ajbilou, R, Ojedaand, F & Arroya, J 1999, Biodiversity of woody species in oak woodland of southern Spain and northern Morocco. Forest Ecology and Management, 115: 147-156.

Mardikis, MG, Kalivas, DP & Kollias, VJ 2005, Comparison of interpolation methods for the prediction of reference evapotranspir –ation — an application in Greece. Water Resources Management, 19: 251-278.

Maroufi, S, Toranjeyan, A & Zare-Abyaneh, H 2009, Evaluation of geostatistical methods for estimating electrical conductivity and pH of stream drained water in Hamedan -Bahar Plain, Iran. Iranian Journal of Water and Soil Conservation, 16: 169-187 (In Persian).

McCoy, RM 2005, Field methods in remote sensing. The Guilford Press, New York, 177 p.

Mir-Abolfathi, M 2013, Outbreaks of charcoal disease on Quercus spp. and Zelkova carpinifolia trees in forests of Zagros and Alborz mountains in Iran. Iranian Journal of Plant Pathology, 49: 257-263 (In Persian).

Mir-Mousavi, SH, Mazidi, A & Khosravi, Y 2010, Determining the best geostatistical method for estimating the distribution of rainfall using GIS (Case study: Isfahan Province, Iran). Iranian Journal of Geographic Space, 10: 105-120 (In Persian).

Mirzaei, J 2005, Investigating the relationship between vegetation with topography and soil in northern forests of Ilam Province, Iran. MSc. thesis, Mazandaran University, Sari, Iran, 75 p.

Misaghi, F & Muhammad, K 2005, Zoning of rainfall data using classical statistics and geostatistics and compared to artificial neural networks. Iranian Journal of Agriculture, 29: 1-13 (In Persian).

Mohamadi, S,  Salajegheh, A, Mahdavi, M & Bagheri, R 2012, An investigation on spatial and temporal variations of groundwater level in Kerman plain using suitable geostatistical method (During a 10-year period). Iranian journal of Range and Desert Research, 19: 60-71 (In Persian).

Montes, F, Hernandez, MJ & Canellas, I 2005, A geostatistical approach to cork production sampling in Quercus suber forests. Canadian Journal of Forest Research, 35: 2787-2796.

Nalder, IA & Wein, RW 1998, Spatial interpolation of Climatic Normals Test of a new method in the Canadian Boreal Forest. Agricultural and Forest Meteorology, 92: 211–225.

Nourzadeh-Haddad, M, Mahdian, MH & Malakouti, MJ 2013, Efficiency comparison of some geostatistical methods for investigating spatial variability of micro nutrients in agricultural lands, Case study: Hamadan Province, Iran. Iranian Journal of Water and Soil Science, 23: 71-81 (In Persian).

Olsen, SR & Sommers, LE 1990, Phosphorous. In: Page AL, Method of soil analysis. Part 2. 2ndAgron Monoger, ASA, Madison, WI. 403- 431.

Omidvar, K, Ebrahimi, R & Rasti, F 2014, Evaluation of geo statistical techniques in drought intensity zoning (Case study: northwest and central regions of Iran). Iranian Journal of Irrigation and Water Engineering, 5: 30-44 (In Persian).

Pourreza, M, Hosseini, SM & Zohrevandi, AA 2012, Spatial variations of diameter of Pistacia atlantica (Desf.) trees in Zagros area (Case study: Pirkashan, Kermanshah).  Iranian Journal of Wood & Forest Science and Technology, 19: 1-20 (In Persian).

Ristaino, JB & Gumpertz, M 2000, New frontiers in the study of dispersal and spatial analysis of epidemics caused by species in the genus Phyophthora. Annual Review of Phytopathology, 38: 541− 576.

Robertson, GP, Coleman, DC, Bledsoe, CS & Solins, P 1999, Standard soil methods for long-term ecological research. Oxford University Press Inc., 480 p.

Sánchez, ME, Caetano, P, Ferraz, J & Trapero, A 2002, Phytophthora disease of Quercus ilex in Southwestern Spain. Forest Pathology. 32: 5-18.

Sarmadian, F & Taghizadeh, R 2008, Comparison of interpolation methods for mapping soil quality, Case study (farm of Agricultural Faculty, Tehran university). Iranian Journal of Soil and Water Research, 40: 157-165.

Thiessen, AH 1911, Precipitation averages for large areas. Monthly Weather Review, 39: 1082–1084.

Tuominen, S, Fish, S & Poso, S 2003, Combining remote sensing, data from earlier inventories and geostatistical interpolation in multi-source forest inventory. Canadian Journal of Forest Research, 33: 624- 634.

Vannini, A, Natili, G, Anselmi, N, Montaghiand, A & Vettraino, AM 2010, Distribution and gradient analysis of Ink disease in chestnut forests. Forest Pathology, 40: 73–86.

Wackernagel, H 2003, Multivariate geostatistics: An introduction with applications. 3rd ed. Springer Verlag, New York, 387 p.

Walkley, A & Black IA 1934, An examination of the Degtjareff method for determining soil organic matter and a proposed modif -ication of the chromic acid titration method. Soil Science, 34: 29-38.

Westerman, REL, 1990, Soil testing and plant analysis. Soil Science Society of America, Madison, Wisconsin, USA, 784 p.

WuJ, Norvell, WA & Welch, RM 2006, Kriging on highly skewed data for DTPA-extractable soil Zn with auxiliary information for pH and organic carbon. Geoderma, 134: 187–199.

Yamamoto, JK 2005, Correcting the smoothing effect of ordinary kriging estimates. Mathematical Geology, 37: 69-94.