Precipitation imputation with probability space-based weighting methods

Journal Pre-proofs Research papers Precipitation Imputation with Probability Space-based Weighting Methods Ramesh S.V. Teegavarapu PII: DOI: Reference:

S0022-1694(19)31182-5 https://doi.org/10.1016/j.jhydrol.2019.124447 HYDROL 124447

To appear in:

Journal of Hydrology

Received Date: Revised Date: Accepted Date:

5 August 2019 4 December 2019 9 December 2019

Please cite this article as: Teegavarapu, R.S.V., Precipitation Imputation with Probability Space-based Weighting Methods, Journal of Hydrology (2019), doi: https://doi.org/10.1016/j.jhydrol.2019.124447

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

Precipitation Imputation with Probability Space-based Weighting Methods

Ramesh S.V. Teegavarapu

Professor Department of Civil, Environmental and Geomatics Engineering, Florida Atlantic University, Boca Raton, Florida, 33431 E-mail: [email protected]

1|Page

Abstract Euclidean distance-based spatial interpolation methods that are conceptually simple are commonly used for the imputation of missing precipitation and other hydroclimatic variable datasets. Improved variants of these methods are essential as Euclidean distance will not always serve as an appropriate surrogate that can quantify the inter-gauge relationships in all hydroclimatic and topographical settings. Probability space-based error measure (linear error in probability space (LEPS)) and three distribution similarity hypothesis test statistic values are proposed as surrogates for distances in weighting methods to address this limitation. A k-fold cross-validation exercise of the imputation of precipitation data at 22 rain gauges in a temperate climatic region is used for the evaluation of methods. Improved imputation of missing data was noted from proposed methods compared to that from a distance-based method and is confirmed with statistical inference testing of multiple error and performance measures. Also, the LEPSbased method provided performance measures that are as good as those from a nonlinear optimization method. A local spatial interpolation approach that uses LEPS and two-sample Kolmogorov-Smirnov hypothesis test results quantified as binary outcomes for objective selection of rain gauges are also evaluated. Better estimates of missing data are obtained using this approach compared to those from a Euclidean distance-based global interpolation method. The proposed new probability space-based distances are conceptually superior alternatives to Euclidean distances when used in weighting methods for spatial interpolation. Keywords: missing precipitation, spatial interpolation, distance weighting method, linear error in probability space (LEPS), hypothesis test statistic, Kentucky

2|Page

1. Introduction

Availability of serially complete, continuous precipitation data from rain gauge measurements is not always possible due to malfunctions of instruments, systematic, random and transcription errors, and extreme weather events that disrupt the measurement process itself. Serially complete precipitation data is needed for climate variability and change studies, environmental assessments, and hydrological frequency analyses (Goly and Teegavarapu, 2014; SerranoNotivoli et al., 2017; Aissia et al., 2017). Therefore, interpolation (spatial or temporal) becomes essential to fill any existing gaps. Missing data estimation at a rain gauge is a point estimation process as opposed to a surface generation (i.e., gridded data) exercise (Chen and Liu, 2012; Yeggina et al., 2019). Temporal interpolation for the imputation of missing precipitation data is possible only if high autocorrelation is evident. At low temporal resolutions (e.g., 15 minutes or one hour), rainfall observations often show strong persistence allowing for the use of the last observation carry forward (LOCF) procedure or baseline observation carry forward (BOCF) or single imputation approach (Buuren, 2012). Imputation of missing precipitation at a rain gauge using spatial interpolation is the focus of this study. Deterministic and stochastic spatial interpolation methods using different weighting schemes employing Euclidean (Simanton and Osborn, 1980) and angular (Shepard, 1968; Yeggina et al., 2019) distances have been used in the past for estimation of missing precipitation data. Precipitation data imputation was addressed by many studies using Euclidean distance-based weighting method (e.g., Inverse distance weighting method (IDWM)) (Wei and McGuinness, 1973; Simonton and Osborn, 1980; ASCE, 1996; Garcia et al., 2008). Several variants of IDWM 3|Page

were also developed to improve the interpolation by including a learned search approach (Hodgson, 1989) that reduces the number of distance calculations and topographical information (Shepard, 1968) and anisotropic inverse distance weighting (Tomczak, 1998). Quadrant method (Teegavarapu and Chandramouli, 2005) uses the closest gauge by Euclidean distance to the base gauge (i.e., gauge at which precipitation data is missing or assumed to be missing) in each quadrant whereas normal ratio method (Paulhus and Kohler, 1952; Tung, 1983) uses nearest gauges with similar magnitudes of average annual precipitation totals. Variants of inverse distance and normal ratio methods are introduced by Suhaila et al. (2008) with the incorporation of the strength of correlation between any two sites, and improvements were noted for precipitation imputation compared to the original versions. A recent survey of spatial and temporal interpolation methods for estimation of missing precipitation data by Teegavarapu (2016) serves as a good reference for understanding the limitations and advantages of these methods. Variants of kriging have been used to impute missing precipitation records (Ly et al., 2011; Teegavarapu and Chandramouli, 2005; Teegavarapu 2007). Several drawbacks of the kriging method include difficulty in the selection of the appropriate semi-variogram, uncertainty associated with the assignment of sill and nugget parameters, and the computational effort required to obtain a missing precipitation value at a site. Imputations of precipitation data using artificial neural networks (ANNs) were reported by multiple studies (Kajornrit et al., 2012 and Teegavarapu and Chandramouli, 2005; Teegavarapu et al., 2017). The semi-variogram modeled using ANN in a kriging approach to estimate missing precipitation data by Teegavarapu (2007) is another example use of ANNs. Kim and Pachepsky (2010) documented the use of two stepmethod for constructing daily precipitation data using regression trees and ANNs. Spatially 4|Page

interpolated precipitation estimates were corrected with the help association rule mining technique by Teegavarapu (2009). Bardossy and Pegram (2014) evaluated several methods (regression, kriging, multiple linear regression, and copula-based methods). They noted that copula-based methods performed better than all others. A copula-based method with support vector machines (SVM) for filling missing precipitation data was reported by Landot et al. (2008). However, Teegavarapu et al. (2017) reported that this method did not provide better estimates of missing precipitation compared to those from linear optimization and ANN-based methods. Simple substitution, linear and polynomial regression, and ranked regression approaches were evaluated by Presti et al. (2010) for the imputation of precipitation data, and they concluded that the polynomial regression and simple substitution method provided better estimates of missing data than other methods. A linear programming method using topographical information was reported by Yoo (2013). Pappas et al. (2014) presented a method that uses nearest neighbors in a temporal interpolation approach for the estimation of missing data. Simolo et al. (2010) developed a probability density function-preserving approach to improve the estimation of missing daily precipitation data. According to Teegavarapu (2016) and Brimicombe (2003), most of the spatial interpolation methods lack: (1) a procedure for objective selection of an optimal number of neighbors (or gauges) and neighborhood size (in reference to local or global interpolation); (2) optimal weights in spatial interpolation schemes. Variants of IDWM with the help of optimization formulations were developed by Teegavarapu (2012), which provided mechanisms to select an optimal number of gauges, identify spatial clusters to group rain gauges and use them in interpolation and use the best available gauge in each quadrant for interpolation.

5|Page

Previous studies using spatial interpolation for missing data imputation have proposed and evaluated different conceptually superior alternatives or surrogates for geographical or Euclidean distances. Ahrens (2006) adopted a semi-mean squared difference in non-zero precipitation values between any two rain gauges as a replacement for geographical distance. Teegavarapu and Chandramouli (2005) used the Pearson correlation coefficient (CC) as a replacement for the distance that led to improvements in estimates. Multiple studies (Marteau, et al., 2011; Westerberg et al., 2009; Kim et al., 2008) have reported better estimates of missing precipitation data using this method compared to those from other deterministic methods. Ten real and binary metrics used in numerical taxonomy were adopted as substitutes for Euclidean distances by Teegavarapu (2014a) for imputing missing data. The study documented substantial improvements in multiple error measures compared to those obtained from over a dozen spatial interpolation methods. A modification of IDWM using elevation information of base and other rain gauges by Barrios et al. (2018) provided improved estimates of missing monthly precipitation records in Chile compared to ANN and multiple regression-based methods. Evaluation of variants of Euclidean distance-based method suggests that by use of proper weighting schemes and selection of neighborhood with rain gauges surrounding a base gauge, conceptually superior and simple interpolation methods that are competitive to computationally complex methods can be developed. Exploring parameters that characterize the inter-gauge data relationships better than previously evaluated ones is attempted in this study for developing new variants of distance-based or weighting methods. The contents of this paper are structured as follows. First, the methodology adopted using probability space and hypothesis test statisticbased measures is described. Use of these measures as substitutes for Euclidean distances, weighting method development, and description of the optimization method, objective selection 6|Page

of rain gauges for spatial interpolation are elaborated next. Applications of methods for imputation of missing data, evaluation of methods, and general remarks are presented next. Finally, conclusions from the study are presented.

