Incorporation of Parameter Uncertainty into Spatial Interpolation Using Bayesian Trans-Gaussian Kriging

Joon Jin SONG,Soohyun KWON,and GyuWon LEE?

1Department of Statistical Science,Baylor University,USA

2Department of Astronomy and Atmospheric Sciences,Research and Training Team for Future Creative Astrophysicists and Cosmologist,Kyungpook National University,Republic of Korea

Incorporation of Parameter Uncertainty into Spatial Interpolation Using Bayesian Trans-Gaussian Kriging

Joon Jin SONG1,Soohyun KWON2,and GyuWon LEE?2

1Department of Statistical Science,Baylor University,USA

2Department of Astronomy and Atmospheric Sciences,Research and Training Team for Future Creative Astrophysicists and Cosmologist,Kyungpook National University,Republic of Korea

Quantitative precipitation estimation(QPE)plays an important role in meteorological and hydrological applications. Ground-based telemetered rain gauges are widely used to collect precipitation measurements.Spatial interpolation methods are commonly employed to estimate precipitation felds covering non-observed locations.Kriging is a simple and popular geostatistical interpolation method,but ithastwo known problems:uncertainty underestimation andviolation of assumptions. This paper tackles these problems and seeks an optimal spatial interpolation for QPE in order to enhance spatial interpolation through appropriately assessing prediction uncertainty and fulflling the required assumptions.To this end,several methods are tested:transformation,detrending,multiple spatial correlation functions,and Bayesian kriging.In particular,we focus on a short-term and time-specifc rather than a long-term and event-specifc analysis.This paper analyzes a stratiform rain event with an embedded convection linked to the passing monsoon front on the 23 August 2012.Data from a total of 100 automatic weather stations are used,and the rainfall intensities are calculated from the difference of 15 minute accumulated rainfall observed every 1 minute.The one-hour average rainfall intensity is then calculated to minimize the measurement random error.Cross-validation is carried out for evaluating the interpolation methods at regional and local levels.As a result, transformation is found to play an important role in improving spatial interpolation and uncertainty assessment,and Bayesian methods generally outperform traditional ones in terms of the criteria.

precipitation,kriging,transformation,Bayesian kriging,detrend,Korea

1. Introduction

Quantitative precipitation estimation(QPE)plays an important role in meteorological and hydrological applications. Rain gaugesare widely used to collect precipitationmeasurements dueto certainadvantages.Forinstance,rain gaugesdirectlymeasurerainfallonthe ground,thusprovidingaccurate ground-level precipitation observations with limited error.

Precipitation varies signifcantly in time and space. Hence,spatial interpolation methods are commonly employedto estimate precipitationin locationslackingmeasurement equipment.Examples of such interpolationmethodsinclude inverse distance weighting(Franke,1982),local polynomial(Yilmaz,2007),and radial basis function(Carlson and Foley,1991).Spatial interpolation methods are typically classifed into two categories:deterministic and stochastic. The deterministic methods provide no assessment of possible errors,while the stochastic methods offer probabilisticestimates(Dirks et al.,1998;Nalder and Wein,1998;Buytaert et al.,2006;Basistha et al.,2008;Ly et al.,2011).In this paper,we focus on geostatistical stochastic methods,in particular ordinary kriging(OK)and its variants are considered as a spatial interpolation tool.

Kriging is a simple and popular geostatistical interpolation method.There are several variants,such as simple kriging,ordinary kriging,universal kriging,and indicator kriging(Cressie,1993;Schabenberger and Gotway,2004).This paper tackles two problems that are often neglected in kriging analysis—uncertainty underestimation and violation of assumptions—in a single framework.

