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

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

Karami O.¹, Fallah A.*¹, Shataei S.H.², Latifi H.³

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

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 m²), 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 (R²). 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.

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 m² 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(x_i)  is the measured sample value at point i, Z(x_i+h) is the measured sample value at position i+h and Z₂(x_i) and Z₂(x_i+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 + C₀), range (R) and nugget semivariance (C₀). 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 x₀ 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 x_i and the i^th 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.

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 C₀/(C+C₀) 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.

Table 3. variogram and covariogram parameters and fit models.

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, R² = 0.52. In contrast, the GPI (average ME percentage = 0.22, average RMSE percentage = 43.41, R² = 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.

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.

                                       %ME: mean error percentage

                                       %RMSE: root mean square error percentage

                                       R²: coefficient of determination

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