2. Methodology The methodology adopted for the development of different weighting methods for missing precipitation data estimation is shown in Figure 1. A series of steps illustrated in Figure 1 will show the division of available precipitation data (N days) at all available rain gauges in a region into calibration data (N-M days) and validation data (M days) and estimation of surrogates for Euclidean distances to be later used in a spatial interpolation method. Any surrogate for Euclidean distance should (1) have an ability to characterize the inter-gauge relationships using data; (2) be conceptually simple to estimate, and (3) be robust in handling outliers. Hypothesis test statistic values from distance- and area-based distribution similarity tests that use cumulative distribution functions (CDFs) of datasets such as Kolmogorov Smirnov (KS) (Dodge, 2010), Chi-Square and AD tests can be ideal candidates for the surrogates of Euclidean distances. Probability space-based distance defined by a linear error in probability space (LEPS) score (Potts et al., 1996) is another potential substitute for distance. Details of these surrogates are explained in the next section. The hypothesis test results from the KS, Chi-square, and AD tests can also be used to help in the selection of rain gauges for spatial interpolation. The calibration data is used to derive LEPS score values and test statistic values and also to obtain the hypothesis test results (i.e., null and alternative hypotheses). Considering ns rain gauges and with N observations at each rain gauge, for each site (or rain gauge) assumed to be having missing data, 7|Page

a total of ns-1 values of LEPS score and KS test statistic are derived using the calibration data. The methodology proposed is applicable for estimation of missing precipitation data at a single site (i.e., rain gauge) and data are assumed to missing completely at random (MCAR).

2.1 Surrogates for Distances in Weighting Methods

Four surrogates of distances defined by LEPS score, KS, Chi-square (𝜒2) and Anderson-Darling (AD) test statistic values used in this study are described, and the procedures for estimation of the same and their use in weighting methods are provided in this section. Linear error in probability space (LEPS) score (Potts et al., 1996) is used in this study to define the mean absolute difference between cumulative frequencies (Deque, 2012) of the rain gauge data at any two sites. The calculation of cumulative frequencies based on precipitation data at a base gauge and any other gauge and the representation of the LEPS score is shown in Figure 2. The score 𝛽𝑏,𝑘 in equation (1) is calculated using precipitation data at the base gauge b and those at other gauge k. Variable N refers to the total number of observations, i is the index number for a specific observation and 𝐹𝑏(𝜃𝑏,𝑖) is the empirical cumulative distribution function (ECDF) of observations, and 𝜃𝑏,𝑖 and 𝜃𝑘,𝑖 are the observed precipitation values at a base and another gauge k, respectively. 𝑁

1

𝑁

𝛽𝑏,𝑘 = ∑𝑖 = 1∆𝜖𝑖 = 𝑁∑𝑖 = 1|𝐹𝑏(𝜃𝑏,𝑖) ― 𝐹𝑏(𝜃𝑘,𝑖) |

∀𝑏, ∀𝑘, ∀𝑖

(1)

The LEPS score is calculated using the ECDF of the base gauge only. A low value of LEPS score is an indication of good agreement between observations from the base rain gauge and any 8|Page

other gauge. Reciprocal of the LEPS score can serve as a weight associated with a gauge when used in a weighting method. Another surrogate of distance, a measure of distributional similarity that is used in this study is the test statistic (𝐷𝑏,𝑘) obtained from a nonparametric hypothesis test, a two-sample, two-sided Kolmogorov Smirnov (KS) test (Dodge, 2010). The statistic 𝐷𝑏,𝑘 is also referred to as KS distance or uniform distance (Guidici and Figini, 2009). The statistic is estimated by equation (2) where the variables 𝐸𝐹𝑏 and 𝐸𝐹𝑘 refer to the empirical cumulative distribution functions (CDFs). 𝐸𝐹𝑏(𝜃) is the proportion of 𝜃𝑏 values less than or equal to 𝜃 at base gauge b and 𝐸𝐹𝑘(𝜃) is the proportion of 𝜃𝑘 values less than or equal to 𝜃 at any other gauge k. 𝐷𝑏,𝑘 = max |𝐸𝐹𝑏(𝜃) ― 𝐸𝐹𝑘(𝜃) |

(2)

∀𝑘 , ∀𝜃

𝜃

An alternative to the KS test statistic as a replacement for distance is the test statistic from the Chi-square two-sample test (Press et al., 1992). The Chi-square (𝜒2) test statistic given by equation (3) is also used and evaluated in this study. Precipitation data at base gauge 𝑏 and another gauge 𝑘 is used to calculate the statistic. The variables 𝑧 and z𝑏 refer to the index for the bin (i.e., class interval) and a total number of bins respectively. The number of precipitation values from base gauge 𝑏 within bin 𝑧 is S𝑏𝑧 and the number of values from gauge 𝑘 in the same bin 𝑧 is S𝑘𝑧. 𝜒2𝑏,𝑘

=

𝑏 𝑘 𝑧𝑏 (S𝑧 ― S𝑧 ) ∑𝑧 = 1 𝑏 𝑘 S𝑧 + S𝑧

2

∀𝑏, ∀𝑘, 𝑘 ≠ 𝑏

( 3)

The Anderson-Darling (AD) test statistic (MIL, 2002; Pettitt, 1976; Scholz and Stephens, 1987) is also used as a surrogate of Euclidean distance in this study. Notation for multiple sample AD 9|Page

test is provided here, although the test is used only with two samples (i.e., observations at a base gauge and any other gauge) in this study. The AD test statistic ( 𝐴𝐷𝑏,𝑘) is given by equation (4).

𝐴𝐷𝑏,𝑘 =

(𝑛𝐺𝑙𝑚 ― 𝑛𝑙𝐼𝑚)2 1 𝐿 𝑛 ― 1 𝑙𝑡 ∑ ∑ ℎ 𝑙 = 1 𝑛𝑙 𝑚 = 1 𝑚𝐼𝑚 (𝑛 ― 𝐼𝑚) ― 𝑛ℎ𝑚/4 𝑛2

[

]

(4)

∀𝑏, ∀𝑖, 𝑘 ≠ 𝑏

Where 𝑙𝑡 is the number of groups (samples); 𝑙𝑡 = 1,2,….𝑙𝑡; 𝑛 is the total number of observations (𝑛 = 𝑛1 + 𝑛2 +…𝑛𝑙𝑡), ℎ𝑚 is the number of values in the combined samples equal to 𝑍𝑚, 𝐼𝑚 is the number of values in the combined samples less than 𝑍𝑚 plus one half the number of values in the combined samples equal to 𝑍𝑚, 𝐺𝑙𝑚 is the number of values in the 𝑙-th group which are less than 𝑍𝑚 plus one half the number of values in the combined samples equal to 𝑍𝑚, 𝑛𝑙 is the number of observations in group 𝑙, 𝑍1 , 𝑍2, …….𝑍𝐿 are the distinct values in the combined data set ordered from smallest to largest. The 𝐴𝐷𝑏,𝑘 is estimated with 𝑙𝑡 =2 using observations at base gauge 𝑏 and any other gauge 𝑘. More details of the AD test statistic can be obtained from MIL (2002), Pettitt (1976) and Scholz and Stephens (1987). 2.2 Weighting Method

The imputed precipitation value (𝜃𝑏,𝑖) at base gauge b, in an interval i is given by equation (5). As all available rain gauges are used, the interpolation is considered global. Variable 𝜃𝑜𝑘,𝑖 is the precipitation value at any other rain gauge, k, 𝑤𝑏,𝑘 is the weight for the gauge k when missing data is estimated at base gauge b, and ns is the total number of rain gauges. 𝑛𝑠 ― 1

𝜃𝑒𝑏,𝑖

=

∑𝑘 = 1 𝑤𝑏,𝑘𝜃𝑜𝑘,𝑖 𝑛𝑠 ― 1

∑𝑘 = 1 𝑤𝑏,𝑘

∀𝑏, ∀𝑖, 𝑘 ≠ 𝑏

(5)

10 | P a g e

The weights ( 𝑤𝑏,𝑘) required in Equation (5) are estimated using Equations (6) and (7). The variable 𝑅𝑏,𝑘 can be the Euclidean distance (𝑑𝑏,𝑘) or the surrogate of the distance (𝛽𝑏,𝑘 or 𝐷𝑏,𝑘 or 𝜒2𝑏,𝑘 or 𝐴𝐷𝑏,𝑘 ) with exponent 𝜑 is referred to as friction distance. Generally, a value of 2 is used for 𝜑. However, it can be varied until the best performance from interpolation is achieved, or optimal value for a given dataset can be obtained by solving a mathematical programming formulation (Teegavarapu et al., 2017). 𝑤𝑏,𝑘 =

1

(𝑅𝑏,𝑘)𝜑

∀𝑏, ∀𝑘, 𝑘 ≠ 𝑏, 𝜑 ∈ {1, 2, …10}

𝑅𝑏,𝑘 = 𝛽𝑏,𝑘 or 𝑅𝑏,𝑘 = 𝐷𝑏,𝑘 𝑜𝑟 𝑅𝑏,𝑘 = 𝑑𝑏,𝑘 𝑜𝑟 𝑅𝑏,𝑘 = 𝜒2𝑏,𝑘𝑜𝑟 𝑅𝑏,𝑘 = 𝐴𝐷𝑏,𝑘 ∀𝑏, ∀𝑘

(6) (7)

Methods using Euclidean distance, LEPS score, KS, Chi-square, and Anderson-Darling test statistic values are referred to as ED, LEPS, KS, Chi-Square, and AD methods respectively in this study. 2.3 Objective Selection of Rain Gauges

The selection of rain gauges objectively for imputation is accomplished with the help of the binary variable𝑠 assigned to the outcome of the KS, AD and hypothesis tests. The Chi-square (𝜒2 ) two-sample test with null and alternative hypotheses (i.e., 𝐻0: the two datasets come from a common distribution) and 𝐻1: the two datasets do not come from a common distribution) can be used for the selection of rain gauges for interpolation similar to the use of the KS test. Also, the two-sample version of the AD test with null and alternative hypotheses (i.e., 𝐻0: the populations from which two or more groups of data were drawn are identical) and 𝐻1: the that the populations from which two or more groups of data were drawn are not identical) can be used for selection of rain gauges for interpolation similar to the use of the KS test. 11 | P a g e