Although kriging is widely used for spatial interpolation, the spatial structure of the underlying process is presumed known,leading to a plug-in or two-stage procedure.Hence, kriging is often performed after estimating the parameters of the spatial structure.In this framework,the uncertainty in spatial interpolation is often underestimated because the procedureignoresuncertaintyintheparameterestimationrelated to spatial structure,as if the parameters were the true values.This is an optimistic assessment of predictive accuracy.To overcome this problem,there are several ways to account for the uncertainty,including bootstrapping(Wang and Wall, 2003)and Bayesian statistics(Diggle et al.,1998;Handcock and Stein,1993).The Bayesian approach allows explicit accounting of uncertainties in the model parameters by treating them as random quantities,rather than as unknown constants as in the classical approach.Within this framework, parameter estimation and prediction can be conducted simultaneously without the separation in the two-stage procedure above.This framework accounts for the uncertainty ignored by model parameterestimation,leading to more realistic spatial prediction,with better uncertainty assessment.In addition,any a priori knowledge about the unknown quantities can be incorporated into the inference.

The second problem is that spatial data sets in practice often violate the assumptions required for kriging,and it is neglected in spatial analysis.Kriging may be used to fnd the best linear unbiased predictor(BLUP),as it fully complies with the validity of all the required assumptions, such as normality,stationarity,and homoscedasticity(Isaaks and Srivastava,1989).Furthermore,the effect of violating the assumption on spatial prediction would be substantial.This paper will investigate suitable remedies for the violations.

The main objective of this study is to obtain an optimal spatial interpolation for QPE to assess prediction uncertainty appropriately and fulfll the required assumptions. To this end,we compare several kriging variants:(1)transformation;(2)detrending;(3)spatial correlation function; and(4)Bayesian kriging.Very few studies have explored these issues simultaneously in a single framework.We focus on a short-term and time-specifc analysis rather than longterm and event-specifc ones.The data and methodology are shown in section 2 and the analysis results are discussed in section 3.

2. Data and method

2.1.Data and pre-processing

The rain gauge data are collected by tipping bucket rain gaugesin an AutomaticWeather System(AWS)with 0.5 mm resolution.Thus,0.5 mm h-1is used as the cut-off threshold to reject dry areas experiencing no rain from the analysis.A total of 100 AWS stations are used over the area of(34.33°–37.03°N,126.84°–128.25°E).The rainfall intensities are calculated from the difference between the AWS-observed 15 minute accumulated rainfall amounts every minute.The one hour average rainfall intensity is calculated to minimize the measurement random error.Figure 1 displays the locations of the rain gauges in the study area.The event to be analyzed is a stratiform precipitation event with embedded convection related to the passage of the monsoon front that occurred on 23 August 2012.The convective line developed over the southeast of the domain associated with the monsoon front(Fig.2).Two lines intensifed from 0600–0900 LST,and the systems became weaker and spread widely.The rain band moved southeast to northwest.The total rainfall amounts were recorded up to a maximum of 173 mm in the analysis region.

The precipitation data are summarized in Table 1.Since only wet stations are used in the analysis,the number of stations considered varies from 33 to 76 over the time steps of thestudyperiod.Therewas intenseprecipitationduringsome time steps,e.g.time steps 4(42 mm h-1),6(59.5 mm h-1), and 13(84.5 mm h-1).Data variation,measured by the standarddeviation(SD),wasnotconstantovertime,rangingfrom 1.402 mm h-1to 10.262 mm h-1.This fnding motivates hourly-specifc spatial analysis rather than aggregated analysis over all time steps in a single event.

2.2.Transformation

Some inherent characteristics of precipitation lead to violation of the normality and constant variance(homoscedasticity)conditionsnecessary for OK to be BLUP.For instance, precipitation data are commonly skewed,due to their intermittency;data transformations have been commonly adopted to remedy this.Although square root and logarithmic transformations have been used(e.g.,Schuurmans et al.,2007; Verworn and Haberlandt,2011),they do not always suffciently account for the departure from the assumptions.An alternative option is the family of power transformations, including the Box–Cox transformation,which is the mostwidely used(Box and Cox,1964).This is given by:

Table 1.Summary of the 24 hour precipitation data.Q1 is quartile 1,Q3 is quartile 3,Max is maximum,and SD is standard deviation.

whereZis the transformed data,Yis the original data larger than the threshold,andλis the transformation parameter.This transformation is typically applied to positive data,which exclude dry areas in this study.The transformation is data-driven because the transformation parameter is determined according to profle log-likelihoods over the value within some ranges.In this paper,we consider a timespecifc rather than an event-specifc transformation,because the same transformation over time in an event seems to be unreasonable due to variation of the precipitation process.

