# A Comparison of Bangladesh Climate Surfaces from the Geostatistical Point of View

Authors: Avit Kumar Bhowmik
Journal: ISRN Meteorology
 (2012)

---

## Abstract

This study analyses the degree and margin of differences among the surfaces of annual total precipitation in wet days (PRCPTOT) and the yearly maximum value of the daily maximum temperature (TXx) of Bangladesh, produced by thin plate spline (TPS), inverse distance weighting (IDW), ordinary kriging (OK), and universal kriging (UK) methods of spatial interpolation. From the surface differences, the maximum and minimum differences are observed between the surfaces produced by TPS and IDW, and OK and UK, respectively. The residual plots from cross-validation depict that IDW and OK methods mostly under predict and TPS and UK methods mostly overpredict the observed climate indices’ values. Both the tendency of methods’ over and underprediction and the surface-differences decrease with the increase in the number of available spatial point observations. Finally, two performance measures—the index of agreement (d) and the coefficient of variation of prediction (ρf)—imply that there is a little difference in the prediction ability of the four different methods. The performance of the spatial interpolation improves with the increase in the number of available spatial points, and eventually the predicted climate surfaces get similar. However, UK shows better interpolation performance in most of the years.


---

## Body

## 1. Introduction

The intrinsic responsibility of the professionals, who deal with climate, is to provide insights regarding climate variability at any place at any time. The crucial part of this responsibility is to predict those climate variables at those places and those times, where observations of the climate elements are not available [1]. This issue has become more critical when it has been proved by the climatologists that the global metric for climate is no longer useful because climate effects are felt locally, and they are region-specific [2]. From this point of view, special skills and knowledge are required to predict the most reliable value for the desired local climate information. As Tveito [1] presents “traditionally this is done by using observed values at neighboring stations which are then adjusted for representativity, terrain and other effects affecting the local climatology. Such estimates have usually been carried out as single point calculations, often including subjective considerations based on local knowledge and experience. Most of these estimates will not be consistently derived and they are thereby not reproducible and cannot be regarded as homogenous. They are therefore of limited value, for example, for advanced climate analysis.”Spatial interpolation methods as geostatistical prediction methods, with their inherent properties and applications, have successfully been implemented to combine different georeferenced climate variables and parameters in such a way that it is now possible to give consistently derived estimates at any place for a certain time [3–8]. This is also because of the fact that the spatial interpolation methods deal with the most important property of climate variables—they have a temporal extent along with a spatial extent [9]. Caruso and Quarta [10] classified these methods according to their fundamental hypotheses and mathematical properties and entitled them as “deterministic method, statistical method, stochastic simulation method, physical model simulation method and combined method.” The application and performance of the classified methods are solely dependent on the characteristics of study areas and algorithms and parameters used. Thus, obtaining a universally appropriate spatial interpolation method for climate variability analysis is impossible; rather locally an application oriented interpolation method is obtainable [11]. Additionally, this locally appropriate spatial interpolation method selection is subject to the qualitative and quantitative analysis of the local spatial data, their exploratory analysis, and different stages of trial and errors with the methods, which is commonly done from cross-validation results. More precisely, the result of the appropriate method needs to be further examined for their accuracy [1, 12]. Studies [7, 10, 13–16] have put utmost importance on comparing the performances of different spatial interpolation methods based on the cross-validation. Thus the methods’ performance in terms of their over- and underprediction is a subject to study in their application of creating high resolution climate surfaces.Among the geostatistical approaches of spatial interpolation developed in light of climate prediction [17–19], the deterministic and stochastic methods have turned out to be the most reliable methods for climate variability analysis. The stochastic methods are the complicated geostatistical approaches in contrast to the simpler deterministic methods, which have been designed for better prediction of the spatial pattern of climate phenomena. But recently it has been explored that for few observation datasets, complicated kriging methods (stochastic) do not show significantly greater predictive skills than simpler deterministic techniques, such as the inverse square distance method [20, 21]. Therefore, it is important to compare different climate surfaces created by the application of different deterministic and stochastic methods. Especially, the difference of the surfaces created by different methods and their overprediction and underprediction tendency require thorough analysis to decide on the most accurate spatial patterns of climate variability.The spatial interpolation techniques have been used for mapping the spatial patterns of climatic fields in several regions of the world, such as France [22], Germany [23], the UK [24], Italy [25], Mexico [26, 27], Portugal [28, 29], and the USA [30]. There have been a few studies conducted on the Bangladesh climate, based on the data from meteorological stations [2, 21, 31–34], but so far no study has been conducted in Bangladesh on comparing the performance and prediction ability of different spatial interpolation methods to produce climate surfaces. This kind of study is very important since many climate indices representing a wide variety of Asian climate aspects are already in the phase of implementation [35].In light of the discussion above, this study compares climate surfaces that have been produced by two deterministic and two stochastic spatial interpolation methods, by evaluating their surface differences, residuals, and predictability. For comparison, it creates the high-resolution continuous surfaces of two climate indices—PRCPTOT and TXx. The PRCPTOT characterizes the annual total precipitation in wet days, and the TXx corresponds to the yearly maximum value of the daily maximum temperature. The indices are recommended by the Joint Project Commission for Climatology/Climate Variability and Predictability (CLIVAR) and Joint WMO/IOC Technical Commission for Oceanography and Marine Meteorology Expert Team on climate change detection and indices. The objectives established to achieve the research goal of obtaining the most accurate spatially interpolated surfaces are (1)  to compile a rainfall and temperature dataset for Bangladesh and thereby compute annually defined climate indices and prepare spatially interpolated surfaces with two alternative deterministic methods—thin plate spline (TPS) and inverse distance weighting (IDW), and with two alternative stochastic methods—ordinary kriging (OK) and universal Kriging (UK); (2) to compute the surface differences between the prepared climate surfaces and analyze the degree and margin of differences in the surfaces; (3) to create the residual plots from cross-validation of the different methods and analyze the tendency of overprediction and underprediction by the methods; (4) to evaluate the interpolation cross-validation results by using two statistical performance measurements—index of agreement (d) and coefficient of variation of prediction (ρf).


## 2. Study Area, Dataset, and Spatial Point  Configuration

Bangladesh, situated in south-east Asia, is one of the most vulnerable countries of the world regarding the adverse impacts of anthropogenic climate change [2, 31–34] (Figure 1). The total area of the country is 147,570 square kilometers [36], approximately one fifth of which consists of low-lying coastal zones within one meter of the high water mark [37]. The Bangladesh Meteorological Department (BMD) is the authorized government organization for all meteorological activities in the country. It measures the daily precipitation and daily temperature with 34 meteorological stations situated in different locations all over the country [38] (Figure 1).Figure 1
Study area—Bangladesh with the spatial point configuration of the 34 meteorological stations from which the data is available in 2007.The dataset used in this study includes 60 years’, that is, 1948–2007, daily precipitation and temperature measurements from the meteorological stations of BMD. The measurements are not available from all stations, that is, spatial points, for the whole study period, that is, temporal points. As described in Figure2, the spatial point configuration of the dataset therefore is different from year to year (temporal point to temporal point). Precipitation and temperature measurements from 10 stations are available for 1948 which gradually increase to the measurements from 34 meteorological stations by 2007. Eventually the average lag distance between the spatial points available in 1948 is 100 km which gradually decreases to 50 km by 2007. Number and distribution of available spatial points are important aspects to consider in spatial interpolation [1, 23, 39]. The variogram analysis and distance decay principle of spatial interpolation need special consideration addressing the spatial point configuration [7, 8, 23, 40] (Figure 2), which is applicable for this study in spite of the irregular spatial point configuration.Figure 2
Spatial point configuration of Bangladesh over the study period of 1948–2007.


## 3. Methodological Framework

The open-source statistical software “R” [41] with its dedicated packages—“spacetime,” “intamap,” “fields,” “gstat,” “maptools,” and “sm”—has been used in this study for the analyses and computation. The following methodology has been followed to obtain the research objectives.


### 3.1. The Climate Indices

Two climate indices—PRCPTOT and TXx—[42–45] have been calculated from the available precipitation and temperature measurements, for each of the years of 1948–2007 for each station. The unavailable measurements from a station for a particular year have been considered as “NA value” in indices calculation.PRCPTOT, the annual total precipitation in wet days, has been calculated by summing the daily precipitation in the wet days of the years [42, 45]. Bangladesh has clearly defined wet days in every year; it is known as “Monsoon”, and the duration is June-September of every year [31, 37, 46–48]. Thus the measured precipitation during June-September has been summed up to calculate the PRCPTOT of each year.TXx, the yearly maximum value of the daily maximum temperature, has been calculated by picking the maximum value of the temperatures available from the daily measurements of each year [42, 43]. Previous studies have proven that the change in temperature of Bangladesh due to climate change is recognizable from the trend of the maximum temperatures [31, 32]. Thus the maximum temperature of each of the years 1948–2007 has been recorded as TXx values for spatial interpolation.