For the KS test, if 𝐹𝑏 and 𝐹𝑘 represent the true (but unknown) cumulative distribution functions (CDFs) of observations at the base and other gauge k respectively, then KS test statistic (equation (2)) can be used to determine if independent random samples (data at the base and another gauge) are taken from the same population. The KS test examines two hypotheses: the null hypothesis (𝐻0: 𝐹𝑏(𝜃) = 𝐹𝑘(𝜃)), suggesting two samples are from the same distribution) or the alternative hypothesis (𝐻1: 𝐹𝑏(𝜃) ≠ 𝐹𝑘(𝜃)), suggesting two samples are from different distributions). Testing of the hypothesis can be carried out either by comparing the test statistic ( 𝐷𝑏,𝑘) with a critical value 𝐷𝑐 as indicated by equation (8) or by comparing p-value with a prespecified significance level 𝛼. The result from KS test (i.e., null hypothesis accepted, or null hypothesis rejected) at a pre-specified significance level can be quantified as a binary value using a variable 𝜆𝑏,𝑘 as indicated by inequality 8. 𝑖𝑓 𝐷𝑏,𝑘 < 𝐷𝑐

𝑜𝑟 𝑝 ― 𝑣𝑎𝑙𝑢𝑒 > 𝛼

𝑡ℎ𝑒𝑛 𝐻0 (𝑁𝑢𝑙𝑙 ℎ𝑦𝑝𝑜𝑡ℎ𝑒𝑠𝑖𝑠 𝑖𝑠 𝑎𝑐𝑐𝑒𝑝𝑡𝑒𝑑) 𝑎𝑛𝑑 𝜆𝑏,𝑘 = 1 (𝐴𝑙𝑡𝑒𝑟𝑛𝑎𝑡𝑖𝑣𝑒 ℎ𝑦𝑝𝑜𝑡ℎ𝑒𝑠𝑖𝑠 𝑖𝑠 𝑎𝑐𝑐𝑒𝑝𝑡𝑒𝑑) 𝑎𝑛𝑑 𝜆𝑏,𝑘 = 0

𝑒𝑙𝑠𝑒 𝐻1

(8)

Inclusion of binary variable (𝜆𝑏,𝑘) variable in equation (9) may exclude a subset of all rain gauges from the weighting scheme. 𝑛𝑠 ― 1

𝜃𝑒𝑏,𝑖

=

∑𝑘 = 1 𝑤𝑏,𝑘 𝜃𝑜𝑘,𝑖 𝜆𝑏,𝑘 𝑛𝑠 ― 1

∑𝑘 = 1 𝑤𝑏,𝑘 𝜆𝑏,𝑘

∀𝑏, ∀𝑖, 𝑘 ≠ 𝑏

(9)

This procedure of objective selection of gauges for interpolation is referred to as the SL approach 𝑛𝑠 ― 1

in this study. For any base site if the ∑𝑘 = 1 𝜆𝑏,𝑘 = 0, then imputed data value is obtained out using equation (9) with an appropriate distance. The weights (i.e., 𝑤𝑏,𝑘 values) can be based on Euclidean distance, LEPS score, or any test statistic discussed in the earlier sections.

12 | P a g e

2.4 Optimal Interpolation Method An optimal spatial interpolation method using non-negative least squares (NNLS) formulation is also developed in this study. The formulation, when solved as a constrained least-squares problem, will provide optimal weights required for the imputation of missing data. The formulation with objective function in equation (10) provides non-negative weights due to the inclusion of a constraint given by inequality (11). Minimize

‖𝑬.𝒙 ― 𝑭‖22

(10)

Subject to: 𝑤𝑜𝑏,𝑘 ≥ 0

∀𝑏, ∀𝑘, 𝑘 ≠ 𝑏

(11)

The variables E, 𝒙, and F are vectors. E is the 𝑀 ― 𝑁 x ns-1 matrix of 𝜃𝑀 ― 𝑁, 𝑛𝑠 ― 1 values, 𝒙 defines the ns-1x 1 of 𝑤𝑜𝑏,𝑘 weight values and F refers to M-N x 1 observed precipitation values of data at the base gauge. Equation (10) serves as an objective function that is minimized. A constraint that enforces non-negative weights is included in the formulation using (inequality (11)). The optimal weights obtained from the solution of NNLS formulation are used in equation (11) for the imputation of precipitation data using equation (12). 𝑛𝑠 ― 1

𝜃𝑒𝑏,𝑖 = ∑𝑘 = 1 𝑤𝑜𝑏,𝑘𝜃𝑜𝑘,𝑖

∀𝑏, ∀𝑖, 𝑘 ≠ 𝑏

(12)

This method has been used in past studies (Teegavarapu, 2012, 2014a) and was proven to provide better estimates than kriging and several deterministic spatial interpolation methods. This method is referred to as the OPT method.

13 | P a g e

3. Case Study Application

The methods described in earlier sections are used for the imputation of daily missing precipitation records at several rain gauges located in Kentucky, USA. Rain gauge locations and the topographic details of the region are shown in Figure 3. The case study region (i.e., Kentucky) has a temperate climate dominated by frontal precipitation. The climate for the region is characterized as fully humid, warm with hot summer conditions, as indicated by Kottek et al. (2006), based on the Köppen-Geiger classification scheme. Long-term precipitation data from the calendar year 1971 to 2016 at 22 rain gauges in this region are used for the evaluation of proposed methods. A chronologically continuous (i.e., gap-free) observation set of 16,802 days of precipitation data for all gauges is obtained from the Kentucky Agricultural Weather Center (KAWC). A holdout procedure in which data are randomly assigned to two datasets for calibration (or training) and validation (or testing) is used. The lengths of calibration and validation datasets are 12,602 and 4,200 days, respectively. The validation dataset (4,200 days) constitutes 25% of the entire available data and is randomly selected. The methodology explained in Figure 1 refers to these two datasets by variable notation N-M and M, respectively. Also, to evaluate the robustness of the methods proposed, 35 random 4200-day blocks of data at each gauge are selected and are assuming to be missing. These missing data are estimated using different methods at all gauges. A k-fold validation procedure in which the training dataset (sample) of 12602 days is randomly partitioned into k unequal sized sub-samples with a constant size of validation dataset (i.e., 4200 days) is adopted in the current study. The sub-samples are used for estimation of parameters for the methods, and the validation data is used for testing.

14 | P a g e

4. Measures for Evaluation of Methods

Error measures derived using observed and imputed precipitation data such as mean absolute error (MAE), root mean squared error (RMSE), and a performance measure, Pearson correlation coefficient (CC), are used to evaluate the methods. The measures are calculated based on the validation data. Statistically significant differences in measures between any two methods or among measures derived from multiple executions of the same method are confirmed using parametric and nonparametric hypothesis tests. Notations for parameters used to explain null and alternative hypotheses in these tests refer to population parameters as inferences are made about them from the available sample data. 4.1.

Evaluation of differences in measures

An unpaired two-sample t-test is used to confirm statistically significant differences between error and performance measures obtained from any two methods. The test uses the null and alternative hypotheses: 𝐻0: 𝜇1 = 𝜇2; and 𝐻1: 𝜇1 ≠ 𝜇2 and helps to assess the differences in the mean values (i.e., 𝜇1, 𝜇2) of the measures. As normality of the samples is a requirement for the ttest, a goodness-of-fit test for the Gaussian distribution, the Kolmogorov-Smirnov (KS) test (Massey, 1951; Sheskin, 2011) is used. The t-test also requires evaluation of variances (i.e., 𝜎21, 𝜎22), for which the Fisher test (F-test) (Sheskin, 2011) with null and alternative hypotheses: 𝐻0: 𝜎21 = 𝜎22 and 𝐻1: 𝜎21 ≠ 𝜎22 , respectively is used. A t-test that employs Satterthwaite’s approximation (Satterthwaite, 1946) for the estimation of effective degrees of freedom for situations with unequal variances is also used in this study. 15 | P a g e

4.2.

Evaluation of bi-variate relationships and trends

Spearman’s rho (𝜌) test (Sheskin, 2011) that uses the rank-order correlation coefficient to assess the monotonic association between any two metrics or evaluate trends is used in this study. The null and alternative hypotheses for the test are 𝐻0: 𝜌 = 0 and 𝐻1: 𝜌 ≠ 0, respectively. Spearman’s rho test is used to evaluate the association between LEPS and Pearson’s correlation coefficient and also trends in extreme precipitation indices in this study.

4.3 Evaluation of differences in multiple measures One-way analysis of variance (ANOVA) (Sheskin, 2011) is used to evaluate differences in mean values of more than two performance measures derived from multiple executions of one single method or from multiple methods. The null and alternative hypotheses are 𝐻0: 𝜇1 = 𝜇2 = 𝜇3 …. = 𝜇𝑛𝑝 and 𝐻1: means (𝜇1, 𝜇2…. 𝜇𝑛𝑝) are not all equal. The variable 𝑛𝑝 is the number of groups of performance measures. Homoskedasticity of variances of the measures is evaluated using Bartlett’s test (Bartlett, 1937; Sheskin, 2011) with null and alternative hypotheses as 𝐻0: 𝜎21 = 𝜎22 = ….𝜎2𝑛𝑝 and 𝐻1: 𝜎2𝑙 ≠ 𝜎2𝑚 for at least one pair (l, m) and 𝑙,𝑚 ∈ 𝑛𝑝. A non-parametric version of the one-way ANOVA test, the Kruskal-Wallis H-test (Kruskal and Wallis, 1952; Sheskin, 2011; Corder and Foreman, 2009) was also employed in this study. The p-value (i.e., probability calculated under the null hypothesis, of having an outcome as extreme as the observed value in the sample (Dodge, 2010)) for each test is also reported. 4.4 Nonparametric visual assessment of measures