2.3.Variogram model

Kriging requires spatial pattern information.An empirical variogram is frst computed and ftted to a theoretical variogram in order to estimate spatial parameters,such as sill, range,and nugget,via

whereN(h)is the set of pairs of observationsZ(si)andZ(sj) such that distance between two locationsd(si,sj)is equal to spatiallaghand|N(h)|is thenumberofthepairs.Variogrambased parameter estimation is generally ineffcient because it is based on the smoothed variogram,which is not the original but a summary of the data.Likelihood-based methods are a general means to make use of the data generating process,but this approach requires a distributional assumption, such as normality.In this study,we adopt a maximum likelihood estimation that is widely implemented in statistical inference.Maximum likelihood estimators require a spatial distribution to construct the likelihood function.In this case,letZ=(Z(s1),...,Z(sn))Tdenote the vector of transformed observations,with a multivariate normal distribution withmeanμ nnn,andcovariancematrixΣ(θ),wherennnis an×1 vector of ones,andθis the vector of spatial parameters such as partial sill andrange.Theresultinglog-likelihoodfunction is given by

The maximumlikelihoodestimators can be obtainedby maximizingthe log-likelihoodfunctionwithrespecttotheparametersμandθ.The resulting estimator has meaningful statistical properties under some mild regularity conditions(Cox and Hinkley,1974).