### 3.2. Spatial Interpolation of Climate Indices

In practice, the spatial interpolation methods that incorporate temporal information with spatial information in the modeling function are most appropriate for interpolating climate variables [49, 50]. The structure of a basic spatial interpolation problem denotes that the dependent variable of interest Z(s0,t) is predicted as an output of the mathematical function of known predictors (si,t)  (i=0,…,M), where the location vectors s are the elements of the given space domain D and t is time. The vector form of the predictors is  ZT(t)=[Z(s1,t),…,Z(sM,t)  [18]. The probability distribution of the climate variables sets up the appropriate interpolation formulae, which include some unknown interpolation parameters. These parameters can be obtained through known functions of certain statistical parameters. Modeling of climate variables with these statistical parameters assumes that the expected values of the variables are changing in space and in time in a similar way [1, 51]. The spatial change in the variables indicates that the climate is different in the regions, whereas temporal change indicates the global climate variability with time [52]. As a result, expected values of climate variables can be obtained by the following linear model of (1) [53, 54]:
(1)E(Z(si,t))=μ(t)+E(si)(i=0,…,M),
where μ(t) is the temporal trend or the climate change signal and E(s) is the spatial trend. Typically, there is only a single realization in time for the modeling of the statistical parameters in spatial interpolation. Therefore only the predictors Z(si,t)(i=0,…,M) constitute the usable information or the sample for the modeling of variability over space [18]. At the linear model (1), the basic statistical parameters can be allocated into two categories—deterministic and stochastic parameters. Thus the spatial interpolation methods can be divided into two groups—deterministic and stochastic spatial interpolation methods [18, 51].


#### 3.2.1. Deterministic Spatial Interpolation

Two deterministic approaches based on linear equation (1) have been fitted as spatial interpolation models—the inverse distance weighting (IDW) approach predicts the dependent variable based on the extent of similarity, whereas the thin plate spline (TPS) approach predicts it based on the degree of smoothing [6]. IDW predicts spatially continuous long range variation by giving more weight to nearby measurements than to distant measurements, and TPS predicts the variation by a suitably continuous smoothing spline function. There are distinguishing features of these two methods, which characterize their performances in estimating continuous surfaces. TPS considers the spatial dependence of the phenomenon as a smoothed spline from the GCV function (Figure 3(a)), whereas IDW considers the spatial dependence as straight lines connecting the measured points without any smoothing (Figure 3(b)) [10]. The nonsmoothing property of IDW sometimes enables it to result in better cross-validation results, since it takes the actual measurement as variability factors. Yet, it is also necessary to consider that it includes the measurement errors into the estimated values, which are unknown in most cases [55]. Smoothing property and taking the long range and short range variability into consideration enables TPS to perform better for the few observation point regions like Bangladesh.Inherent interpolation models of (a) inverse distance weighting and (b) thin plate spline methods.
(a)
(b)


#### 3.2.2. Stochastic Spatial Interpolation

Stochastic spatial interpolation methods exploit the statistical properties of the measured points and quantify the spatial autocorrelation among the measured points. Furthermore, they account for the spatial configuration of the sample points around the prediction location [51, 54, 56]. These methods are based on the covariance or the variance between the dependent variables and the predictors. The principle of stochastic methods for which they are distinguished from the deterministic methods is that they do not design the variography exactly through the measured points. Rather they describe a smoothed variogram over the area of interest and thus try to fit a smoothed line, which describes the continuous variables’ behavior over the region of interest [51] (Figure 4). Two stochastic models based on linear equation (1) have been fitted. The first approach is ordinary kriging (OK), which accounts for any scale-dependent trend present in the indices behavior through the local system and also considers anisotropy in the indices interpolation. The final approach is the universal kriging (UK) which explicitly models the trend component as a linear combination of functions of the spatial coordinates [16].Figure 4
Inherent interpolation model based on variography of ordinary kriging and universal kriging methods.


#### 3.2.3. Variogram for Stochastic Spatial Interpolation

As discussed before, in case of irregular spatial point configuration and highly increased average lag distance between the spatial points, the covariance function may not be defined, and the number of pairs of measurements separated by the lag distance may not be enough to analyze the variogram. The temporal measurements in that case can be used to increase the number of pairs of measurements within the lag distance for variogram analysis. This approach is known as the “pooled variogram” analysis discussed by Bivand et al. [57] and also supported by the idea of “increasing features” in analysis by Raudys and Jain [58]. In pooled variogram analysis the temporal domain is included by stacking the temporal measurements on top of each other. The vertical (temporal) distance from pooled measurements is scaled by “dX parameter” [57] to avoid the temporal influence in the stacks so that a mean variogram can be created out of the stack. Many geostatistical studies have applied pooled variogram approach in case of irregular spatial point configuration over temporal points with increased lag [59, 60]. 12 mean variograms have been created following the pooled variogram approach for the temporal periods of 1948–1975, 1976–1992, and 1993–2007 to create the stochastically interpolated climate surfaces of PRCPTOT and TXx. The temporal periods for creating the pooled variogram have been selected based on the availability of the number of spatial points in a way, so that every temporal period has more or less the same number of total spatial points.The mathematical formulations of the four spatial interpolation methods and the pooled variograms are available in the provided literature. The formulae have been adopted in R and applied to produce the climate surfaces. The surfaces of the climate indices have been interpolated in 15 m× 15 m continuous grids distributed over the spatial extent of the study area. The grid size has been determined in a way so that it provides the highest possible spatial resolution in climate indices’ prediction and at the same time is competent enough according to the computation efficiency of the tools used.


### 3.3. Surface-Differences

The surface-differences between the climate surfaces produced by the four different interpolation methods have been calculated and created in R using the raster-subtraction method. Thus the differences between the produced surfaces by TPS and IDW, TPS and OK, TPS and UK, IDW and OK, IDW and UK, and OK and UK have been computed.


### 3.4. Residual Plots

The residual plots have been created to analyze the tendency of overprediction and underprediction by the spatial interpolation methods using the “fictitious-point method,” which is popularly known as “cross-validation.” It was first proposed by Seaman [56] and subsequently applied by geostatisticians in a whole range of studies for evaluating the performances of the spatial interpolation techniques [10, 20, 49–54, 61, 62]. Validation is a statistical method of evaluating and comparing learning algorithms by dividing data into two segments: one segment is used to learn or train a model and the other one is used to validate the model. In typical cross-validation, the training and validation sets must cross over in successive rounds so that each data point has a chance of being validated against. The basic form of cross-validation is k-fold cross-validation [63]. In geostatistics, generally, the training set is the measured values of all sample points, which are validated using the appointed spatial interpolation technique. It is typically known as “leave-one-out” cross-validation, since each sampled point is taken out successively and then predicted using the remaining measured points of the same sample set [57]. An ideal cross-validation plot should give the points plotted along an equidistant line from the axes defined by the measured and predicted values. In practice, the predicted values differ from the observed values and the difference in known as the residuals. If the predicted value is higher than the measured value, it is overpredicted, and if the predicted value is lower than the measured value, it is under predicted. Thus, the residuals determine the tendency of over- and underprediction by the methods and therefore the quality of prediction [64]. The residual plots for each of the years of 1948–2007 of each interpolation method have been created by plotting the residuals applying “bubble plot” in the gstat package of R and then compared.


### 3.5. Performance Measures

Two performance measures have been computed from cross-validation to analyze the predictability of the spatial interpolation methods in describing climate variability—the index of agreement (d) [65] and the coefficient of variation of prediction (ρf) [64]. d  has been computed by the formulae of (2):
(2)d=1-N*RMSE2PE,
where
(3)PE=∑i=1N(|Pi-O-|+|Oi-O-|)2,(4)RMSE=[N-1∑i=1N(Pi-Oi)2]1/2.Pi and Oi are predicted and observed indices’ values at ith location, and O- is the mean of the Oi s. d varies between 0.0 and 1.0. 1.0 conveys the perfect agreement and thus the better performance of the spatial interpolation technique in prediction, and 0.0 conveys complete disagreement [20]. The root mean square error RMSE, is the error measure commonly computed in geographic applications, which computes the distribution of the predicted and observed values.The coefficient of variation of prediction (ρf) has been computed by (5):
(5)ρf=RMSE(1/n)∑i=inOi.ρf is expressed in percentage. It is the measure of spatial uncertainty of prediction [64] and is also used to compare the performance of interpolation schemes for different integration times. Thus, the lower is the percentage of ρf  resulted in the cross-validation of a spatial interpolation technique, the better is the performance of that technique and vice versa.


## 3.1. The Climate Indices

Two climate indices—PRCPTOT and TXx—[42–45] have been calculated from the available precipitation and temperature measurements, for each of the years of 1948–2007 for each station. The unavailable measurements from a station for a particular year have been considered as “NA value” in indices calculation.PRCPTOT, the annual total precipitation in wet days, has been calculated by summing the daily precipitation in the wet days of the years [42, 45]. Bangladesh has clearly defined wet days in every year; it is known as “Monsoon”, and the duration is June-September of every year [31, 37, 46–48]. Thus the measured precipitation during June-September has been summed up to calculate the PRCPTOT of each year.TXx, the yearly maximum value of the daily maximum temperature, has been calculated by picking the maximum value of the temperatures available from the daily measurements of each year [42, 43]. Previous studies have proven that the change in temperature of Bangladesh due to climate change is recognizable from the trend of the maximum temperatures [31, 32]. Thus the maximum temperature of each of the years 1948–2007 has been recorded as TXx values for spatial interpolation.


## 3.2. Spatial Interpolation of Climate Indices

In practice, the spatial interpolation methods that incorporate temporal information with spatial information in the modeling function are most appropriate for interpolating climate variables [49, 50]. The structure of a basic spatial interpolation problem denotes that the dependent variable of interest Z(s0,t) is predicted as an output of the mathematical function of known predictors (si,t)  (i=0,…,M), where the location vectors s are the elements of the given space domain D and t is time. The vector form of the predictors is  ZT(t)=[Z(s1,t),…,Z(sM,t)  [18]. The probability distribution of the climate variables sets up the appropriate interpolation formulae, which include some unknown interpolation parameters. These parameters can be obtained through known functions of certain statistical parameters. Modeling of climate variables with these statistical parameters assumes that the expected values of the variables are changing in space and in time in a similar way [1, 51]. The spatial change in the variables indicates that the climate is different in the regions, whereas temporal change indicates the global climate variability with time [52]. As a result, expected values of climate variables can be obtained by the following linear model of (1) [53, 54]:
(1)E(Z(si,t))=μ(t)+E(si)(i=0,…,M),
where μ(t) is the temporal trend or the climate change signal and E(s) is the spatial trend. Typically, there is only a single realization in time for the modeling of the statistical parameters in spatial interpolation. Therefore only the predictors Z(si,t)(i=0,…,M) constitute the usable information or the sample for the modeling of variability over space [18]. At the linear model (1), the basic statistical parameters can be allocated into two categories—deterministic and stochastic parameters. Thus the spatial interpolation methods can be divided into two groups—deterministic and stochastic spatial interpolation methods [18, 51].


### 3.2.1. Deterministic Spatial Interpolation

Two deterministic approaches based on linear equation (1) have been fitted as spatial interpolation models—the inverse distance weighting (IDW) approach predicts the dependent variable based on the extent of similarity, whereas the thin plate spline (TPS) approach predicts it based on the degree of smoothing [6]. IDW predicts spatially continuous long range variation by giving more weight to nearby measurements than to distant measurements, and TPS predicts the variation by a suitably continuous smoothing spline function. There are distinguishing features of these two methods, which characterize their performances in estimating continuous surfaces. TPS considers the spatial dependence of the phenomenon as a smoothed spline from the GCV function (Figure 3(a)), whereas IDW considers the spatial dependence as straight lines connecting the measured points without any smoothing (Figure 3(b)) [10]. The nonsmoothing property of IDW sometimes enables it to result in better cross-validation results, since it takes the actual measurement as variability factors. Yet, it is also necessary to consider that it includes the measurement errors into the estimated values, which are unknown in most cases [55]. Smoothing property and taking the long range and short range variability into consideration enables TPS to perform better for the few observation point regions like Bangladesh.Inherent interpolation models of (a) inverse distance weighting and (b) thin plate spline methods.
(a)
(b)


### 3.2.2. Stochastic Spatial Interpolation

Stochastic spatial interpolation methods exploit the statistical properties of the measured points and quantify the spatial autocorrelation among the measured points. Furthermore, they account for the spatial configuration of the sample points around the prediction location [51, 54, 56]. These methods are based on the covariance or the variance between the dependent variables and the predictors. The principle of stochastic methods for which they are distinguished from the deterministic methods is that they do not design the variography exactly through the measured points. Rather they describe a smoothed variogram over the area of interest and thus try to fit a smoothed line, which describes the continuous variables’ behavior over the region of interest [51] (Figure 4). Two stochastic models based on linear equation (1) have been fitted. The first approach is ordinary kriging (OK), which accounts for any scale-dependent trend present in the indices behavior through the local system and also considers anisotropy in the indices interpolation. The final approach is the universal kriging (UK) which explicitly models the trend component as a linear combination of functions of the spatial coordinates [16].Figure 4
Inherent interpolation model based on variography of ordinary kriging and universal kriging methods.


### 3.2.3. Variogram for Stochastic Spatial Interpolation

As discussed before, in case of irregular spatial point configuration and highly increased average lag distance between the spatial points, the covariance function may not be defined, and the number of pairs of measurements separated by the lag distance may not be enough to analyze the variogram. The temporal measurements in that case can be used to increase the number of pairs of measurements within the lag distance for variogram analysis. This approach is known as the “pooled variogram” analysis discussed by Bivand et al. [57] and also supported by the idea of “increasing features” in analysis by Raudys and Jain [58]. In pooled variogram analysis the temporal domain is included by stacking the temporal measurements on top of each other. The vertical (temporal) distance from pooled measurements is scaled by “dX parameter” [57] to avoid the temporal influence in the stacks so that a mean variogram can be created out of the stack. Many geostatistical studies have applied pooled variogram approach in case of irregular spatial point configuration over temporal points with increased lag [59, 60]. 12 mean variograms have been created following the pooled variogram approach for the temporal periods of 1948–1975, 1976–1992, and 1993–2007 to create the stochastically interpolated climate surfaces of PRCPTOT and TXx. The temporal periods for creating the pooled variogram have been selected based on the availability of the number of spatial points in a way, so that every temporal period has more or less the same number of total spatial points.The mathematical formulations of the four spatial interpolation methods and the pooled variograms are available in the provided literature. The formulae have been adopted in R and applied to produce the climate surfaces. The surfaces of the climate indices have been interpolated in 15 m× 15 m continuous grids distributed over the spatial extent of the study area. The grid size has been determined in a way so that it provides the highest possible spatial resolution in climate indices’ prediction and at the same time is competent enough according to the computation efficiency of the tools used.


## 3.2.1. Deterministic Spatial Interpolation

Two deterministic approaches based on linear equation (1) have been fitted as spatial interpolation models—the inverse distance weighting (IDW) approach predicts the dependent variable based on the extent of similarity, whereas the thin plate spline (TPS) approach predicts it based on the degree of smoothing [6]. IDW predicts spatially continuous long range variation by giving more weight to nearby measurements than to distant measurements, and TPS predicts the variation by a suitably continuous smoothing spline function. There are distinguishing features of these two methods, which characterize their performances in estimating continuous surfaces. TPS considers the spatial dependence of the phenomenon as a smoothed spline from the GCV function (Figure 3(a)), whereas IDW considers the spatial dependence as straight lines connecting the measured points without any smoothing (Figure 3(b)) [10]. The nonsmoothing property of IDW sometimes enables it to result in better cross-validation results, since it takes the actual measurement as variability factors. Yet, it is also necessary to consider that it includes the measurement errors into the estimated values, which are unknown in most cases [55]. Smoothing property and taking the long range and short range variability into consideration enables TPS to perform better for the few observation point regions like Bangladesh.Inherent interpolation models of (a) inverse distance weighting and (b) thin plate spline methods.
(a)
(b)


## 3.2.2. Stochastic Spatial Interpolation

Stochastic spatial interpolation methods exploit the statistical properties of the measured points and quantify the spatial autocorrelation among the measured points. Furthermore, they account for the spatial configuration of the sample points around the prediction location [51, 54, 56]. These methods are based on the covariance or the variance between the dependent variables and the predictors. The principle of stochastic methods for which they are distinguished from the deterministic methods is that they do not design the variography exactly through the measured points. Rather they describe a smoothed variogram over the area of interest and thus try to fit a smoothed line, which describes the continuous variables’ behavior over the region of interest [51] (Figure 4). Two stochastic models based on linear equation (1) have been fitted. The first approach is ordinary kriging (OK), which accounts for any scale-dependent trend present in the indices behavior through the local system and also considers anisotropy in the indices interpolation. The final approach is the universal kriging (UK) which explicitly models the trend component as a linear combination of functions of the spatial coordinates [16].Figure 4
Inherent interpolation model based on variography of ordinary kriging and universal kriging methods.


## 3.2.3. Variogram for Stochastic Spatial Interpolation

As discussed before, in case of irregular spatial point configuration and highly increased average lag distance between the spatial points, the covariance function may not be defined, and the number of pairs of measurements separated by the lag distance may not be enough to analyze the variogram. The temporal measurements in that case can be used to increase the number of pairs of measurements within the lag distance for variogram analysis. This approach is known as the “pooled variogram” analysis discussed by Bivand et al. [57] and also supported by the idea of “increasing features” in analysis by Raudys and Jain [58]. In pooled variogram analysis the temporal domain is included by stacking the temporal measurements on top of each other. The vertical (temporal) distance from pooled measurements is scaled by “dX parameter” [57] to avoid the temporal influence in the stacks so that a mean variogram can be created out of the stack. Many geostatistical studies have applied pooled variogram approach in case of irregular spatial point configuration over temporal points with increased lag [59, 60]. 12 mean variograms have been created following the pooled variogram approach for the temporal periods of 1948–1975, 1976–1992, and 1993–2007 to create the stochastically interpolated climate surfaces of PRCPTOT and TXx. The temporal periods for creating the pooled variogram have been selected based on the availability of the number of spatial points in a way, so that every temporal period has more or less the same number of total spatial points.The mathematical formulations of the four spatial interpolation methods and the pooled variograms are available in the provided literature. The formulae have been adopted in R and applied to produce the climate surfaces. The surfaces of the climate indices have been interpolated in 15 m× 15 m continuous grids distributed over the spatial extent of the study area. The grid size has been determined in a way so that it provides the highest possible spatial resolution in climate indices’ prediction and at the same time is competent enough according to the computation efficiency of the tools used.


## 3.3. Surface-Differences

The surface-differences between the climate surfaces produced by the four different interpolation methods have been calculated and created in R using the raster-subtraction method. Thus the differences between the produced surfaces by TPS and IDW, TPS and OK, TPS and UK, IDW and OK, IDW and UK, and OK and UK have been computed.


## 3.4. Residual Plots

The residual plots have been created to analyze the tendency of overprediction and underprediction by the spatial interpolation methods using the “fictitious-point method,” which is popularly known as “cross-validation.” It was first proposed by Seaman [56] and subsequently applied by geostatisticians in a whole range of studies for evaluating the performances of the spatial interpolation techniques [10, 20, 49–54, 61, 62]. Validation is a statistical method of evaluating and comparing learning algorithms by dividing data into two segments: one segment is used to learn or train a model and the other one is used to validate the model. In typical cross-validation, the training and validation sets must cross over in successive rounds so that each data point has a chance of being validated against. The basic form of cross-validation is k-fold cross-validation [63]. In geostatistics, generally, the training set is the measured values of all sample points, which are validated using the appointed spatial interpolation technique. It is typically known as “leave-one-out” cross-validation, since each sampled point is taken out successively and then predicted using the remaining measured points of the same sample set [57]. An ideal cross-validation plot should give the points plotted along an equidistant line from the axes defined by the measured and predicted values. In practice, the predicted values differ from the observed values and the difference in known as the residuals. If the predicted value is higher than the measured value, it is overpredicted, and if the predicted value is lower than the measured value, it is under predicted. Thus, the residuals determine the tendency of over- and underprediction by the methods and therefore the quality of prediction [64]. The residual plots for each of the years of 1948–2007 of each interpolation method have been created by plotting the residuals applying “bubble plot” in the gstat package of R and then compared.


## 3.5. Performance Measures

Two performance measures have been computed from cross-validation to analyze the predictability of the spatial interpolation methods in describing climate variability—the index of agreement (d) [65] and the coefficient of variation of prediction (ρf) [64]. d  has been computed by the formulae of (2):
(2)d=1-N*RMSE2PE,
where
(3)PE=∑i=1N(|Pi-O-|+|Oi-O-|)2,(4)RMSE=[N-1∑i=1N(Pi-Oi)2]1/2.Pi and Oi are predicted and observed indices’ values at ith location, and O- is the mean of the Oi s. d varies between 0.0 and 1.0. 1.0 conveys the perfect agreement and thus the better performance of the spatial interpolation technique in prediction, and 0.0 conveys complete disagreement [20]. The root mean square error RMSE, is the error measure commonly computed in geographic applications, which computes the distribution of the predicted and observed values.The coefficient of variation of prediction (ρf) has been computed by (5):
(5)ρf=RMSE(1/n)∑i=inOi.ρf is expressed in percentage. It is the measure of spatial uncertainty of prediction [64] and is also used to compare the performance of interpolation schemes for different integration times. Thus, the lower is the percentage of ρf  resulted in the cross-validation of a spatial interpolation technique, the better is the performance of that technique and vice versa.


## 4. Results and Discussions

### 4.1. Justifying the Choice of Spatial Interpolation Methods

The calculated climate indices—PRCPTOT and TXx—show spatial trend and anisotropic behavior over the study area. The correlation of PRCPTOT with longitude (0.55) is higher than the correlation of PRCPTOT with latitude (−0.42), which indicates that the spatial trend is more dominant in the west-east direction than the north-south direction. On the other hand, the correlation of TXx with longitude (−0.52) is again higher than the correlation of TXx with latitude (0.34), which indicates that the spatial trend is more dominant in the west-east direction than the north-south direction. The anisotropic behavior of the indices is not very significant since the ratio of the major and minor axes varies from 0.7 to 1.0 (Table1), but it exists in the indices behavior. The results of the analysis of the temporal trend of the calculated climate change indices indicate that for PRCPTOT the range of the trend value is −18 to 36 for 60 years. Most of the stations show an increasing trend of PRCPTOT over time. Especially the stations in the south-east region of the country experience the highest increasing trend, where the highest values of PRCPTOT are also experienced. On the other hand, for TXx a range of −0.09 to 0.33 for the trend value has been obtained. Though most of the stations show an increasing trend of TXx, almost all the stations in the midregion, including capital Dhaka, show a decreasing trend.Table 1
Variogram models and parameters for OK interpolation of PRCPTOT and TXx.
Index
Temporal period
Variography
Model
Sill
Range
Nugget
Anisotropy parameters
Angle
Ratio of major and minor axis
PRCPTOT
1948–1975
Spherical
700000
400
45000
2.96
0.83
1976–1992
Spherical
630000
550
205000
47.28
0.97
1993–2007
Spherical
780000
550
70000
83.94
0.87
TXx
1948–1975
Spherical
15
430
0
2.96
0.71
1976–1992
Spherical
10
530
0
0.96
0.71
1993–2007
Spherical
5.8
410
1
176.72
0.73As  Isaaks and Srivastava [5], and Atkinson and Tate [66] discovered, that in case of uncertain spatial point configuration, deterministic methods should be the preferable option for spatial interpolation, two deterministic methods—thin plate spline (TPS) and inverse distance weighting (IDW)—have been applied for interpolation of the indices. These methods have mostly been applied in climate modeling by the geostatisticians. Based on the very low temporal trend of the indices, the pooled variogram approach is justified since a mean variogram of several temporal points should describe the spatial variability of the indices more accurately. Finally, after modeling the pooled variograms, ordinary kriging (OK) has been applied taking the anisotropic behavior of the indices into account, and universal kriging (UK) has been applied by taking the discussed spatial trend of the indices into account.


### 4.2. The Pooled Variograms

12 pooled variograms have been designed for modeling OK and UK for interpolating PRCPTOT and TXx. The dX parameter has been set to 0.55 to rescale the temporal distance of 1 year to the spatial distance of 550 km (the maximum distance obtained in the spatial point pairs) within the stacks. The total numbers of spatial point collected in each stack of the 1948–1975, 1976–1992, and 1993–2007 are 483, 494, and 503, respectively. The variogram parameters have been presented in Tables1 and 2. The anisotropy has been analyzed and modeled in the variograms of OK. The fundamental principle behind designing the UK variograms is that both PRCPTOT and TXx indices behavior over the study region are functions of longitude and latitude, which indicates that there is a spatial trend in the indices behavior. The UK variograms for interpolating PRCPTOT and TXx result in higher nugget values than the OK variograms, which represents that the spatial trend parameter models the long range spatial variability of the indices and therefore results in abrupt changes in the short range variability.Table 2
Variogram models and parameters for UK interpolation of PRCPTOT and TXx.
Index
Temporal period
Dependent parameters
Variography
Model
Sill
Range
Nugget
PRCPTOT
1948–1975
Longitude and Latitude
Spherical
300000
400
115000
1976–1992
Longitude and Latitude
Spherical
200000
550
260000
1993–2007
Longitude and Latitude
Spherical
270000
500
136000
TXx
1948–1975
Longitude and Latitude
Spherical
10
420
1
1976–1992
Longitude and Latitude
Spherical
2.1
480
1
1993–2007
Longitude and Latitude
Spherical
3.9
450
1.4OK variograms for TXx result in no nugget for the temporal periods of 1948–1975 and 1976–1992, whereas for 1993–2007 there is a nugget value of 1, which is high for this temporal period. The OK variogram for TXx of the temporal period of 1993–2007 also represents a lower sill value, while the OK variogram for PRCPTOT of the 1976–1992 period represents a higher sill value. These values show that the TXx values vary to a lower spatial extent during 1993–2007, while the PRCPTOT values vary to a higher spatial extent during 1976–1992 over the study region. The UK variograms for TXx of the temporal periods of 1976–1992 and 1993–2007 represent lower sill values, which indicate that TXx values did not vary to a great extent in the direction of spatial trend during the mentioned temporal periods. The ranges of all variograms are similar to the extent of the whole study region, which represents that the spatial variability of the indices is more accurately describable in the long range of the whole study area than in the short range of the neighboring stations.


### 4.3. The Interpolated Surfaces of the Climate Indices

The inherent methodologies of the four different spatial interpolation methods have been reflected in their produced surfaces of the climate indices. The search neighborhood [67, 68] has been set to the maximum number of available station observations in the corresponding years. The produced climate surfaces have been represented in Figures 5 and 6. To ease the comparison and visualization, surfaces from four different years have been selected as representatives of the three temporal periods for which the pooled variograms have been computed—1949, 1975, 1987, and 2006. And the following difference surfaces and residual plots will also be represented for these four different years to provide the readers with better understanding of the result. The PRCPTOT values range from 0 to 4500 mm, and TXx values range from 20° to 65°C.Figure 5
Representative and comparative PRCPTOT surfaces of different years produced by the four different spatial interpolation methods.Figure 6
Representative and comparative TXx surfaces of different years produced by the four different spatial interpolation methods.The effect of the property of the global information carried by the stationary mean in case of few observations [12] has been represented by the IDW surfaces for both PRCPTOT and TXx indices. Instead of representing the long range variability of the indices, IDW represents more variability close to the sample location in a circular manner and shows almost no variation in distance. This is identified by the “Bull’s eyes” that appeared in the surfaces in every station location, which are expected from the modeling features of IDW described in Section 3.2.1. Consequently, no clear spatial trend has been recognized from the surfaces, but rather some higher values have been recognized close to the areas of the representative stations.TPS surfaces, as a result of the inherent smoothing properties, are smoothed and present the long range variability of the phenomena. From the PRCPTOT surfaces, the spatial trend of the index is apparent, in all years the highest PRCPTOT has been experienced in the south-eastern part of the country. It is hard to predict any temporal trend in the surfaces for PRCPTOT. As described in the beginning of Section4, the spatial trend of TXx is opposite to the one of PRCPTOT, which has also been represented by the spline surfaces. The surfaces represent a spatial pattern of higher TXx in the north-western part of the country, which is dominant in almost all the years.OK method has been applied using the variography described in Section4.2. The surfaces describe the indices behavior significantly—both in terms of scale-dependent trend and anisotropy. Again it is difficult to predict any temporal trend from the surfaces. Like PRCPTOT surfaces, TXx surfaces produced by ordinary kriging do not describe any temporal trend in particular, but the spatial trend to the north-western part of the country is clearly represented by them.The tilted surfaces of UK have been created using the variograms designed in Section4.2. The basic principle, which has been applied to create these surfaces, is the fact that there is a spatial trend in indices behavior, which is also clearly visible in the created surfaces. The spatial trend is more dominant in the east-west direction than the north-south as expected from the correlation with longitude and latitude.From the climate surfaces presented in Figures5 and 6, a clear difference is visible in the surfaces in terms of the spatial pattern of climate variability. The difference occurs because of the different principles followed by different spatial interpolation methods apart from the fact that they interpolate the same observations from the stations.


### 4.4. The Surface-Differences of the Climate Indices

To analyze the degree and the margin of differences between the interpolated climate surfaces of PRCPTOT and TXx generated through the four different methods, the surface-difference has been computed. The computed difference surfaces are presented in Figures7 and 8. The surface-differences for both PRCPTOT and TXx show significant differences in the spatial pattern of the climate variability interpolated by different spatial interpolation methods. The range of differences observed among the PRCPTOT surfaces is 1000 mm to −1000 mm and among the TXx surfaces is 8°C to −8°C. The variation among the PRCPTOT surfaces is more diverse than the variation among the TXx surfaces and has therefore been reflected in their difference surfaces. This is a result of the difference between the spatial pattern of PRCPTOT and TXx; precipitation is not as continuous in time and space as temperature. In most cases, precipitation field represents a bounded continuum in its spatial and temporal pattern, and therefore the surfaces induce more noise.Representative difference surfaces showing the difference between the PRCPTOT surfaces produced by the four different spatial interpolation methods.Representative difference surfaces showing the difference between the PRCPTOT surfaces produced by the four different spatial interpolation methods.The maximum difference has been detected between the surfaces created by the two deterministic methods (TPS and IDW) for both PRCPTOT and TXx. The reason of the maximum difference between the deterministic surfaces is clear from their spatial modeling concept described in Section3.2.1. The dispute between the smoothing and nonsmoothing properties and therefore modeling short range and long range variability cause significant differences in the surfaces. The minimum difference has been revealed between the surfaces created by the two stochastic methods (OK and UK) for both PRCPTOT and TXx. The stochastic surfaces have minimum differences because the interpolators are identical and characterize the spatial variability through the variogram in similar ways. Thus, the stochastic methods are similar in their interpolation principle but model the variables based on the variable behavior (anisotropy and spatial trend) over space.The IDW method is unique among the methods as it creates the bull’s eyes in the surfaces, and therefore in the computed difference-surfaces between the IDW and others, significant differences have been observed in the areas within the close buffer of the station locations. The most variation in the differences has also been occurred in the areas nearby station locations. Therefore, the tendency of the IDW method being restricted to interpolate short range variability and including measurement errors is apparent from this analysis.The surfaces are more different in the areas with extreme indices values like TPS-IDW and TPS-OK surfaces of PRCPTOT of 1987. In the northern and eastern parts of the country, the differences in these surfaces exceed 1000 mm (Figure7). This indicates that different spatial interpolation methods predict the climate variables more differently when the climate behavior deviates than the trend. Another important fact is that the differences among the surfaces have been decreasing over years with the increasing number of spatial points and decreasing spatial lag, which indicates that the performance of the different spatial interpolation techniques becomes similar if there are more point observation available. This is an important property to be considered in decision making on the part of spatial interpolation, because the decision is more complicated where there are few point observationss available. If enough point observations are available, the interpolated surfaces by different spatial interpolation methods are mostly similar.


### 4.5. The Residual Plots from Cross-Validation

The residual plots of prediction have been created from the cross-validation of the four different spatial interpolation methods for creating the surfaces of PRCPTOT and TXx. The residual plots are presented in Figures9 and 10. Due to the irregular spatial point configuration, the extents of the residual plots are not the same in all years. But the over- and underpredicted values have carefully been scaled while plotted to maintain the conformity in comparison.Figure 9
Residual plots from cross-validation of the four different spatial interpolation methods with their tendency of over- and under-prediction of PRCPTOT pattern. The relative sizes of the circles depict the relative degree of over- or underprediction in the same mm scale.Figure 10
Residual plots from cross-validation of the four different spatial interpolation methods with their tendency of over- and underprediction of TXx pattern. The relative sizes of the circles depict the relative degree of over- or underprediction in the same °C scale.The PRCPTOT and TXx residuals generated through IDW prove that the IDW method has typically underpredicted the indices values in spatial interpolation. Unlike IDW, TPS method has typically overpredicted the values of PRCPTOT and TXx in spatial interpolation. Among the stochastic methods, OK has mostly under predicted the PRCPTOT and TXx values in spatial interpolation, whereas UK has mostly overpredicted the indices values in spatial interpolation. The similarity in the pattern of over- and underprediction among the methods is also describable from their interpolation methodologies. On one hand, both the TPS and UK interpolate the long range variability of the climate over space and therefore overpredict the values than the observed ones in light of the higher values observed to create smoothed surfaces to the whole extent of study area. On the other hand, IDW interpolates the short range variability and OK, though interpolates the long range variability, takes anisotropy into account, which considers the short range variability of climate. Thus, taking the short range variability into account, both mostly underpredict the neighboring values of the low observations than the observed ones.Hence, the prediction performance of the stochastic methods is better than the one of the deterministic methods, since the size of the circle of difference between measured and predicted values is smaller in case of the stochastic methods’ residuals. Between the two deterministic methods, TPS shows better conformity to the observed values than the IDW. Importantly, all the methods’ conformity to the observed climate indices’ values increases over time with the increase in the number of the available spatial point observations, which can be proven by the gradually decreased size of the residual circles over time. The conformity of all the methods to the observed climate indices’ values gets also similar to each other with the increase in the available spatial points, as the comparative size of the circle represents.


### 4.6. The Trend of the Performance Measures

Two performance measures—index of agreement (d) and coefficient of variation of prediction (ρf)—have been calculated from the cross-validation (observed and predicted climate indices’ values) of the four different spatial interpolation methods. The d value approximately varies between 0 and 1 for both of the indices; the ρf value for PRCPTOT varies between 15 and 65 and for TXx between 0 and 40, respectively. The comparative linear trends of d and ρf measures of the different methods for interpolating PRCPTOT and TXx are represented in Figure 11. The figure reveals that there is a little difference in the prediction ability and the performance of the four different spatial interpolation methods, though they follow inherently different spatial interpolation methodology. The margin of the differences among the performance measures of the different methods is minimal in all the years, which proves the fact of the different spatial interpolation methods being equally accurate in climate prediction. In 1980 all the four spatial interpolation methods provide the least accurate surfaces of PRCPTOT with d close to 0 and ρf close to 65, and in 1967 all the four spatial interpolation methods provide the least accurate surfaces of TXx with d close to 0 and ρf close to 40. The better performance of all the methods is available in the years after 2000 with higher d values and lower ρf  values.Trend of the index of agreement (d) and the coefficient of variation of prediction (ρf) during the study period for the four different spatial interpolation methods in interpolating (a) PRCPTOT and (b) TXx.
(a)
(b)The trend ofd for PRCPTOT interpolation show that the UK interpolated PRCPTOT values shows the better agreement with the observed PRCPTOT values from the beginning until the beginning of the 1980s. It then loses in a hard competition with better performed OK and TPS; though OK performs best for the remaining years. For some of the years in between, UK again shows better performance but not as dominantly as before.The trend ofρf  for PRCPTOT indicates that IDW only performs better for two starting years followed by UK. OK performs equally well as UK at the very end of the study period by showing the less variation between the measured and predicted values. But in most of the years, UK method shows the least variation between the measured and predicted values, which is conformal to the d values that prove that UK method shows better agreement between the measured and predicted values.The trend ofd for TXx shows that UK and TPS perform equally well from 1948 until the beginning of the 1980s in terms of the agreement between the measured and predicted climate values. Then for the remaining years OK performs better. The trend of ρf for TXx indicates a very high competition of the methods in terms of their predictability; there is hardly any better method. Out of the very close performances, IDW performs better from the beginning until the mid 1960s and then UK performs better until the end in terms of the variability between the measured and predicted values. Like in case of PRCPTOT interpolation, UK method performs better in most of the years in terms of predictability.Two important observations have been made from the trend of the performance measures. The first observation is that the average trend values ofd of the cross-validation of the four different methods to interpolate PRCPTOT and TXx are 0.0052 and 0.0045, respectively. And the average trend values of ρf  of the cross-validation of the four different methods to interpolate PRCPTOT and TXx are −0.28 and −0.05, respectively. These results indicate that the agreement between the measured and predicted values in interpolation for all the four methods increases over the study period and the deviation between the measured and predicted values decreases. In a nutshell, the predictability of the spatial interpolation methods improves with time, and as discussed in Section 2, the number of spatial point observations increases with time.The conclusion is that the performance of the spatial interpolation improves with the increase in the available spatial points, which is a significant point to consider in spatial interpolation.The second as well the final observation includes a few important features of the applied four spatial interpolation methods, which are universally considerable in climate prediction. Table3 represents a comparative overview of the four different spatial interpolation methods in light of the four important performance parameters in the interpolation of PRCPTOT and TXx. Though the margin of differences among the performance measures of the four different spatial interpolation methods is minimal, the differences among the measures decrease further with time with the increase in the number of spatial points. Nevertheless, though modeled by different methodologies, the predictability of the four different spatial interpolation methods becomes similar if there are a certain number of spatial point observations available.Table 3
Different spatial interpolation methods.
Thin plate spline
Inverse distance weighting
Ordinary kriging
Universal kriging
Description of the climate variability
Smoothed spline exactly through the measured spatial points
Straight line exactly through the measured spatial points
Smoothed line that best fits through the measured semivariance considering anisotropy
Smoothed line that best fits through the measured semivariance considering spatial trend
Performance parameters
Performance in spite of the changing spatial point configuration
Performs comparatively worse than the stochastic methods and shows improved performance in cross-validation with the increased number of spatial point observations
Performs comparatively better than other methods when the number of spatial point observations is too little but in general shows improved performance in cross-validation with the increased number of spatial point observations
Performs comparatively better than the deterministic methods and shows improved performance in cross-validation with the increased number of spatial point observations
Performs best of all the methods applied and shows improved performance in cross-validation with increased number of spatial point observations
Measurement errors’ inclusion
Includes measurement errors partially in results because of the minimal smoothing
Includes measurement errors in results totally
Excludes measurement errors in result by smoothing
Excludes measurement errors in result by smoothing
Short versus long range variability
Describes mostly long range variability by removing the noise of short range variability
Describes only short range variability
Adjusts between the short and long range variability by incorporating anisotropy and semivariance
Describes mostly long range variability in light of the spatial trend


## 4.1. Justifying the Choice of Spatial Interpolation Methods

The calculated climate indices—PRCPTOT and TXx—show spatial trend and anisotropic behavior over the study area. The correlation of PRCPTOT with longitude (0.55) is higher than the correlation of PRCPTOT with latitude (−0.42), which indicates that the spatial trend is more dominant in the west-east direction than the north-south direction. On the other hand, the correlation of TXx with longitude (−0.52) is again higher than the correlation of TXx with latitude (0.34), which indicates that the spatial trend is more dominant in the west-east direction than the north-south direction. The anisotropic behavior of the indices is not very significant since the ratio of the major and minor axes varies from 0.7 to 1.0 (Table1), but it exists in the indices behavior. The results of the analysis of the temporal trend of the calculated climate change indices indicate that for PRCPTOT the range of the trend value is −18 to 36 for 60 years. Most of the stations show an increasing trend of PRCPTOT over time. Especially the stations in the south-east region of the country experience the highest increasing trend, where the highest values of PRCPTOT are also experienced. On the other hand, for TXx a range of −0.09 to 0.33 for the trend value has been obtained. Though most of the stations show an increasing trend of TXx, almost all the stations in the midregion, including capital Dhaka, show a decreasing trend.Table 1
Variogram models and parameters for OK interpolation of PRCPTOT and TXx.
Index
Temporal period
Variography
Model
Sill
Range
Nugget
Anisotropy parameters
Angle
Ratio of major and minor axis
PRCPTOT
1948–1975
Spherical
700000
400
45000
2.96
0.83
1976–1992
Spherical
630000
550
205000
47.28
0.97
1993–2007
Spherical
780000
550
70000
83.94
0.87
TXx
1948–1975
Spherical
15
430
0
2.96
0.71
1976–1992
Spherical
10
530
0
0.96
0.71
1993–2007
Spherical
5.8
410
1
176.72
0.73As  Isaaks and Srivastava [5], and Atkinson and Tate [66] discovered, that in case of uncertain spatial point configuration, deterministic methods should be the preferable option for spatial interpolation, two deterministic methods—thin plate spline (TPS) and inverse distance weighting (IDW)—have been applied for interpolation of the indices. These methods have mostly been applied in climate modeling by the geostatisticians. Based on the very low temporal trend of the indices, the pooled variogram approach is justified since a mean variogram of several temporal points should describe the spatial variability of the indices more accurately. Finally, after modeling the pooled variograms, ordinary kriging (OK) has been applied taking the anisotropic behavior of the indices into account, and universal kriging (UK) has been applied by taking the discussed spatial trend of the indices into account.


## 4.2. The Pooled Variograms

12 pooled variograms have been designed for modeling OK and UK for interpolating PRCPTOT and TXx. The dX parameter has been set to 0.55 to rescale the temporal distance of 1 year to the spatial distance of 550 km (the maximum distance obtained in the spatial point pairs) within the stacks. The total numbers of spatial point collected in each stack of the 1948–1975, 1976–1992, and 1993–2007 are 483, 494, and 503, respectively. The variogram parameters have been presented in Tables1 and 2. The anisotropy has been analyzed and modeled in the variograms of OK. The fundamental principle behind designing the UK variograms is that both PRCPTOT and TXx indices behavior over the study region are functions of longitude and latitude, which indicates that there is a spatial trend in the indices behavior. The UK variograms for interpolating PRCPTOT and TXx result in higher nugget values than the OK variograms, which represents that the spatial trend parameter models the long range spatial variability of the indices and therefore results in abrupt changes in the short range variability.Table 2
Variogram models and parameters for UK interpolation of PRCPTOT and TXx.
Index
Temporal period
Dependent parameters
Variography
Model
Sill
Range
Nugget
PRCPTOT
1948–1975
Longitude and Latitude
Spherical
300000
400
115000
1976–1992
Longitude and Latitude
Spherical
200000
550
260000
1993–2007
Longitude and Latitude
Spherical
270000
500
136000
TXx
1948–1975
Longitude and Latitude
Spherical
10
420
1
1976–1992
Longitude and Latitude
Spherical
2.1
480
1
1993–2007
Longitude and Latitude
Spherical
3.9
450
1.4OK variograms for TXx result in no nugget for the temporal periods of 1948–1975 and 1976–1992, whereas for 1993–2007 there is a nugget value of 1, which is high for this temporal period. The OK variogram for TXx of the temporal period of 1993–2007 also represents a lower sill value, while the OK variogram for PRCPTOT of the 1976–1992 period represents a higher sill value. These values show that the TXx values vary to a lower spatial extent during 1993–2007, while the PRCPTOT values vary to a higher spatial extent during 1976–1992 over the study region. The UK variograms for TXx of the temporal periods of 1976–1992 and 1993–2007 represent lower sill values, which indicate that TXx values did not vary to a great extent in the direction of spatial trend during the mentioned temporal periods. The ranges of all variograms are similar to the extent of the whole study region, which represents that the spatial variability of the indices is more accurately describable in the long range of the whole study area than in the short range of the neighboring stations.


## 4.3. The Interpolated Surfaces of the Climate Indices

The inherent methodologies of the four different spatial interpolation methods have been reflected in their produced surfaces of the climate indices. The search neighborhood [67, 68] has been set to the maximum number of available station observations in the corresponding years. The produced climate surfaces have been represented in Figures 5 and 6. To ease the comparison and visualization, surfaces from four different years have been selected as representatives of the three temporal periods for which the pooled variograms have been computed—1949, 1975, 1987, and 2006. And the following difference surfaces and residual plots will also be represented for these four different years to provide the readers with better understanding of the result. The PRCPTOT values range from 0 to 4500 mm, and TXx values range from 20° to 65°C.Figure 5
Representative and comparative PRCPTOT surfaces of different years produced by the four different spatial interpolation methods.Figure 6
Representative and comparative TXx surfaces of different years produced by the four different spatial interpolation methods.The effect of the property of the global information carried by the stationary mean in case of few observations [12] has been represented by the IDW surfaces for both PRCPTOT and TXx indices. Instead of representing the long range variability of the indices, IDW represents more variability close to the sample location in a circular manner and shows almost no variation in distance. This is identified by the “Bull’s eyes” that appeared in the surfaces in every station location, which are expected from the modeling features of IDW described in Section 3.2.1. Consequently, no clear spatial trend has been recognized from the surfaces, but rather some higher values have been recognized close to the areas of the representative stations.TPS surfaces, as a result of the inherent smoothing properties, are smoothed and present the long range variability of the phenomena. From the PRCPTOT surfaces, the spatial trend of the index is apparent, in all years the highest PRCPTOT has been experienced in the south-eastern part of the country. It is hard to predict any temporal trend in the surfaces for PRCPTOT. As described in the beginning of Section4, the spatial trend of TXx is opposite to the one of PRCPTOT, which has also been represented by the spline surfaces. The surfaces represent a spatial pattern of higher TXx in the north-western part of the country, which is dominant in almost all the years.OK method has been applied using the variography described in Section4.2. The surfaces describe the indices behavior significantly—both in terms of scale-dependent trend and anisotropy. Again it is difficult to predict any temporal trend from the surfaces. Like PRCPTOT surfaces, TXx surfaces produced by ordinary kriging do not describe any temporal trend in particular, but the spatial trend to the north-western part of the country is clearly represented by them.The tilted surfaces of UK have been created using the variograms designed in Section4.2. The basic principle, which has been applied to create these surfaces, is the fact that there is a spatial trend in indices behavior, which is also clearly visible in the created surfaces. The spatial trend is more dominant in the east-west direction than the north-south as expected from the correlation with longitude and latitude.From the climate surfaces presented in Figures5 and 6, a clear difference is visible in the surfaces in terms of the spatial pattern of climate variability. The difference occurs because of the different principles followed by different spatial interpolation methods apart from the fact that they interpolate the same observations from the stations.


## 4.4. The Surface-Differences of the Climate Indices

To analyze the degree and the margin of differences between the interpolated climate surfaces of PRCPTOT and TXx generated through the four different methods, the surface-difference has been computed. The computed difference surfaces are presented in Figures7 and 8. The surface-differences for both PRCPTOT and TXx show significant differences in the spatial pattern of the climate variability interpolated by different spatial interpolation methods. The range of differences observed among the PRCPTOT surfaces is 1000 mm to −1000 mm and among the TXx surfaces is 8°C to −8°C. The variation among the PRCPTOT surfaces is more diverse than the variation among the TXx surfaces and has therefore been reflected in their difference surfaces. This is a result of the difference between the spatial pattern of PRCPTOT and TXx; precipitation is not as continuous in time and space as temperature. In most cases, precipitation field represents a bounded continuum in its spatial and temporal pattern, and therefore the surfaces induce more noise.Representative difference surfaces showing the difference between the PRCPTOT surfaces produced by the four different spatial interpolation methods.Representative difference surfaces showing the difference between the PRCPTOT surfaces produced by the four different spatial interpolation methods.The maximum difference has been detected between the surfaces created by the two deterministic methods (TPS and IDW) for both PRCPTOT and TXx. The reason of the maximum difference between the deterministic surfaces is clear from their spatial modeling concept described in Section3.2.1. The dispute between the smoothing and nonsmoothing properties and therefore modeling short range and long range variability cause significant differences in the surfaces. The minimum difference has been revealed between the surfaces created by the two stochastic methods (OK and UK) for both PRCPTOT and TXx. The stochastic surfaces have minimum differences because the interpolators are identical and characterize the spatial variability through the variogram in similar ways. Thus, the stochastic methods are similar in their interpolation principle but model the variables based on the variable behavior (anisotropy and spatial trend) over space.The IDW method is unique among the methods as it creates the bull’s eyes in the surfaces, and therefore in the computed difference-surfaces between the IDW and others, significant differences have been observed in the areas within the close buffer of the station locations. The most variation in the differences has also been occurred in the areas nearby station locations. Therefore, the tendency of the IDW method being restricted to interpolate short range variability and including measurement errors is apparent from this analysis.The surfaces are more different in the areas with extreme indices values like TPS-IDW and TPS-OK surfaces of PRCPTOT of 1987. In the northern and eastern parts of the country, the differences in these surfaces exceed 1000 mm (Figure7). This indicates that different spatial interpolation methods predict the climate variables more differently when the climate behavior deviates than the trend. Another important fact is that the differences among the surfaces have been decreasing over years with the increasing number of spatial points and decreasing spatial lag, which indicates that the performance of the different spatial interpolation techniques becomes similar if there are more point observation available. This is an important property to be considered in decision making on the part of spatial interpolation, because the decision is more complicated where there are few point observationss available. If enough point observations are available, the interpolated surfaces by different spatial interpolation methods are mostly similar.


## 4.5. The Residual Plots from Cross-Validation

The residual plots of prediction have been created from the cross-validation of the four different spatial interpolation methods for creating the surfaces of PRCPTOT and TXx. The residual plots are presented in Figures9 and 10. Due to the irregular spatial point configuration, the extents of the residual plots are not the same in all years. But the over- and underpredicted values have carefully been scaled while plotted to maintain the conformity in comparison.Figure 9
Residual plots from cross-validation of the four different spatial interpolation methods with their tendency of over- and under-prediction of PRCPTOT pattern. The relative sizes of the circles depict the relative degree of over- or underprediction in the same mm scale.Figure 10
Residual plots from cross-validation of the four different spatial interpolation methods with their tendency of over- and underprediction of TXx pattern. The relative sizes of the circles depict the relative degree of over- or underprediction in the same °C scale.The PRCPTOT and TXx residuals generated through IDW prove that the IDW method has typically underpredicted the indices values in spatial interpolation. Unlike IDW, TPS method has typically overpredicted the values of PRCPTOT and TXx in spatial interpolation. Among the stochastic methods, OK has mostly under predicted the PRCPTOT and TXx values in spatial interpolation, whereas UK has mostly overpredicted the indices values in spatial interpolation. The similarity in the pattern of over- and underprediction among the methods is also describable from their interpolation methodologies. On one hand, both the TPS and UK interpolate the long range variability of the climate over space and therefore overpredict the values than the observed ones in light of the higher values observed to create smoothed surfaces to the whole extent of study area. On the other hand, IDW interpolates the short range variability and OK, though interpolates the long range variability, takes anisotropy into account, which considers the short range variability of climate. Thus, taking the short range variability into account, both mostly underpredict the neighboring values of the low observations than the observed ones.Hence, the prediction performance of the stochastic methods is better than the one of the deterministic methods, since the size of the circle of difference between measured and predicted values is smaller in case of the stochastic methods’ residuals. Between the two deterministic methods, TPS shows better conformity to the observed values than the IDW. Importantly, all the methods’ conformity to the observed climate indices’ values increases over time with the increase in the number of the available spatial point observations, which can be proven by the gradually decreased size of the residual circles over time. The conformity of all the methods to the observed climate indices’ values gets also similar to each other with the increase in the available spatial points, as the comparative size of the circle represents.


## 4.6. The Trend of the Performance Measures

Two performance measures—index of agreement (d) and coefficient of variation of prediction (ρf)—have been calculated from the cross-validation (observed and predicted climate indices’ values) of the four different spatial interpolation methods. The d value approximately varies between 0 and 1 for both of the indices; the ρf value for PRCPTOT varies between 15 and 65 and for TXx between 0 and 40, respectively. The comparative linear trends of d and ρf measures of the different methods for interpolating PRCPTOT and TXx are represented in Figure 11. The figure reveals that there is a little difference in the prediction ability and the performance of the four different spatial interpolation methods, though they follow inherently different spatial interpolation methodology. The margin of the differences among the performance measures of the different methods is minimal in all the years, which proves the fact of the different spatial interpolation methods being equally accurate in climate prediction. In 1980 all the four spatial interpolation methods provide the least accurate surfaces of PRCPTOT with d close to 0 and ρf close to 65, and in 1967 all the four spatial interpolation methods provide the least accurate surfaces of TXx with d close to 0 and ρf close to 40. The better performance of all the methods is available in the years after 2000 with higher d values and lower ρf  values.Trend of the index of agreement (d) and the coefficient of variation of prediction (ρf) during the study period for the four different spatial interpolation methods in interpolating (a) PRCPTOT and (b) TXx.
(a)
(b)The trend ofd for PRCPTOT interpolation show that the UK interpolated PRCPTOT values shows the better agreement with the observed PRCPTOT values from the beginning until the beginning of the 1980s. It then loses in a hard competition with better performed OK and TPS; though OK performs best for the remaining years. For some of the years in between, UK again shows better performance but not as dominantly as before.The trend ofρf  for PRCPTOT indicates that IDW only performs better for two starting years followed by UK. OK performs equally well as UK at the very end of the study period by showing the less variation between the measured and predicted values. But in most of the years, UK method shows the least variation between the measured and predicted values, which is conformal to the d values that prove that UK method shows better agreement between the measured and predicted values.The trend ofd for TXx shows that UK and TPS perform equally well from 1948 until the beginning of the 1980s in terms of the agreement between the measured and predicted climate values. Then for the remaining years OK performs better. The trend of ρf for TXx indicates a very high competition of the methods in terms of their predictability; there is hardly any better method. Out of the very close performances, IDW performs better from the beginning until the mid 1960s and then UK performs better until the end in terms of the variability between the measured and predicted values. Like in case of PRCPTOT interpolation, UK method performs better in most of the years in terms of predictability.Two important observations have been made from the trend of the performance measures. The first observation is that the average trend values ofd of the cross-validation of the four different methods to interpolate PRCPTOT and TXx are 0.0052 and 0.0045, respectively. And the average trend values of ρf  of the cross-validation of the four different methods to interpolate PRCPTOT and TXx are −0.28 and −0.05, respectively. These results indicate that the agreement between the measured and predicted values in interpolation for all the four methods increases over the study period and the deviation between the measured and predicted values decreases. In a nutshell, the predictability of the spatial interpolation methods improves with time, and as discussed in Section 2, the number of spatial point observations increases with time.The conclusion is that the performance of the spatial interpolation improves with the increase in the available spatial points, which is a significant point to consider in spatial interpolation.The second as well the final observation includes a few important features of the applied four spatial interpolation methods, which are universally considerable in climate prediction. Table3 represents a comparative overview of the four different spatial interpolation methods in light of the four important performance parameters in the interpolation of PRCPTOT and TXx. Though the margin of differences among the performance measures of the four different spatial interpolation methods is minimal, the differences among the measures decrease further with time with the increase in the number of spatial points. Nevertheless, though modeled by different methodologies, the predictability of the four different spatial interpolation methods becomes similar if there are a certain number of spatial point observations available.Table 3
Different spatial interpolation methods.
Thin plate spline
Inverse distance weighting
Ordinary kriging
Universal kriging
Description of the climate variability
Smoothed spline exactly through the measured spatial points
Straight line exactly through the measured spatial points
Smoothed line that best fits through the measured semivariance considering anisotropy
Smoothed line that best fits through the measured semivariance considering spatial trend
Performance parameters
Performance in spite of the changing spatial point configuration
Performs comparatively worse than the stochastic methods and shows improved performance in cross-validation with the increased number of spatial point observations
Performs comparatively better than other methods when the number of spatial point observations is too little but in general shows improved performance in cross-validation with the increased number of spatial point observations
Performs comparatively better than the deterministic methods and shows improved performance in cross-validation with the increased number of spatial point observations
Performs best of all the methods applied and shows improved performance in cross-validation with increased number of spatial point observations
Measurement errors’ inclusion
Includes measurement errors partially in results because of the minimal smoothing
Includes measurement errors in results totally
Excludes measurement errors in result by smoothing
Excludes measurement errors in result by smoothing
Short versus long range variability
Describes mostly long range variability by removing the noise of short range variability
Describes only short range variability
Adjusts between the short and long range variability by incorporating anisotropy and semivariance
Describes mostly long range variability in light of the spatial trend


## 5. Conclusion

This paper analyzes the degree and margin of differences among the surfaces of the PRCPTOT and the TXx of Bangladesh produced by the TPS, IDW, OK, and UK methods of spatial interpolation. The indices have been calculated from the daily precipitation and temperature measurements of the time series between 1948 and 2007 distributed among the 34 meteorological stations. The results indicate that the performances of the four different spatial interpolation methods become less uncertain and more similar with the increase in the available spatial data points. Though it is hard to forecast the exact number of required spatial points for the spatial interpolation methods to provide exactly same surfaces, according to Kelley [69], to interpolate with 95% confidence where the desired full confidence interval width is 10% (ω=0.10) and desired degree of assurance is 99% (γ=0.99), the required number of spatial points for the PRCPTOT is 16,233 and for the TXx is 199. From this point of view, the maximum available number of spatial points (34) is too low to obtain acceptable accuracy. The required number of spatial points to interpolate with acceptable accuracy and to compute variogram with less uncertainty is subject to further studies. Additionally, a brief discussion of the temporal trends in the climatic indices and their spatial distribution has been included in the paper within the scope. It would be interesting for future studies to perform temporal trend analysis on the spatially interpolated data and examine the effects of different interpolation methods on the results.The research has some inherent limitations which cannot be overlooked but have been treated with utmost importance. In designing pooled variograms for the stochastic methods, the spherical model has been chosen for all variograms of all the climate indices. This is mainly related to the semivariances of the indices over space, which shows spherical spatial dependence. Variogram models were fitted iteratively to the experimental variograms as recommended by Goovaerts [12], thus the variogram parameters were estimated subjectively by giving more relevance to the first lags. Furthermore, other suitable variogram models have also been applied, but the spherical model also proved to be appropriate in the cross-validation. The weighted regression method and the cross-validations have been considered in the determination of the variogram parameters. This method is indeed similar to the weighted least squares method (WLS). However there are scopes to try more probabilistic procedures such as maximum likelihood or restricted maximum likelihood.Importantly, the nugget effect under OK is small, while under UK it is very high. This clearly indicates that the spatial trend is affecting the short-range variability of the indices. The reason for this effect can be seen in the spatial trend for the indices occurring in long range (across the whole country). Additionally, there is not enough point observations to model the short-range variability and to describe the gradual change in the indices’ values in the trend direction. Given that the smallest accessible lag of the empirical variogram is large, the nugget effect is likely to be overestimated. Additionally, the nugget value and the semivariogram behavior at the origin cannot be cross-validated because variogram model values for lags smaller than the shortest sampling interval do not intervene in interpolation algorithms.At last but not least, selection of the appropriate spatial interpolation method to interpolate a particular climate phenomenon is a hard decision to make. From a very concise view, to interpolate the spatial pattern of the PRCPTOT and TXx in Bangladesh, UK method is the best choice to make. This decision is highly influenced by the spatial trend in the indices behavior and the inadequately available spatial point observations. And this decision is highly specific for these two climate indices and for interpolating the surfaces of Bangladesh despite the universal features of the applied four spatial interpolation methods.The research has broader perspective in meteorology and climatology. Comparing the spatial interpolation methods and eventually determining the locally appropriate method will contribute to preparing appropriate climate dataset and understanding their spatial variability in local scale. Moreover, these high-resolution local continuous surfaces of the climate indices can be input in several climate models to forecast climate change phenomenon and related consequences on the developing countries like Bangladesh, where climate studies are suffering from unavailability of high-resolution climate data recording and monitoring. To conclude, choosing appropriate methodology of spatial interpolation as ideal meteorological (climate) method will result in appropriate climate forecast and climate resilience activities.



---

*Source: 353408-2012-11-29.xml*