16 | P a g e

Kernel density estimates (KDEs) (Bowman and Azzalini, 1997) as nonparametric representations of the probability density functions (PDFs) and as replacements for histograms are used to aid comparative visual assessments of error and performance measures. A Gaussian kernel smoothing function with optimal bandwidth that controls the smoothness is used in the current study. 5. Results and Analysis

The possibility of using temporal interpolation for the imputation of missing values was initially evaluated in this study by checking the autocorrelation function (ACF) of precipitation data at several temporal lags at each rain gauge. Figure 4 shows the autocorrelation values at several lags for data at all rain gauges. The upper and lower confidence bounds (95% confidence limits) ( ≈±

1.96 𝑁

) of 0.0178 and -0.0178, respectively, are also shown. Almost all autocorrelation values

at several lags are between these bounds. The ACFs suggest that daily precipitation data constitutes a white noise with serially uncorrelated random values. Based on this observation it can be concluded that temporal interpolation is not an option for the imputation of missing data due to lack of strong autocorrelation at several temporal lags. Imputation of missing data is initially carried out at 22 stations using 4200 days of missing data and later with 35 ensembles of 4200 days randomly selected from the training data without replacement. The friction distance parameter 𝜑 is assigned a value of 2 in this study. A total of 462 LEPS scores (22*21) and Pearson correlation coefficients (CCs) based on data from 22 gauges are obtained. The calculations of LEPS scores are based on n P2 permutations (22 P2) as there are 22 gauges and two gauges are used for one LEPS score estimation. The total number of LEPS scores obtained is equal to 462. It important to note that for LEPS score 𝛽1,2 ≠ 𝛽2,1 as 17 | P a g e

ECDFs of gauge 1 and gauge 2 are used to calculate 𝛽1,2 and 𝛽2,1 respectively. The joint variation of LEPS scores and correlation coefficients is shown in Figure 5. This variation is evaluated as correlations based on observations reflect the inter-gauge relationships. A leastsquares regression line provided points to evidence of a strong linear relationship between LEPS and CC. Statistically significant correlation of -0.963 was observed and was confirmed by the alternative hypothesis (i.e., the existence of non-zero correlation) from Spearman’s’ rho test with a p-value less than 0.0001.

Error and performance measures from the applications of the ED, LEPS and KS methods for imputation of missing data comprising of 4200 days are provided in Table 1. A review of the error and performance measures suggests that LEPS method performs better than the ED method. To confirm that the differences in the error and performance measures from two methods are statistically significant, a parametric hypothesis test, a t-test is used. The normality of RMSE, MAE and CC values are confirmed using a goodness-of-fit test, one-sample KS test, with null hypothesis accepted with p-values of 0.91, 0.46; 0.94, 0.99 and 0.86 and 0.98 and respectively. The F-test confirmed the equality of variances for RMSE values with a p-value of 0.19, inequality of variances for MAE and CC values with p-values of 0.015 and 0.001, respectively. Two types of two-sample t-tests were conducted for equal and unequal variances at a 5% significance level (𝛼). The test results indicate the acceptance of the alternative hypothesis suggesting the mean values of RMSE, MAE, and CC from the ED and LEPS methods are unequal with p-values close to zero. The mean values of RMSE, MAE, and CC for the ED and LEPS methods are 7.46 mm and 6.10 mm, 3.16 mm and 2.59 mm and 0.54 and 0.70 respectively. 18 | P a g e

Error and performance measures based on KS are also provided in Table 1. To confirm statistically significant differences in measures between ED and KS methods, hypothesis tests (i.e., two-sample t-tests) are used. The normality checks for RMSE, MAE and CC values are confirmed using a goodness-of-fit test, one-sample KS test, with the null hypothesis being accepted with p-values of 0.91, 0.93; 0.94, 0.53 and 0.86 and 0.91 and respectively. The F-test confirmed the equality of variances for RMSE values with a p-value of 0.10, inequality of variances for MAE and CC values with p-values of 0.005 and 0.003, respectively. Two types of unpaired two-sample t-tests were conducted for equal and unequal variances at a 5% significance level (𝛼). The test results indicated the acceptance of the alternative hypothesis suggesting the mean values of RMSE, MAE, and CC from the ED and KS methods are unequal with p-values close to zero. The mean values of RMSE, MAE, and CC for the ED and KS methods are 7.46, 6.18, 3.16, and 2.56 mm and 0.54 and 0.69, respectively. Imputation of precipitation data using the Chi-square method obtained from training data, with 9 and 15 bins and a value of 2 for the exponent (𝜑) in equation (5) is carried out. Results from these experiments using two different bin sizes are provided in Table 2, confirming the sensitivity of the statistic to the bin size. Also, the error and performance measures improved as the number of bins increased from 9 to 15. Measures from bin size of 15 are used for further comparison. Performance measures using the Chi-square method were better than those from the ED method. A two-sample unpaired t-test (after conducting normality and Fisher tests) with alternative hypotheses being accepted with p-values close to 0 confirmed that the performance measures from both the methods are different at a statistical significance level of 5%. The error measures from Chi-square and LEPS methods are also compared using t-tests. The t-tests results noted the null hypothesis being true with p-values of 0.97, 0.31, and 0.97 for RMSE, MAE and 19 | P a g e

CC measures respectively, suggesting the Chi-square method provided similar results with no statistically significant differences. Results based on the AD method are also provided in Table 2. The performance measures obtained from this method were inferior to those from LEPS method and are statistically different for RMSE and CC with p-values from t-tests equal to 0.003 and 0.02, respectively. A two-sample, two-sided KS test is used for two specific purposes in this study: (1) to obtain the test statistic (𝐷𝑏,𝑘) and (2) to obtain hypothesis test result with one of the hypotheses (𝐻0 or 𝐻1) being accepted. The two-sided test assumes that two distributions are unequal (i.e., 𝐹𝑏(𝜃) ≠ 𝐹𝑘 (𝜃)). A significance level (𝛼) of 5% is used for hypothesis testing. In this study, the decision to reject the null hypothesis is based on a comparative evaluation of p-value and significance level (α) as suggested by conditional inequality (3) and not based on comparison of the test statistic with the critical value. In the current study, the asymptotic p-value (a p-value that is calculated using an approximation to the true distribution) is used. The use of p-value is mainly due to the algorithmic implementation of test available in the proprietary software used in this study. Since the sample size (i.e., length of training data (N): 12802 days) is large, the asymptotic p-value and exact p-value (a value calculated using the true distribution of the sample) are similar and can provide defensible hypothesis test results. The asymptotic p-value used is considered accurate and valid as the condition: (N1 * N2)/ (N1+ N2) ≥ 4 or (N2/2N) ≥ 4 if N1 =N2, is satisfied with N value equal to the length of the training data (i.e., N-M, 12802 days) for any two rain gauges. The KS test is also used to confirm the distribution similarity between observations at a base rain gauge and any other gauge at a significance level of 5%. For each rain gauge ascertained as a base rain gauge, training data is used to identify the gauges which have similar probability distribution. Table 3 provides the list of gauges with similar distribution as the base gauge based 20 | P a g e

on the KS test. The average number of gauges based on the data is 4, with a maximum of 9 and a minimum of zero. No gauge data with similar distribution as that of gauge K22 was identified based on the two-sample KS test. The groups of gauges having similar distributions as base gauges are also clearly demarcated along with the topographic variations in the study region. A k-means clustering approach (Tan et al., 2014) using a squared Euclidean distance measure is used to isolate four clusters with the help of calibration data. The four clusters identified by the k-means approach are: # 1 (K3, K5, K9, K17); # 2(K4, K8, K10, K20; K12, K16, K18, K19), # 3(K1, K6, K11, K13, K15 and K22) and # 4 (K2, K7, K14 and K21). These clusters seem to align reasonably well according to topographical variations of the region. For each base gauge, one or more gauges identified by distribution similarity using the KS test belong to the appropriate cluster with a base gauge. Figure 6 shows the common gauges identified from the KS test and cluster analysis. Performance measures based on ED, LEPS and KS methods using the SL approach are provided in Table 3. Objective selection of gauges is also possible using Chi-square and AD test results. In the case of the Chi-square test, the similarities in distributions confirmed between the base and any other rain gauge observations were bin size-dependent. Considering the arbitrariness associated with bin size specification, the use of Chi-square test results for a selection of rain gauges is not recommended. Also, the use of results from Chi-square and AD tests did not result in improved estimates of rainfall than those from the KS test in this study. Hence these test results are not used for rain gauge selection in this study. The NNLS mathematical programming formulation in the OPT method is solved using a leastsquares methodology suggested by Lawson and Hanson (1974). The LEPS method provided estimates of missing data as good as those provided by the OPT method evaluated in this study. 21 | P a g e