Exponential and spherical functions have often been used as spatial correlation functions in variogram modeling and kriging analysis(Chil`es and Delfner,1999).In addition to thesefunctions,thecircularcorrelationfunctionisconsidered in this study,a frst in the study of rain.Table 2 presents the forms of the functions.Surprisingly,as shown in the results below,the new correlation function is the most frequentlyselected as the optimal spatial correlation function over time. For every time step,three correlation models are ftted,and one of them is selected according to a model selection criterion,known as the Bayesian Information Criterion(BIC). Therefore,a different correlation model is ftted to each time step.

The mean of the function often spatially varies over the region of interest,while one of the required assumptions for BLUP is a constant mean.In this case,detrending can deal with the non-constantmean problem.A trend surface is commonly modeled using spatial coordinates or available covariate information,and the residuals between the observations and the ftted trends delineate the spatial structure.As a result,this detrending,or removal of trends,can reduce the variability of the predictive distribution.We examine the effect of detrending on spatial prediction for spatial coordinates,longitude and latitude.

2.4.Interpolation method

OK is a linear interpolation method that is unbiased and minimizes the variance of the observations.The weightings of the linear interpolator are found by solving a system of equations with some constraints in orderto achieve the sound properties.Several variants of OK are considered in this study.Trans-Gaussian OK(TOK)is a variant of kriging with a transformed Gaussian random feld when the transformation is known(Cressie,1993).Application of the Box–Cox transformation is assumed to transform non-normally distributed data into a normal distribution.Bayesian OK(BOK) and Bayesian trans-Gaussian OK(BTOK)perform OK and TOK from a Bayesian perspective.The kriging variants with detrending based on spatial coordinates are also considered.

Bayesian analysis requires estimation of prior distributions for unknown parameters(θ).Combining the likelihood functionL(θ|Z)with a prior distributionP(θ)leads to an expression for the posterior distributionp(θ|Z)via normalized Bayes’Theorem:

The posterior distribution provides a probability statement about the parameters and allows for uncertainty in allparameters.Similarly,the Bayesian predictive distributionp(Z(s0)|Z)for an arbitrary and unobserved locations0can be obtained as

Table 2.The three spatial correlation functions;?is a range parameter that varies for each function,andhis spatial lag.

In this study,we choose non-informative priors due to a lack of a priori information about the parameters.

2.5.Validation

Cross-validationis carriedout to evaluate the infuenceof data transformation,detrending,spatial autocorrelation,and different interpolation methods(i.e.classical and Bayesian) on interpolation performance.LetZ(i,k)and?Z(i,k)denote the observed and predicted values from the leave-one-out cross-validation at theith monitoring station in thekth time step,respectively.The mean absolute error(MAE)is then employed to evaluate the quality of the interpolations,

wherenkis the number of wet stations in thekth time step.This criterion measures the unbiasedness of the crossvalidation prediction.A correlation coeffcient is used to evaluate the agreement between the observed and predicted values,

To compare the prediction uncertainty,two measures are employed,the length of the prediction interval and the coverage probability.The prediction intervalLkis formed from the predicted value and its prediction error,and a wider interval represents greater uncertainty.The former is the average length of the cross-validation prediction interval over all wet stations at a given time step,

whereU(i,k)andL(i,k)are the upper and lower bounds of the prediction intervals at stationiat a given time stepk,respectively.Thecoverageprobabilityis computedbycounting how many times the observed values fall in the prediction intervals, whereIA(x)is the indicator function,which is 1 forx∈A, and 0 otherwise.It is expected that the coverage probability is close to the nominal value(e.g,95%).

The methods considered in this study are evaluated on two spatial scales,regional and local(Xie et al.,2011).The regional-scale evaluation is performed over the entire study area for each time step,while the local-scale evaluation compares the performance of the different spatial interpolation methods at each station over all time steps.

3. Results

3.1.Transformation

Data transformation is employed to examine the assumptions required for the optimal spatial interpolation methods. Most statistical interpolation methods assume data are normally distributed.Figure 3 describes the distributions of raw and transformed precipitation datasets over 24 time steps using boxplots.The raw precipitation data include outliers in all time steps,but not the transformed data.This highlights the advantage of data transformation because outliers easily lead to non-normal data distributions.Robust distributions are more appropriate for data with outliers,such as thetdistribution.

Other statistical measures and testing are further used to study the impact of the power transformation on distribution normality(Table 3).Skewness is a measure of the extent of symmetry of a distribution,and kurtosis is a descriptor of the shape of a distribution,measuring the peakedness or fatness of a distribution.Positive kurtosis indicates a peaked distribution,while negative kurtosis corresponds to a fat distribution.The normal distribution has zero skewness and kurtosis by defnition.Table 3 reports that positive skewness and large kurtosis are found in raw data over all time steps,indicating asymmetric and positively skewed distributions.This is expected for rainfall data,which are typically skewed toward heavy rain.After transformation,such data have signifcantly reduced kurtosis.For testing normality quantitatively,the Kolmogorov–Smirnov test is performed,and the results are presented in Table 3,before and after transformation.For raw data,none of the stations’data have a normal distribution,although two stations fulfll that assumption at the 0.05 signifcance level,after applying the transformation.Although the remaining stations havep-values less than 0.05,the symmetry is much enhanced in terms of both skewness andkurtosis.Figure4 presentsnormalquantile–quantile plots of the raw and Box–Cox transformed datasets over two timesteps withtransformationparameterestimates.Itis clear that the transformation improves the normality.A signifcant deviation from normality is found in the tails of the distribution of the raw data,though this is greatly reduced in the transformed data.

3.2.Variogram ftting

Table 3.Distribution measures and Kolmogorov–Smirnov normality test metrics of the raw and transformed data.The symbol*indicates statistical non-signifcance at the 0.05 level.

To optimize ftting of the variogram,we consider two potential factors that can improve spatial prediction.The frst is the spatial correlation function,which models spatial structure,and thereforeleads to poorspatial predictionif specifed incorrectly.This effect is more pronounced for classical estimation,whichassumes thattheselectedfunctionis trueinthe predictionstage.Surprisingly,the spatial correlationfunction variessignifcantlyintime.We fndthattheunder-usedcircular correlationfunctionis most suitable(Table4).Detrending does not signifcantly affect the selection,whereas transformation does.

The signifcance of the trend(i.e.the variation of the mean in space)is investigated at each time step,and Table 5 summarizes the result.Over 67%of raw and 83%of transformed time steps have a signifcant trend effect.This motivates modeling of the trend surface in order to characterize the spatial structure for spatial interpolation.In particular, latitude is an important factor for the modeling trend in this analysis.Figure 1 illustrates the signifcant latitudinal variation of the total rain amount.

Table 4.Frequency table for the correlation functions selected for every time step.“Raw”is raw data;“Trans”is transformed data;“Trend”is constant mean;“Detrend”is spatially varying mean.

Table 5.Frequency table for the signifcant trend effects of raw and transformed data for every time step.“None”indicates no variation of mean,while“Lon”or“Lat”indicates the variation mean in the longitudinal or latitudinal direction.

3.3.Evaluation

Thecomparisonofdifferentspatial interpolationmethods is shown in terms of MAE in Fig.5 as a box-plot.The BTOK approach results in the smallest median,whereas the TOKD approach has the largest.In general,the Bayesian methods outperform the traditional methods in terms of MAE.The distributions of the MAE values of BTOK and TOK are narrower than those of the other methods.There are no outlying values(circles in Fig.5)of MAE in BTOK and BTOKD.

The BTOK has the largest correlation coeffcient and the smallest variation.Detrending does not improve the correlation.Particularly,BOKD and BTOKD have low correlation coeffcient values.Overall,the Bayesian methods perform better than the classical methods,in terms of both criteria. In contrast to detrending,transformation enhances the unbiasedness of the spatial interpolation.

Local-scale evaluationat each rain gaugeis shown in Fig. 6.Similar to the regional scale,transformation-based methods without detrending have,in general,a smaller MAE on the local scale,and Bayesian methods overall outperform the classical methods.TOK and BTOK perform best in terms of the local MAE.As for correlation coeffcients,TOK and BTOK have larger values than the other methods.Detrending has a negative effect on both MAE and the correlation coeffcient,whereas transformation and the Bayesian methods provide some improvement.

Two aspects of uncertainty estimation are compared. First,the average coverage probabilities for each time step are comparedfor the classical and Bayesianmethodsforeach kriging variant.Second,prediction intervals are constructed with 95%nominal probability.Hence,a good prediction interval has a coverage probability close to the nominal probability.

Overall,the Bayesian methods generally have reasonable ranges of coverage probability,while the classical methods have some extreme probabilities(Fig.7).The coverageprobabilities for the Bayesian methods(y-axis)are reasonably spread with some consistency,while those of the classical methods range too widely with extreme values.

The average lengths of the prediction intervals based on the spatial interpolation methods are shown in Fig.8.As addressed earlier,the width of the intervals represents uncertainty in the prediction.As expected,Bayesian methods have generally longer widths than the classical ones,and this indicates that the Bayesian approach evaluates a more realistic uncertaintyin spatial prediction.There is no systematic effect of transformation and detrending on the length.The scatterplot of prediction intervals(Fig.9)ensures that Bayesian methods have longer intervals.There are some signifcantlengthdifferencesbetweenBTOKandTOK.Overall, the Bayesian approaches have realistic and larger uncertainty estimation.Figure 10 presents the spatial precipitation estimated by BTOK for four time steps.

4. Summary and discussion

In this paper,we investigate two problems commonly foundand oftenneglectedin the spatial interpolationof QPE: uncertainty underestimation and violation of assumptions, along with their effects on spatial prediction and uncertainty estimation.The methods addressed are considered to resolve the problems,and implemented in a single framework.The proposed method is illustrated with a rain gauge dataset consisting of 100 AWS stations in South Korea.A stratiform precipitation event is analyzed with one hour average rainfall intensity in order to minimize measurement error.The proposed kriging variants are applied to the dataset and compared with several criteria at regional and local levels.

Overall,the methods improve spatial interpolation.For instance,we fnd that transformation plays an important role in improving spatial interpolation and is not constant over each time and event.Hence,it is challenging to fnd an appropriate transformation for each time step and event.The same transformation is commonly applied to all time steps in the event.Time-specifc transformation is recommended, even though this may demand slightly more computational resources.In summary,a desirablemethodfor a reliable QPE must be fexible and sophisticated enough to account for dynamicprecipitationprocesses,as opposedtoastaticapproach that assumesthe sametransformationorspatial structureover time in a given event.

As used in this study,parametric spatial functions typically model spatial structure.However,they are often too restrictive to explain the structure adequately.A fexible spatial structure function is necessary to characterize the spatial structure of precipitation,e.g.by using a nonparametric approach(Yao,1998).As a parametric approach,the Mat′ern class of covariance functions is a reasonable alternative(Handcock and Stein,1993).

Detrending is used as a means of correcting the nonstationarity of the underlying process.In this study,the trend is modeled using only the spatial coordinates,longitude and latitude.Although detrending appears to infuence the selection of the spatial correlation function,this method barely impactsontheevaluationofspatialinterpolation.However,itis possible that detrending can be a powerful tool if the trend is modeled adequately with the other potential variables.For the sake of convenience,only spatial coordinate variables are considered in this study.Elevation and weather radar data could help to characterize the precipitation feld and supplement spatial variation.

A Bayesian approach is introduced to assess prediction uncertainty related to the uncertaintyof parameter estimation in variogram analysis.There are some obstacles when conducting Bayesian prediction practically,to the frst of which is eliciting prior information.When there is little a priori information about the parameters,a non-informativeor diffuse prior distribution is often chosen regularly.In this case,the subsequentresult is dominatedby likelihoodor data information when the sample size is suffciently large.Furthermore, although the sample size is large in spatial prediction,the effective sample size is much smaller than the original sample size due to the high correlation between the spatial data.It is clear that prior information enhances not only parameter estimation but also the spatial prediction.The second concern is that the Bayesian approach usually requires intensive computationfor evaluatingthe integralrelated to fndingposterior and predictive distribution.If some convenient prior data are used,a conjugate prior distribution,the evaluationbecomes simple because they can be evaluated analytically. Otherwise,more feasible evaluation strategies are necessary, such as Markov chain Monte Carlo(MCMC)methods.For large sample sizes,the problem becomes more serious due to high dimensionality leading to inversion of the high dimension matrix.

Temporal processes are often modeled to account for the temporal variation of the precipitation process.It is obvious that temporal processes can enhance the modeling and prediction of the underlying process if the level of temporal variation is signifcant and the corresponding temporal process can take account of the variation appropriately.There are a variety of spatial–temporal processes available for this purpose(Guttorpet al.,1994;Cressie andHuang,1999;Genton,2007;Ma,2008).However,it is necessary to be aware that these do not always improve the spatial prediction because complicated models can worsen it.The complexity of the processes also increases the computation time,which can be an important concern in practical operations.

In QPE,it is usual to observe outliers in the precipitation dataset.Classical variogram estimation is sensitive to outliers,which can propagate in spatial predictions.TOK-and BTOK-transformeddatasets have no outliers.As such,transformation is helpful to not only comply with the required assumptions,but also to deal with outliers.Transformation can also remedy heteroscedasticity,required in kriging for fulflling second-order stationarity(Erdin and Frei,2012).

Acknowledgements.This work was funded by the Korea Meteorological Administration Research and Development Program (Grant No.CATER 2013-2040).JJS’s research was supported by the Brain Pool program of the Korean Federation of Science and Technology Societies(KOFST)(Grant No.122S-1-3-0422).

REFERENCES

Basistha,A.,D.S.Arya,and N.K.Goel,2008:Spatialdistribution of rainfall in Indian Himalayas—A case study of Uttarakhand Region.Water Resour.Manag.,22,1325–1346.

Box,G.E.P.,andD.R.Cox,1964:Ananalysisof transformations.Journal of the Royal Statistical Society(B),26,211–252.

Buytaert,W.,R.Celleri,P.Willems,B.de Bi`evre,and G.Wyseure, 2006:Spatial and temporal rainfall variabilityin mountainous areas:A case study from the south Ecuadorian Andes.J.Hydrol.,329,413–421.

Carlson,R.E.,and T.A.Foley,1991:The parameter R2in multiquadric interpolation.Computers&Mathematics with Applications,21,29–42.

Chil`es,J.P.,and P.Delfner,1999:Geostatistics:Modeling Spatial Uncertainty.John Wiley and Sons,New York,731 pp.

Cox,D.R.,and D.V.Hinkley,1979:Theoretical Statistics.Chapman and Hall,London,528 pp.

Cressie,N.,1993:Statistics for Spatial Data.Wiley-Interscience, 928 pp.

Cressie,N.,and H.-C.Huang,1999:Classes of nonseparable, spatio-temporal stationary covariance function.Journal of the American Statistical Association,94,1330–1340.

Diggle,P.J.,J.A.Tawn,and R.A.Moyeed,1998:Model-based geostatistics(with discussion).Applied Statistics,47,299–350.

Dirks,K.N.,J.E.Hay,C.D.Stow,and D.Harris,1998:Highresolution of rainfall on Norfolk Island,Part II:Interpolation of rainfall data.J.Hydrol.,208,187–193.

Erdin,R.,and C.Frei,2012:Data transformation and uncertainty in geostatistical combination of radar and rain gauges.J.Hydrometeor.,13,1332–1346.

Franke,R.,1982:Scattered data interpolation:Test of some methods.Mathematics of Computations,33,181–200.

Genton,M.G.,2007:Separable approximations of space-time covariance matrices.Environmetrics,18,681–695.

Guttorp,P.,W.Meiring,and P.D.Sampson,1994:A space-time analysis of ground-level ozone data.Environmetrics,5,241–254.

Handcock,M.S.,and M.L.Stein,1993:A Bayesian analysis of kriging.Technometrics,35,403–410.

Isaaks,E.H.,and R.M.Srivastava,1989:Introduction to Applied Geostatistics.Oxford University Press,Oxford,561 pp.

Ly,S.,C.Charles,and A.Degr′e,2011:Geostatistical interpolation of daily rainfall at catchment scale:The use of several variogram models in the Ourthe and Ambleve catchments,Belgium.Hydrology and Earth System Sciences,15,2259–2274.

Ma,C.S.,2008:Recent developments on the construction of spatio-temporal covariance models.Stochastic Environmental Research and Risk Assessment,22,S39–S47.

Nalder,I.A.,and R.W.Wein,1998:Spatial interpolation of climatic Normals:Test of a new method in the Canadian boreal forest.Agricultural and Forest Meteorology,92,211–225.

Schabenberger,O.,and C.A.Gotway,2004:Statistical Methods for Spatial Data Analysis.Chapman and Hall,512 pp.

Schuurmans,J.M.,M.F.P.Bierkens,E.J.Pebesma,and R.Uijlenhoet,2007:Automatic prediction of high-resolution daily rainfall felds for multiple extents:The potential of operational radar.J.Hydrometeor.,8,1204–1224.

Verworn,A.,and U.Haberlandt,2011:Spatial interpolation of hourly rainfall-effect of additional information,variogram inference and storm properties.Hydrology and Earth System Sciences,15,569–584.

Wang,F.J.,and M.M.Wall,2003:Incorporating parameter uncertainty into prediction intervals for spatial data modeled via a parametric variogram.Journal of Agricultural,Biological, and Environmental Statistics,8,296–309.

Xie,H.,X.Zhang,B.Yu,and H.Sharif,2011:Performance evaluation of interpolation methods for incorporating rain gauge measurements intoNEXRADprecipitationdata:Acasestudy in the upper Guadalupe river basin.Hydrological Processes, 25,3711–3720.

Yao,T.,1998:Automatic modeling of(cross)covariance tablesusing fast Fourier transform.Mathematical Geology,30,589–615.

Yilmaz,H.M.,2007:The effect of interpolation methods in surfacedefnition:Anexperimental study.EarthSurface Process and Landforms,32,1346–1361.

:Song,J.J.,S.Kwon,and G.W.Lee,2015:Incorporation of parameter uncertainty into spatial interpolation using Bayesian trans-Gaussian kriging.Adv.Atmos Sci.,32(3),413–423,

10.1007/s00376-014-4040-4.

(Received 04 March 2014;revised 17 June 2014;accepted 22 July 2014)

?Corresponding author:GyuWon LEE

Email:gyuwon@knu.ac.kr