Performance measures related to these two methods are provided in Table 4. This was confirmed using statistical tests like those that were conducted for comparison of ED and LEPS methods and ED and KS methods. The differences in the mean values of the error and performance measures (i.e., RMSE, MAE, and CC) are not statistically significant and were confirmed using two-sample t-test results. Different values of friction distance (𝜑) were evaluated, and a value of 8 was found to be best for the LEPS method. An optimal value of friction distance can be obtained by using a nonlinear mathematical programming formulation (Teegavarapu, 2013) using the calibration data. A one-way ANOVA test conducted at a 5% significance level is used to evaluate differences in mean values of performance measures at all gauges (given in Tables 1 and 2) from LEPS, KS, Chi-square and AD methods. Measures from 15 bin Chi-square test are used for this evaluation. The test results point to an acceptable alternative hypothesis for measures RMSE, CC, and null hypothesis for MAE with p-values of 0.0002, 0.007 and 0.065, respectively. The mean and median values for RMSE and MAE are higher, and mean and median values for CC for the AD method are lower than those from other three methods. Also, when performance measures from AD are not used, the test resulted in acceptable null hypothesis with p-values of 0.884, 0.587 and 0.872 respectively. Similar results with null (𝐻0: all samples (measures) come from the same distribution) and alternative (𝐻1: not all samples come from the same distribution) hypotheses were also obtained with Kruskal-Wallis tests. Assessments of numerical values of measures and parametric and nonparametric tests used to confirm differences in measures suggests that (1) no statistically significant differences were noted in measures obtained from methods using LEPS, KS and Chi-square statistic values; (2) the KS method provided almost similar performance as that of LEPS method and (3) the median values of performance measures of Chi-square and AD 22 | P a g e

methods were inferior to LEPS and KS methods and (4) LEPS and KS methods provided the best measures among all the methods considering the median values of the measures. The inter-gauge data relationships characterized by correlations and LEPS scores may suggest the interpolation method using correlations as weights (Teegavarapu and Chandramouli, 2005) and a kriging approach can be used for imputation of missing data. Limited experiments using ordinary kriging (OK) for imputation with visual assessment and confirmation of semivariogram fits and with no consideration of anisotropy were conducted. Exponential, Gaussian, and spherical semi-variogram models for OK were evaluated. Results based on CCWM and OK methods for different gauges are provided in Table 5 with one rain gauge chosen from each cluster (shown in Figure 6). In comparison with the correlation coefficient weighting method (CCWM) (Teegavarapu and Chandramouli, 2005; Teegavarapu et al., 2009) and ordinary kriging (OK), the LEPS method proposed in this study performed better when error and performance measures are evaluated. After an exhaustive evaluation of different performance measures, the method based on LEPS was found to be best among all the methods in this study. The KS method was the second-best method. Further evaluations of methods in this study and discussions in the paper are limited mostly to LEPS and in one case KS method. Improved error and performance measures based on the LEPS method in comparison with the ED method were also confirmed using an ensemble of 35 random missing precipitation data blocks of 4200 days at each gauge are used for the imputation of data that are assumed to be missing. A total of 770 values (35 values of each measure x 22 rain gauges) from the ED and LEPS methods are compared using kernel density estimates (KDEs) that use Gaussian kernels and optimal bandwidths. The KDEs for RMSE, MAE, and Pearson’s correlation coefficient (CC) shown in Figure 7 point to improvements in all three measures considering the mean and spread 23 | P a g e

of the values. Shifts in probability densities towards lower magnitudes of RMSE and MAE and higher values for CC confirm the superiority of the LEPS method over the ED method. Similar KDE-based comparisons of performance measures from the ED and KS methods are shown in Figure 8. The length of the historical data may affect the performance of the LEPS method in imputation. The variability in the performances of the LEPS method based on length on the data used to obtain LEPS scores is evaluated by using different lengths of calibration data. Figure 9 shows the variations in RMSE, MAE, and CC values with changes in the length of the data (from 10% to 100% of the training/calibration data). Parametric (one-way analysis of variance (ANOVA)) test and nonparametric hypothesis test (Kruskal-Wallis) are used to assess differences in error and performance measures. Initially, normality of values for each measure is evaluated using onesample KS tests, and the equality of variance (or homogeneity of variance or homoskedasticity) for a sample drawn from normally distributed population is checked using Bartlett’s tests and then the ANOVA test is carried out. The p-values for one sample KS test obtained are 0.73, 0.61, 0.63, 0.49, 0.49, 0.48, 0.47, 0.45, 0.46 for RMSE, and 0.94, 0.68, 0.80, 0.99, 0.99, 0.99, 0.95, 0.97, 0.99 for MAE and 0.99, 0.99, 0.99, 0.97, 0.98, 0.99, 0.96, 0.96, 0.98 for CC respectively. The Bartlett’s test resulted in a null hypothesis being true with p-values of 0.09, 0.06, 0.31, confirming the equality of variances for RMSE, MAE, and CC, respectively. Levene and BrownForsythe tests (Sheskin, 2011) used also confirmed the homogeneity of variances. ANOVA pvalues are equal to 1. The Kruskal-Wallis test also provided a p-value equal to 1 for all the performance measures. Lack of any statistically significant differences in the performance measures based on different lengths of data may be attributed to (1) the representativeness of the data used for calculation of LEPS to the data surrounding temporal intervals in which data is 24 | P a g e

assumed to be missing; (2) stationarity of precipitation time series; (3) time-invariant inter-gauge data relationships and spatial configuration of gauge network. Lack of any trends in WMO (2009) defined eight extreme precipitation indices (viz. Rx1day, Rx5day, CWD, CDD, R10mm, R20mm, SDII, PRCPTOT) derived using all the available historical data at most of the rain gauge sites is also noted. A non-parametric Spearman’s rho test was used at a 5% significance level to assess the trends in these indices. Statistically significant trends are noted at rain gauges 4, 5, 15 for Rx1day, 5, 6, and 15 for Rx5day, 2 and 5 for CWD, 21 for CDD, 15 for R20mm, 5, 8, 13, 15 and 20 for SDII. The minimum and maximum p-values are 0.0006 and 0.99, respectively. Evaluation of ED, LEPS, and KS methods based on complete observed precipitation series and gap-filled series (imputed time series) at 22 gauges using correlation coefficients as shown in Figure 10 indicates that probability space- and KS statistic-based methods perform better than distance-based method at most of the gauges. The median CC values for EP, LEPS and KS methods obtained were 0.92, 0.95, and 0.94, respectively. Observed and estimated values of precipitation at rain gauge K13 using three different methods for 4200 days are shown in Figure 11. While all the three methods, in general, underestimate higher-end precipitation values and overestimate lower-end values, the LEPS-based model provides a lower number of under- and overestimations compared to other methods. In the current study, approximately 25% of the complete dataset is assumed to be missing. In general, site-specific data or distributional characteristics do not change if the threshold value (i.e., the amount of missing data) is below 5%. Appropriate post-correction procedures (i.e., association rule mining, single best estimator and variants of quantile matching approaches) reported in earlier studies by Teegavarapu (2009, 2014b) can be adopted for improved estimates 25 | P a g e

of missing precipitation values and preserve statistical properties of the imputed time series. To reduce the number of overestimations of the lower-end values, especially zero (for no-rain conditions) values, dry-day correction using the nearest neighbor gauge identified by lowest LEPS values derived based on observations at other gauges and the base gauge is adopted. Corrections for spatially interpolated estimates at any base gauge are applied using the condition specified by equation (13). 𝑛𝑛𝑔

𝑖𝑓 ∑𝑛𝑘 = 1𝜃𝑜𝑛𝑘,𝑖 = 0 𝑡ℎ𝑒𝑛 𝜃𝑒𝑐 𝑏,𝑖 = 0, 𝑒𝑙𝑠𝑒

𝑒 𝜃𝑒𝑐 𝑏,𝑖 = 𝜃𝑏,𝑖 , ∀𝑏,∀𝑖, 𝑛𝑘 ≠ 𝑏

(13)

The variables 𝑛𝑛𝑔, 𝑛𝑘 and 𝜃𝑒𝑐 𝑏,𝑖 refer to the maximum number of nearest neighbor gauges used for correction, nearest neighbor gauge, and the corrected spatially interpolated value. In the case of rain gauge K13, the three nearest neighbors identified are K15, K1, and K6 with LEPS values of 0.114, 0.134 and 0.136 respectively. Dry-day corrections when applied to estimates at gauge K13, improved performance measures (i.e., RMSE, MAE and CC values of 4.95 mm, 2.03 mm and 0.774 respectively) were noted, and overestimations are reduced as shown in Figure 12. This correction is not without a limitation, where higher-end extremes are also underestimated. Similar performance measures (i.e., RMSE, MAE and CC values of 4.94 mm, 2.08 mm and 0.774 respectively) were obtained when gauges (K1, K6, K11, K15, and K22) that belong to the same cluster as K13 are used. Quantile-based corrections using historical data at the base gauge is recommended over these corrections even though improved measures are noted. The use of proposed methods to impute missing data at different temporal resolutions is possible. An experiment using LEPS and ED methods was conducted to estimate monthly precipitation data. An improvement of 50%, 45% and 36% in median values of RMSE, MAE and CC respectively for 22 stations were noted for the LEPS method compared to those from the ED

26 | P a g e

method. An exhaustive evaluation of methods proposed in this study needs to be conducted for imputation of missing data at different temporal scales before any general statements can be made about the utility of the proposed methods. 5. Discussion The proposed probability space-based methods in this study eliminate limitations associated with the Euclidean distance-dependent weights in spatial interpolation methods. The use of weights that are representative of the inter-gauge data relationships is beneficial to any spatial interpolation method. LEPS scores and KS, Chi-square and AD test statistic values are sensitive to the variations in inter-gauge data relationships due to climate variability and change. Appropriately chosen temporal windows surrounding the time intervals of missing data for calculation of these scores and statistic values will improve the imputation of data if nonstationarity in precipitation time series is suspected or confirmed. LEPS score calculations require a chronological pairing of data between any two rain gauges, whereas the other three test statistic estimations do not. The latter is less restrictive than the former and is recommended when chronological pairing of data is not possible. The Chi-square test statistic is a competitive alternative to the LEPS score and KS test statistic. One limitation of the Chi-Square test statistic is the subjectivity related to the number of bins used for the calculation of the statistic. While the Anderson-Darling test statistic calculation is objective, it provided inferior performance compared to that of other proposed methods evaluated in this study. Use of probability-space based methods that require LEPS score, KS, Chi-square, and AD test statistic calculations need more computational effort than for the distance-based or optimization methods. The time required for LEPS calculations depends on the length of the historical data and the number of gauges used in interpolation. The methods evaluated in this study are not 27 | P a g e

immune to several limitations of any spatial interpolation approach which include: (1) under- and overestimation of higher-end and lower-end extremes respectively (Teegavarapu, 2013); (2) variance inflation, alteration in the first four statistical moments and autocorrelation, changes in two-state transition probabilities (i.e., probabilities associated with dry or wet state transitions) and extreme precipitation indices when data is imputed using any spatial interpolation method. Statistical characteristics of the precipitation data are altered if the imputed data record length exceeds a specific threshold value. Objective selection of control points for interpolation in moving from global to local interpolation scheme is an issue addressed by Teegavarapu (2013) using mathematical programming models with binary variables. In the current study, an exercise using KS test results as criteria for objective identification of gauges as part of local interpolation has helped in the improvement of the estimates compared to distance-based global interpolation. Conceptually simple interpolation methods of imputation, such as those proposed and evaluated in this study, are easy to implement with minimal computational effort and provide strong support for operational hydrology applications. While the probability-based methods evaluated in this study can be considered as stochastic approaches, uncertainty in the precipitation estimates can only be obtained through the process of multiple imputations (MIs) of missing data using these methods. MI in this context will involve an ensemble of estimates obtained by varying the lengths of the historical data used for deriving LEPS scores or test statistic values (i.e., weights) or local interpolations using sub-sets of all the rain gauges obtained through the evaluations of distributional similarities through two-sample KS tests. In a stationary climate, the uncertainty in the imputed value obtained with weights derived with LEPS or KS statistic will be lower than that from local variants of the proposed methods. The performances of probabilityspace methods can be improved by optimizing the exponents attached to the weights. Also, these 28 | P a g e

methods, when used with appropriate post-imputation correction procedures for estimates, will address most of the limitations of spatial interpolation methods used for imputation of missing precipitation data.

6. Conclusions

Spatial interpolation methods using four new surrogates for Euclidean distances are evaluated for the imputation of precipitation records in this study. The methods using LEPS, two-sample KS, Chi-square (𝜒2), and Anderson-Darling (AD) test statistic values as replacements for distances provided better estimates of missing precipitation records than those from a distance-based method. The LEPS-based method with an appropriate exponent performed as good as an optimal interpolation method when multiple error and performance measures are assessed in this study. This method also outperformed all other three methods that use different test statistic values. Exhaustive assessments using parametric and non-parametric hypothesis tests confirm the superiority of all the proposed methods over inverse distance method and other deterministic, and stochastic interpolation methods with substantial improvements noted in multiple error and performance measures. Objective selection of control points (i.e., rain gauges) for interpolation evaluated using a binary transformation of hypothesis test results from the KS test in one of the methods provided better estimates compared to those from a distance-based method. The proposed data-dependent imputation methods provide alternatives to computationally intensive stochastic interpolation methods. Acknowledgments

29 | P a g e

The data used for analysis in this study is available in the public domain from the Kentucky Agricultural Weather Center (KAWC), University of Kentucky. The author thanks the two anonymous reviewers for their objective criticism of the paper. References Ahrens, B. (2006), Distance in spatial interpolation of daily rain gauge data. Hydrology and Earth System Sciences, 10, 197–208. Aissia, M-A. B., F. Chebana, and T. B. M. J. Ouarda (2017), Multivariate missing data in hydrology- review and applications, Advances in Water Resources, 110, 299-309. ASCE (1996), Hydrology Handbook, 824 pp., Amer. Soc. of Civ. Eng., N. Y Bardossy, A., and G. Pegram (2014), Infilling missing precipitation records – A comparison of a new copula-based method with other techniques, Journal of Hydrology, 519, 1162-1170. Barrios, A., G. Trincado, R. Garreaud (2018), Alternative approaches for estimating missing climate data: application to monthly precipitation records in south-central Chile, Forest Ecosystems, 5(28), https://doi.org/10.1186/s40663-018-0147-x. Bartlett, M.S. (1937), Some examples of statistical methods of research in agriculture and applied biology. Journal of Royal Statistical Society (Supplement) 4, 137–183. Bowman, A. W., and A. Azzalini (1997), Applied Smoothing Techniques for Data Analysis: The Kernel Approach with S-Plus Illustrations. Oxford University Press, Oxford, UK. Brimicombe, A. (2003), GIS, Environmental Modeling and Engineering. Taylor and Francis, London, UK.

30 | P a g e

Buuren, V. S. (2012), Flexible Imputation of Missing Data, CRC Press, Boca Raton, Florida, USA. Chen, F.W. and C.W. Liu (2012), Estimation of the spatial rainfall distribution using inverse distance weighting (IDW) in the middle of Taiwan. Paddy and Water Environment, 10(3), 209-222. Corder, G. W., and D. I. Foreman (2009), Nonparametric Statistics: A Step-by-Step Approach, Wiley. Hoboken, New Jersey. Deque, M. (2012), Deterministic Forecasts of Continuous Variables, in Forecast Verification: A Practitioner’s Guide in Atmospheric Science, I. T. Jolliffe, D. B. Stephenson, editors, 77-94. Wiley. Dodge, Y. (2010), The Concise Encyclopedia of Statistics, Springer. Garcia M, C. D., Peters-Lidard, D. C. Goodrich (2008), Spatial interpolation in a dense gauge network for monsoon storm events in the South western United States. Water Resources Research 44: W05S13, DOI: 10.1029/2006WR005788. Guidici, P., and S. Figini (2009), Applied Data Mining for Business and Industry, Wiley. Guillermo, Q., G. Q. Tabios III, J. D. Salas (1985), A comparative analysis of techniques for spatial interpolation of precipitation. Journal of American Water Resources Association, 21(3), 365-380. Hodgson, M. E., (1989), Searching methods for rapid grid interpolation, Professional Geographer, 411, 51-61.

31 | P a g e

Jolliffe, I. T., and D. B. Stephenson (2012), Forecast Verification: A Practitioner’s Guide in Atmospheric Science, Wiley, Chichester, U.K. Kajornrit, J., W. K. Wong, and C. C. Fung (2012), Estimation of missing precipitation records using modular Artificial Neural Networks (editors) Huang T et al., ICONIP 2012, Part IV, LNCS 7666, pp. 52–59. Kim, J-W., and Y A. Pachepsky (2010) Reconstructing missing daily precipitation data using regression trees and artificial neural networks for SWAT streamflow simulation, Journal of Hydrology, 394(3-4), 305-314. Kottek, M, J. Grieser, C. Beck, B. Rudolf, and F. Rubel (2006), World map of Köppen-Geiger climate classification updated, Meteorologische Zeitschrift, 15(3): 259-263. Kruskal, W. H., and W. A. Wallis (1952), Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621 and errata, 48, 907–911. Kuligowski, R. J., and A. P. Barros (1998), Using artificial neural networks to estimate missing rainfall data. Journal of American Water Resources Association, 34(6), 1437-1447. Landot, T., S. Sgellari, C. Lima, and U. Lall (2008), In-filling missing historical daily rainfall data study, Final Report submitted to South Florida Water Management District, Columbia University, New York, USA. Lawson, C. L., and R. J. Hanson (1974), Solving Least Squares Problems. Prentice-Hall: New Jersey. Ly, S., C. Charles, and A. Degre (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 32 | P a g e

Marteau, R., B. Sultan, V. Moron, A. Alhassane, C. Baron, S. B. Traore (2011), The onset of rainy season and farmer’s sowing strategy for pearl millet cultivation in Southwest Niger, Agricultural and Forest Meteorology, 151, 1356-1369. Massey, F. J. (1951), The Kolmogorov-Smirnov test for goodness-of-fit. Journal of the American Statistical Association. 46(253), 68–78. MIL (2002) Composite Materials Handbook Series, MIL-HDBK-17, Department of Defense Handbook, Army Research Laboratory, Weapons and Material Research Directorate. O’Sullivan, D., and D. Unwin (2010), Geographic Information Analysis, Wiley. Pappas, C., S. M. Papalexiou and D. Koutsoyiannis (2014), A quick gap filling of missing hydrometeorological data, Journal of Geophysical Research: Atmospheres, 1920-1930.

Paulhus, J. L. H. and M. A. Kohler, (1952), Interpolation of missing precipitation records. Monthly Weather Review, 80(8), 129-133. Pettitt, A.N. (1976), A two-Sample Anderson-Darling rank statistic. Biometrika, 63,161-168. Potts, J. M., C. K. Folland, I. T. Jolliffe, and D. Sexton, 1996: Revised “LEPS” scores for assessing climate model simulations and long-range forecasts. J. Climate, 9, 34–53 Presti, R. L., E. Barca, and G. Passarella (2010), A methodology for treating missing data applied to daily rainfall data in the Candelaro River Basin (Italy). Environmental Monitoring and Assessment, 160 (1-4), 1-22. Press, W. H., S. A. Teukolsky, W. T. Vetterlling, and B. P. Flannery (1992), Numerical Recipes in Fortran: The Art of Scientific Computing," Second Edition, Cambridge University Press, pp. 614-622.

33 | P a g e

Satterthwaite, F. E. (1946), An approximate distribution of estimates of variance components, Biom. Bull., 2, 110–114. Scholz, F. W. and M. A. Stephens (1987), K-Sample Anderson-Darling Tests. Journal of the American Statistical Association, 82,918-924. Serrano-Notivoli, R., M. D. Luis, S. Begueria (2017), An R package for daily precipitation climate series reconstruction, Environmental Modelling and Software, 89, 190-195. Sheskin, D. J. (2011), Handbook of Parametric and Nonparametric Statistical Procedures, Chapman and Hall/CRC Press, Boca Raton. Shepard, D., (1968), A two-dimensional interpolation function for irregularly spaced data. Proceedings of the Twenty-third National Conference of the Association for Computing Machinery, 512-524. Simanton, J. R., and H. B. Osborn (1980), Reciprocal-distance estimate of point rainfall. Journal of Hydraulic. Engineering, 106, 1242-1246. Simolo, C., M. Brunetti, M. Maugeri, and T. Nanni (2010), Improving estimation of missing values in daily precipitation series by a probability density function-preserving approach. International Journal of Climatology, 30, 1564-1576. Suhaila, J., M. D. Sayang, A. A. Jemain (2008), Revised spatial weighting methods for estimation of missing rainfall data, Asia-Pacific Journal of Atmospheric Sciences, 44(2), 93104. Tan, P-N., M. Steinbach, and V. Kumar (2014), Introduction to Data Mining, Pearson Education Limited: Essex, UK.

34 | P a g e

Teegavarapu, R. S. V. (2009), Estimation of missing precipitation records integrating surface interpolation techniques and spatio-temporal association rules, Journal of HydroInformatics, 11(2), 133–146, 2009. Teegavarapu, R. S. V. (2007), Use of universal function approximation in variance-dependent interpolation technique: An application in Hydrology, Journal of Hydrology, 332, 16-29. Teegavarapu, R. S. V., M. Tufail, L. Ormsbee (2009), Optimal functional forms for estimation of missing precipitation data, 374(1-2), Journal of Hydrology, 106-115. Teegavarapu, R. S. V. (2012), Spatial interpolation using non-linear mathematical programming models for estimation of missing precipitation records, Hydrological Sciences Journal, 57(3), 383-406. Teegavarapu, R. S. V. (2014a), Missing precipitation data estimation using optimal proximity metric-based imputation, nearest neighbor classification and cluster-based interpolation methods, Hydrological Sciences Journal, 59(11), 2009-2026. Teegavarapu, R. S. V. (2014b), Statistical corrections of spatially interpolated precipitation estimates, Hydrological Processes, 28(11), 3789–3808. Teegavarapu, R. S. V. (2016), Spatial and temporal estimation and analysis of precipitation, Handbook of Applied Hydrology, Edit. V. P. Singh, 2016, McGraw Hill. Teegavarapu, R. S. V., A. Aly, C. S. Pathak, J. Alquist, H. Fuelberg, and J. Hood (2017), Infilling missing precipitation records using variants of spatial interpolation and data-driven methods: use of optimal weighting parameters and nearest neighbor-based corrections, International Journal of Climatology, 38(2), 776-793.

35 | P a g e

Teegavarapu, R. S. V., and V. Chandramouli (2005), Improved weighting methods, deterministic and stochastic data-driven models for estimation of missing precipitation records, Journal of Hydrology, 312(1-4), 191-206. Tomczack, M. (2008), Spatial interpolation and its uncertainty using automated anisotropic inverse distance weighting (IDW) – cross-validation/jackknife approach. Journal of Geographic Information and Decision Analysis, 2, 18-30. Tung, Y. K. (1983), Point rainfall estimation for a mountainous region. Journal of Hydraulic Engineering, 109, 1386-1393 Wei, T. C., and J. L. McGuinness (1973), Reciprocal distance squared method, a Computer Technique for Estimating Area Precipitation, Technical Report ARS-Nc-8., U.S. Agricultural Research Service, North Central Region, Ohio. Westerberg, I., Walther, A., Guerrero, J-L., Coello, Z., Halldin, S., Xu, C-Y., Chen, D., Lundin L-C. (2009). Precipitation data in a mountainous catchment in Honduras: quality assessment and spatiotemporal characteristics. Theoretical and Applied Climatology 101(3-4): 381–296. WMO, World Meteorological Organization. (2009), Guidelines on Analysis of Extremes in a Changing Climate in Support of Informed Decisions for Adaptation. Geneva. Yeggina, S., Teegavarapu, R. S. V., S. Muddu (2019), A conceptually superior variant of Shepard’s method with modified neighborhood selection for precipitation interpolation, International Journal of Climatology, https://doi.org/10.1002/joc.6091 Yoo, J. (2013), Linear programming method considering topographical factors used for estimating missing precipitation. Journal of Hydrologic Engineering, 18(5), 542-551.

36 | P a g e

Figures Precipitation data at all rain gauges

� � � 𝜃����,� 𝜃����,� 𝜃����,� .

� . . 𝜃����,� ��

� 𝜃�,�

Weights Selected gauges for local interpolation

Imputed values at base gauge � 𝜃�,�

......

� 𝜃��,�

. . .

..

..

.. .. ..

� 𝜃��,�

� � 𝜃��,��� �� 𝜃��,� ��

� 𝜃����,�

Spatial Interpolation

� � 𝜃�,� 𝜃�,�

� 𝜃��,�

� � 𝜃�,��� �� 𝜃�,���

...

� 𝜃�,� ��

. . .

� 𝜃�,�

� 𝜃�,���

..

. . .

..

� 𝜃�,�

..

� 𝜃�,�

..

..

� 𝜃�,�

M Observations

.. .. ..

without base gauge data

� 𝜃�,�

Calibration Data � 𝜃�,� �� ��

...

Validation Data

� 𝜃�,�

...

(without replacement)

. . .

� 𝜃�,�

... ... ...

N-M Observations

� 𝜃�,�

..

Random Sampling

..

� 𝜃�,�

...

� � 𝜃��,� �� 𝜃��,�

..

� 𝜃�,�

...

. . .

� 𝜃�,���

...

� 𝜃��,�

: Total number of Observations : Total number of rain gauges M : Number of random samples � 𝜃�,� : Observations assumed to be missing at base gauge � 𝜃�,� :Imputed values using interpolation LEPS : Linear Error in Probability Space KS : Kolmogorov Smirnov 𝜒 2 : Chi-square AD : Anderson-Darling ns

..

.. ...

.. .. .. . . .

� 𝜃�,�

N � 𝜃�,�

� 𝜃�,���

... ... ...

� � 𝜃��,� 𝜃��,�

..

.. � 𝜃�,�

. . .

� 𝜃�,�

...

...

� 𝜃�,�

� 𝜃�,�

...

..

� 𝜃�,�

� 𝜃�,� ��

• Probability space-based distance (LEPS). • Distribution similarity measure : (KS, 𝜒 � and AD test statistic). • Selection of rain gauges using KS test.

� 𝜃�,�

Figure 1. Methodology for imputation of missing data using probability space-based distances and distribution similarity measures

37 | P a g e

1 Δε

Cumulative Probability

ECDF [𝐹𝑏 𝜃𝑏 ] of precipitation data at base gauge

𝜃𝑏: Precipitation at base gauge 𝜃𝑘 : Precipitation at any gauge

0 𝜃𝑘

𝜃𝑏

Precipitation Depth

Figure 2. Schematic representation of error (∆𝜖) calculated for estimation of LEPS using precipitation data from the base and other rain gauges

38 | P a g e

Figure 3. Location of rain gauge stations and topography of the study region.

39 | P a g e

Figure 4. Autocorrelations at different lags for precipitation data at all the gauges.

40 | P a g e

Figure 5. Joint variation of LEPS and correlation coefficients.

41 | P a g e

Figure 6. Clusters of rain gauges based on k-means clustering approach and rain gauges (listed in parenthesis) with similar distribution as that of base rain gauge identified by two-sample KS test.

42 | P a g e

Figure 7. Variations of error and performance measures using KDEs based on two methods based on 35 random missing data blocks of 4200 days each from all the rain gauges. 43 | P a g e

Figure 8. Variations of error and performance measures using KDEs based on two methods based on 35 random missing data blocks of 4200 days each from all the rain gauges. 44 | P a g e

Figure 9. Variations of different measures based on the size of training data used for deriving LEPS scores in LEPS method. 45 | P a g e

Figure 10. Variation of correlation coefficients based on application of different methods using complete observed and filled time series

46 | P a g e

LEPS method

ED method

KS method

Figure 11. Observed and estimated precipitation depths from three methods for rain gauge K13

47 | P a g e

LEPS method

Figure 12. Observed and estimated precipitation depths for rain gauge K13 after multi-gauge dry-day correction

48 | P a g e

Table 1. Error and performance measures based on the application of different methods Method ED LEPS KS Rain RMSE MAE RMSE MAE RMSE MAE CC CC Gauge (mm) (mm) (mm) (mm) (mm) (mm) K1 6.80 2.88 0.627 5.22 2.20 0.791 5.24 2.18 K2 7.49 3.34 0.490 5.07 2.28 0.758 5.15 2.17 K3 6.93 2.78 0.660 6.36 2.68 0.719 6.74 2.75 K4 8.70 4.05 0.289 6.10 2.84 0.614 5.87 2.65 K5 7.53 3.15 0.574 5.54 2.20 0.792 5.25 1.99 K6 7.24 3.08 0.533 6.40 2.75 0.607 7.03 3.01 K7 8.60 3.82 0.375 6.45 2.85 0.632 6.14 2.66 K8 8.49 3.64 0.384 6.15 2.51 0.662 6.76 2.65 K9 7.45 3.05 0.607 5.72 2.37 0.796 5.49 2.24 K10 8.64 3.95 0.293 6.25 2.84 0.588 6.20 2.77 K11 6.99 2.88 0.625 6.42 2.67 0.679 6.57 2.64 K12 6.32 2.52 0.696 6.50 2.69 0.666 6.30 2.57 K13 6.97 2.87 0.603 5.51 2.26 0.756 6.14 2.48 K14 8.09 3.37 0.471 5.95 2.51 0.717 5.98 2.41 K15 6.81 2.70 0.640 5.91 2.43 0.732 6.25 2.55 K16 6.24 2.59 0.747 6.79 2.92 0.691 6.41 2.61 K17 7.55 3.33 0.527 5.85 2.60 0.711 5.96 2.61 K18 7.70 3.09 0.641 8.22 3.39 0.575 7.58 2.97 K19 5.66 2.24 0.774 6.20 2.56 0.720 6.50 2.60 K20 8.55 3.90 0.363 6.00 2.70 0.671 6.33 2.82 K21 8.24 3.40 0.477 6.02 2.41 0.732 5.86 2.31 K22 7.05 3.01 0.555 5.59 2.39 0.704 6.26 2.69 RMSE: root mean squared error, AE: absolute error, CC: correlation coefficient

CC 0.787 0.753 0.676 0.649 0.810 0.527 0.671 0.596 0.808 0.601 0.669 0.692 0.684 0.715 0.695 0.730 0.699 0.656 0.689 0.636 0.748 0.619

49 | P a g e

Table 2. Error and performance measures using Chi-square (𝜒2) method with different bin sizes and AD method

Method 2

Rain Gauge K1 K2 K3 K4 K5 K6 K7 K8 K9 K10 K11 K12 K13 K14 K15 K16 K17 K18 K19 K20 K21 K22

RMSE (mm) 5.45 4.89 6.59 5.57 5.18 6.17 6.13 6.67 5.44 5.91 6.44 6.59 6.4 6.98 6.1 6.07 6.27 7.53 5.82 6.98 6.11 6.14

Chi-square (𝜒 ) Number of bins 9 Bins 15 Bins MAE RMSE MAE CC (mm) (mm) (mm) 2.27 0.766 5.27 2.19 2.1 0.777 4.95 2.15 2.66 0.694 6.54 2.65 2.47 0.687 5.49 2.48 1.91 0.813 5.21 1.96 2.58 0.638 6.4 2.7 2.64 0.674 6.23 2.66 2.64 0.602 6.48 2.59 2.17 0.81 5.46 2.19 2.62 0.631 5.79 2.55 2.59 0.681 6.55 2.65 2.54 0.664 6.6 2.59 2.57 0.654 6.1 2.48 2.74 0.628 5.92 2.48 2.44 0.711 6.1 2.46 2.44 0.762 6.3 2.55 2.64 0.676 6.17 2.61 2.93 0.663 7.49 2.91 2.25 0.759 6.48 2.61 2.76 0.619 6.43 2.77 2.39 0.721 5.99 2.36 2.59 0.635 6.42 2.71 AD: Anderson Darling

AD CC 0.783 0.771 0.698 0.694 0.812 0.608 0.662 0.622 0.810 0.648 0.669 0.664 0.689 0.720 0.711 0.741 0.685 0.667 0.691 0.620 0.734 0.600

RMSE (mm) 6.12 6.01 7.75 6.53 5.64 7.11 6.22 8.45 6.06 6.24 7.24 6.44 8.37 6.24 6.59 8.30 7.10 8.27 6.86 7.14 6.07 7.04

MAE (mm) 2.44 2.41 3.06 2.75 2.03 2.98 2.57 3.17 2.24 2.74 2.87 2.44 3.14 2.50 2.70 2.92 2.97 2.92 2.65 3.11 2.38 3.02

CC 0.727 0.686 0.583 0.600 0.786 0.524 0.672 0.462 0.767 0.594 0.611 0.694 0.470 0.689 0.662 0.629 0.591 0.631 0.663 0.558 0.728 0.523

50 | P a g e

Table 3. Error and performance measures based on the use of rain gauge selection using KS test in different methods

Rain Gauge K1 K2 K3 K4 K5 K6 K7 K8 K9 K10 K11 K12 K13 K14 K15 K16 K17 K18 K19 K20 K21 K22

RMSE (mm) 5.41 4.98 7.40 5.79 5.28 7.53 6.13 6.69 5.72 6.59 6.77 9.12 5.66 6.05 6.36 8.45 6.61 8.45 6.99 6.77 5.95 -

ED MAE (mm) 2.20 2.13 2.89 2.50 1.98 3.09 2.58 2.64 2.34 2.83 2.65 3.42 2.26 2.47 2.48 2.96 2.73 2.96 2.70 2.89 2.31 -

CC 0.770 0.768 0.619 0.668 0.806 0.496 0.679 0.614 0.785 0.571 0.657 0.478 0.740 0.706 0.697 0.623 0.646 0.623 0.656 0.608 0.739 -

Method LEPS RMSE MAE (mm) (mm) 5.37 2.19 4.98 2.12 7.17 2.83 5.59 2.40 5.24 1.96 7.26 2.91 6.08 2.55 6.48 2.54 5.56 2.26 6.35 2.73 6.49 2.56 9.12 3.42 5.43 2.15 5.89 2.41 5.91 2.29 8.45 2.96 6.27 2.66 8.45 2.96 7.03 2.72 6.62 2.80 5.87 2.27 -

CC 0.774 0.769 0.637 0.690 0.809 0.541 0.684 0.636 0.799 0.597 0.680 0.478 0.762 0.724 0.739 0.623 0.677 0.623 0.654 0.624 0.747 -

RMSE (mm) 5.37 5.43 7.40 5.73 5.24 7.43 6.27 7.05 5.55 6.43 6.92 9.12 6.22 6.32 6.86 8.45 6.33 8.45 7.32 6.61 5.88 -

KS MAE (mm) 2.18 2.22 2.91 2.46 1.91 3.04 2.62 2.72 2.21 2.76 2.72 3.42 2.46 2.49 2.62 2.96 2.71 2.96 2.78 2.80 2.28 -

CC 0.774 0.732 0.616 0.678 0.809 0.505 0.662 0.573 0.798 0.588 0.645 0.478 0.678 0.684 0.656 0.623 0.665 0.623 0.630 0.626 0.745 -

51 | P a g e

Table 4. Error and performance measures based on the LEPS method with an exponent of 8 and OPT method

Method Rain Gauge K1 K2 K3 K4 K5 K6 K7 K8 K9 K10 K11 K12 K13 K14 K15 K16 K17 K18 K19 K20 K21 K22

RMSE (mm) 5.08 4.84 5.96 5.10 5.07 6.00 6.01 5.60 5.13 5.56 6.18 6.17 5.28 5.49 5.53 5.94 5.67 7.34 5.56 5.12 5.47 5.20

LEPS* MAE (mm) 2.08 2.10 2.39 2.15 1.87 2.34 2.56 2.19 2.05 2.36 2.47 2.38 2.02 2.23 2.10 2.36 2.48 2.81 2.11 2.20 2.08 2.13

CC 0.800 0.785 0.758 0.749 0.823 0.681 0.691 0.730 0.836 0.689 0.708 0.712 0.778 0.766 0.770 0.774 0.733 0.683 0.782 0.773 0.785 0.751

RMSE (mm) 5.04 4.70 5.90 4.91 5.12 5.48 5.89 5.44 4.98 5.34 6.09 5.91 5.12 5.22 5.31 5.92 5.51 7.27 5.43 5.03 5.32 4.98

OPT MAE (mm) 2.10 2.01 2.35 2.04 1.89 2.02 2.39 2.07 2.02 2.14 2.37 2.23 1.99 2.08 2.08 2.38 2.45 2.68 2.13 2.08 1.99 1.94

CC 0.805 0.798 0.763 0.762 0.818 0.721 0.705 0.746 0.841 0.704 0.719 0.735 0.793 0.794 0.790 0.776 0.751 0.693 0.797 0.783 0.797 0.773

*with an exponent value (friction distance) of 8 for the weights, OPT: Optimization method

52 | P a g e

Table 5. Error and performance measures based on correlation coefficient weighting method and ordinary kriging

Method Rain Gauge K3 K19 K15 K2

RMSE (mm) 6.48 6.30 6.04 5.15

CCWM MAE (mm) 2.75 2.61 2.51 2.33

CC 0.706 0.700 0.718 0.750

RMSE (mm) 7.24 6.36 6.44 5.07

OK MAE (mm) 2.94 2.64 2.58 2.13

CC 0.635 0.702 0.687 0.754

CCWM: Coefficient of Correlation Weighting Method; OK: Ordinary Kriging

53 | P a g e

•

New surrogates for Euclidean distances in weighting methods for imputation of precipitation data

•

Probability space-based error and hypothesis test statistic measures in spatial interpolation

•

Objective selection of rain gauges for local interpolation of precipitation data

54 | P a g e