Ordinary kriging and genetic programming for spatial estimation of rainfall in the Middle Yarra River catchment, Australia

Rainfall maps with continuous gridded data provide an important input for many hydrological modelling and water resources management tasks including water budget analysis, flood forecasting, climate change studies, reservoir operation, irrigation scheduling and management, etc. (Moral 2010). The accuracy of most hydrological analysis highly depends on the correct estimation of the spatial distribution of rainfall. Improved estimates of rainfall distribution can be achieved by a dense network of raingauges. Unfortunately, the raingauge network is often sparse and the number of raingauges in the networks is limited due to financial, logistics and geological factors (Abo-Monasar & Al-Zahrani 2014). Consequently, point rainfalls are available from a finite number of raingauges in space. Moreover, the spatial extent of rainfall may not be arranged on a regular grid. Therefore, spatial interpolation methods are important in such cases to estimate rainfall at unsampled locations using the available sampled rainfall values at surrounding locations.

Different deterministic, stochastic and data-driven interpolation methods are generally used to estimate the spatial distribution of rainfall from the ground-based point data. Deterministic methods use the existing configuration of the observed data points to generate a surface, which usually include traditional weighting techniques. Inverse distance weighting (IDW) is a typical representative of deterministic methods, which have been extensively used for spatial interpolation of rainfall (ASCE 1996). Although attempts have been made to upgrade the deterministic techniques over time, limitations of those techniques continue to exist and are discussed in Teegavarapu & Chandramouli (2005).

Kriging-based stochastic interpolation techniques offer important advantages over the deterministic methods, which have been extensively applied for rainfall interpolation (e.g., Moral 2010; Delbari et al. 2013; Abo-Monasar & Al-Zahrani 2014). Kriging has been increasingly preferred because it takes into account the spatial correlation between neighbouring observations and desired location where the estimation is to be made. A noted advantage of kriging over simpler methods such as IDW is that it is capable of utilizing correlated dense secondary variables such as topographic information (e.g., Goovaerts 2000; Feki et al. 2012) to improve the estimation accuracy of sparsely sampled primary variable (i.e., rainfall).

Ordinary kriging (OK) is one of the most preferred stochastic interpolation methods for spatial rainfall estimation from the point rainfall data. However, the accuracy of traditional OK highly depends on the variogram model that generally demonstrates the spatial autocorrelation structure of the underlying stochastic process. Therefore, extra attention should be paid to selecting appropriate variogram models and finding optimal variogram parameters (i.e., nugget, sill, range). The time-consuming and tedious task for calibrating the variogram parameters of the standard variogram models (i.e., spherical, exponential, Gaussian) is a major difficulty associated with the traditional OK.

In the recent past, researchers have demonstrated the application of kriging with different data-driven methods namely, artificial neural networks (ANN) (e.g., Teegavarapu 2007), neuro-fuzzy system and support vector machine (e.g., Kisi & Sanikhani 2015), fuzzy theory and genetic algorithm (GA) (e.g., Chang et al. 2005), combined ANN and GA (Asghari & Nasseri 2015), genetic programming (GP) (e.g., Adhikary et al. 2015b) for spatial interpolation. The major advantage of kriging with data-driven methods is that they usually yield better results than those obtained by traditional deterministic (e.g., IDW) and stochastic interpolation methods (e.g., traditional kriging). However, a major drawback of the ANN- and GA-based methods is that no analytical equations of functional relationships between variables are obtained in these methods, and it is usually not easy to interpret the network weights given by the ANN models (Muttil & Lee 2005).

GP (Koza 1992) is an evolutionary data-driven modelling technique, which can be applied to obtain the variogram model to be used with kriging. The output from GP is an empirical model to approximate the function of the underlying process. GP inferred models exhibit the advantages of creating simple expressions and thereby offering some likely interpretations of the underlying process (Muttil & Lee 2005). The main advantage of GP as compared to traditional models is that it does not assume any a priori functional form of the solution. Globally, a number of applications of GP have been reported in the hydrology and water resources field (e.g., Khu et al. 2001; Sivapragasam et al. 2010, 2012; Yilmaz & Muttil 2014). To the best knowledge of the authors, kriging with a GP-based data-driven method has not been explicitly used as an interpolator in the past to generate gridded rainfall maps. Therefore, the objective of this study is to apply kriging with genetic programming (KGP) for spatial interpolation of monthly and annual mean rainfall to generate the gridded rainfall data set for the Middle Yarra River catchment in Australia. GP is used to get the GP-based non-parametric variogram model to be used in the KGP method, which replaces the standard parametric variogram models in the traditional OK.

The rest of the paper has been organized as follows. First, details of the methodology are presented, which is followed by a brief description of the study area and data sets used. The results are summarized next and, finally, conclusions are drawn.

The methodology of this study consists of three main parts:

• data preparation for normal distribution;

• spatial interpolation of rainfall using deterministic and stochastic interpolation methods;

• assessment of interpolation methods for spatial estimation of rainfall.

The normal distribution of data is a basic requirement of the kriging-based geostatistical approach. Kriging leads to an optimum estimator and yields best results for the normally distributed data (Adhikary et al. 2015a, 2015b). Abo-Monasar & Al-Zahrani (2014) have indicated that variograms developed with the non-normal data are unreliable for kriging interpolation. If data are not normally distributed, appropriate transformations should be implemented to normalize the data. In this study, log-transformation is used for the transformation of non-normal or skewed data. Finally, the widely used Kolmogorov–Smirnov (K-S) test is used to confirm the normality condition of the transformed rainfall data for a 5% significance level. The K-S test is a simple and straightforward statistical test to examine the normality condition of a data set, which has been explained in detail in McCuen (2003). The resulting normalized rainfall data sets are used for interpolation and map production.

In this study, two stochastic interpolation methods (i.e., OK and KGP), and a deterministic interpolation method (i.e., IDW), are used for interpolation of rainfall and preparation of rainfall maps. The basic conceptual difference between traditional OK and KGP methods is presented in Figure 1(a) and 1(b). As can be seen from the figures, the standard parametric variogram models (e.g., spherical, exponential, Gaussian) are used in the traditional OK method, whereas the GP-based non-parametric variogram model is included in the KGP method. A brief description of these methods is presented in the following sub-sections.

The GP-based variogram model is obtained using GPKernel, a software package developed by DHI Water and Environment (Babovic & Keijzer 2000). The GPKernel is a command line-based tool to obtain the underlying functions (in mathematical form) on data. The model with the lowest RMSE and MAE values and the highest CC value is chosen as the best GP model. The KGP method is implemented through a customized spreadsheet application developed in the Microsoft Excel platform for rainfall interpolation to produce gridded rainfall maps.

As p approaches zero and the weights become more similar, IDW estimates approach the simple average of the surrounding observations. However, the effect of the farthest observations on the estimated value is decreased with the increase of ‘p’. The most commonly used value of p is 2 in the IDW method (Isaaks & Srivastava 1989), which is used in this study. IDW method is implemented using GS+ software (Robertson 2008) for rainfall interpolation to generate gridded rainfall maps.

The performance of three different spatial interpolation methods (i.e., OK, KGP and IDW) are evaluated and compared based on the cross-validation test. The cross-validation is a leave-one-out procedure that involves eliminating the rainfall data set individually one by one from the observed data sets and then re-estimating each rainfall data set by using the remaining data sets. Cross-validation test statistics provides important evidence of the performance measures of the interpolation methods. In this study, RMSE, average standardized error (ASE) and root mean square standardized error (RMSS) are used (Adhikary et al. 2015a, 2015b) for cross-validation. They are briefly described below:

• RMSE is used for checking the estimation accuracy between the observed and estimated data. RMSE value close to zero indicates the higher accuracy in estimation.

• ASE is used for checking the unbiasedness condition of the interpolator. ASE value close to zero indicates unbiased estimation of the variable.

• RMSS is used to test the correct assessment of the variability of prediction, which should be close to 1 for correct assessment of prediction variability of the interpolator. RMSS value close to 1 specifies that the variance is consistent and the variability of the estimated data by interpolation is correctly assessed.

There is remarkable variation in the rainfall patterns through different segments of the Yarra River catchment. The mean annual rainfall varies across the catchment from about 1,100 mm in the Upper Yarra segment to 600 mm in the Lower Yarra segment. The Middle Yarra segment (the case study area) covers an area of 1,511 km². There are 19 (numbers 1 to 19 as shown in Figure 2) raingauges in the study area presently operated by the Australian Bureau of Meteorology (BoM). Historical daily rainfall records of all 19 raingauges for 1980–2012 have been collected from the SILO (http://www.longpaddock.qld.gov.au/silo/) climate database and compiled to generate monthly and annual data, which are then used for the analysis. The annual mean rainfall during the aforementioned periods varies between 508.1 mm and 1,913.9 mm, with the highest rainfall in the southern and south-eastern part and lowest rainfall in the north-western part of the study area.

Summary statistics of monthly and annual rainfall data are given in Table 1. September is the wettest month (rainfall amount equals 112.5 mm) with the highest variation in rainfall. The driest month is February (rainfall amount equals 56.4 mm) with the second highest variation. The skewness coefficient computed for monthly and annual rainfall data is presented in Table 1. The data usually fit normal distribution, if the skewness is very small or close to zero (Robertson 2008). Small skewness values in Table 1 indicate that all data sets are almost normally distributed. Furthermore, it is evident from the K-S test results (Table 1) that the null hypothesis cannot be rejected at 5% significance level because the computed K-S test statistic values for all data sets are smaller than the corresponding critical value (i.e., 0.30143). In this study, the null hypothesis is defined as the condition that indicates the data follow normal distribution. Thus, the rainfall data sets are found normally distributed at 5% significance level and can be directly used for variogram modelling and spatial interpolation.

For variogram modelling, initially an experimental variogram is computed for each monthly and annual rainfall data set based on the variogram cloud (Wackernagel 2003). The variogram cloud is a plot in which the semivariance calculated using Equation (3) for any two data points is plotted against the distance between them (Robertson 2008). For 19 raingauges, the variogram cloud for each data set consists of 171 (¹⁹C₂ = 171) data pairs in all directions. However, the variogram clouds are difficult to interpret due to the data points being congested. Hence, the variogram cloud is clustered into different bins to achieve the experimental variogram. A bin is a classification of lags where all lags with similar separation distance and directions are considered in the same bin (Adhikary et al. 2015a). The average lag distance and semivariance for all data pairs in each individual bin gives a discrete point on the experimental variogram plot. The experimental variogram thus consists of a finite number of discrete points with common distance and direction.

In traditional OK, three standard parametric variogram models (i.e., exponential, Gaussian, spherical) are used to fit the experimental variogram for all data sets. The variogram parameters that result in the best fitting performance based on RMSE, MAE and CC measures are chosen as the best fitted model in each case. Figure 3 shows the plots of fitted variogram models for September (taken as an example from all months, which is also the wettest month with high rainfall variation) and annual rainfall. The optimal variogram parameters and fitting performance of the standard variogram models for all monthly and annual rainfall are presented in Table 2. As can be seen from Table 2, the ratio of nugget coefficient to sill is small for all data sets indicating that monthly and annual rainfall is spatially correlated over the study area. This supports the use of kriging methods (i.e., OK, KGP), which inherently consider the spatial correlation. It is also evident from Table 2 that the CC values for the fitted models range from 0.886 to 0.989. Eleven months have CC values higher than 0.9 that indicates strong correlation between monthly rainfall and spatial distribution of raingauges. It can be seen that the Gaussian variogram model is the best fitted model for 6 months (i.e., February, March, April, July, November and December). For the remaining 6 months, the exponential variogram model is the best model for January, June and October while the spherical variogram model gives the best results for May, August and September.

It is important to note that the best fitted variogram model may not always give the best results in interpolation by kriging if it cannot satisfy the cross-validation criteria (Wackernagel 2003). Cross-validation test is often used in kriging to choose the ultimate form of the standard variogram models. Hence, the ultimate form of these models is chosen based on the cross-validation test and the adopted model is finally used to generate the rainfall map. In order to satisfy the cross-validation test criteria, the variogram parameters (i.e., nugget, sill, range) of the standard parametric variogram models are iteratively changed and interpolation results are re-evaluated. The cross-validation results are given in Table 3, which indicates that all adopted variogram models satisfy the cross-validation criteria and are thus appropriate for use in kriging interpolation and rainfall map production.

In the KGP method, the GP-based variogram model is used with traditional kriging. GP is used to approximate the experimental variogram to achieve the GP-based variogram model. Available simple function sets in the standard variogram models (+, −, *, /, x², exp) are used for the GP runs to obtain simple GP models. The evolved GP-based models with the aforementioned function sets can exhibit a similar functional structure to that of the standard variogram models. The best fitted GP model is obtained by using the objective function of minimizing RMSE. In this study, a population size of 500, crossover rate of 0.95, mutation rate of 0.05 and the number of generations as 500 are adopted for the evolution of the best GP-based models. The performance measures (RMSE, MAE, CC) of the GP-based variogram models for all data sets are presented in Table 2. It can be seen from Table 2 that the GP-based variogram model outperforms the standard variogram models for 7 months (i.e., January, May, June, July, September, October, December) and also shows better performance than the other two standard variogram models for another 3 months (i.e., March, April, November) and annual rainfall when all the performance measures are considered. Figure 3 shows the plots of fitted GP-based variogram models for September and annual rainfall. Teegavarapu (2007) indicated that the nugget coefficient of the variogram model significantly impacts the kriging computation and thus should be correctly estimated for best results by kriging interpolation. The nugget coefficient of the GP-based variogram model is determined by extending the function to touch the ordinate for zero distance in the abscissa. The nugget coefficients can then be read from the ordinate, which are also presented for all data sets in Table 2. It can be seen from Figure 3 that the GP-based variogram model has a similar structure to that of the standard variogram models (i.e., spherical, exponential, Gaussian). It is also evident from Figure 3 that the structure of the GP-based variogram model is composed of two standard models. Thus, the GP-based variogram model can be referred to as a combined variogram model of standard structure. A major advantage of the GP-based variogram model is that it does not have any pre-defined structure and parameters unlike the standard parametric variogram models.

Isaaks & Srivastava (1989) mentioned that any form of new variogram function can be fitted to the experimental variogram. They also indicated that developing and selecting such new variogram models (for example, the GP-based variogram model in this study) other than the standard ones may not result in a unique and stable solution for the kriging weights if the ‘positive definiteness condition’, an essential requirement to be satisfied by any variogram model, is not satisfied. The ‘positive definiteness condition’ for a variogram model demonstrates that the variance obtained by kriging interpolation remains always positive (Wackernagel 2003). It is worth mentioning that the standard variogram models always satisfy the ‘positive definiteness condition’ (Teegavarapu 2007) and thus results in a unique solution for the kriging weights in traditional OK. Therefore, the cross-validation test alone is adequate to select the best standard variogram model. However, both cross-validation and positive definiteness criteria must be satisfied simultaneously for the GP-based variogram model so that it can generate a unique solution for the kriging weights in the KGP method. The cross-validation results based on the GP-based variogram model are given in Table 3. It is obvious from Table 3 that all the cross-validation criteria are satisfied by the GP-based variogram model for all data sets. At the same time, the estimated kriging variance (results not shown) for all data sets is found positive in all stations and hence satisfies the positive definiteness condition. Therefore, the obtained GP-based variogram model for all data sets is able to provide a unique and stable solution for the kriging weights and can be used in the KGP method for rainfall interpolation and subsequent map generation.

Different stochastic (i.e., OK, KGP) and deterministic (e.g., IDW) methods are used in the current study to estimate the monthly and annual rainfall patterns in the Middle Yarra River catchment in Australia. In stochastic interpolation framework, kriging with three standard parametric models, namely, OK with spherical model (OK_S), OK with exponential model (OK_E), OK with Gaussian model (OK_G) as well as kriging with GP-based non-parametric variogram model (KGP) are implemented.

The RMSE and cross-validation plot results have been widely used to indicate how accurately an interpolator estimates the observed values. A smaller error in RMSE indicates better prediction. In the case of a cross-validation plot, the perfect estimation is that wherein all scattered points lie on the 45 ° line in the plot with the maximum coefficient of determination (R²) between estimated and measured data.

Table 3 presents the cross-validation results of monthly and annual rainfall estimates over the study area. As can be seen from the table, stochastic interpolation techniques (i.e., OK and KGP) yield better estimation than the deterministic technique (i.e., IDW) for almost all cases considering RMSE values. Similar findings have been reported in many studies (e.g., Goovaerts 2000; Delbari et al. 2013), where IDW has been found to be less accurate than kriging methods for the monthly and annual rainfall estimation.

In this study, the aim is to apply two kriging-based stochastic interpolation methods (i.e., traditional kriging and KGP) and a deterministic interpolation method (i.e., IDW) for the estimation of monthly and annual rainfall distribution in the Middle Yarra River catchment, Australia. The objective is to compare the performance of the aforementioned methods for generating a high quality gridded rainfall data set for the study area in the form of a rainfall map. This study introduces KGP as a new kriging-based interpolation method for interpolation of point rainfalls to generate the gridded rainfall map. In the KGP method, the GP-based non-parametric variogram model is used with kriging instead of using the standard parametric variogram models (i.e., spherical, exponential, and Gaussian). A major advantage of the GP-based variogram model is that it does not require a pre-defined functional form and parameters to fit the experimental variogram unlike the standard parametric variogram models in traditional kriging.

The following points are the major conclusions of this study:

• The cross-validation results indicate that the kriging-based stochastic interpolation methods (i.e., OK and KGP) clearly outperform the deterministic method (i.e., IDW) for spatial interpolation of rainfall.

• Among all the kriging-based methods, OK with the spherical variogram model (OK_S) produces the lowest prediction error and highest coefficient of correlation for rainfall estimation. Thus, OK_S is found to be the most accurate interpolator for estimating the monthly and annual rainfalls over the study area.

• The KGP method yields almost identical prediction error and coefficient of correlation for rainfall estimation to that given by the OK_S for most of the months. Also, the KGP method gives the lower prediction error and is found to be a better interpolator than the OK with the exponential or Gaussian variogram models for most of the months.

• Although for some months, the KGP method cannot produce the lowest prediction error, the method gives reasonable estimates, which are in line with the OK_S (best selected interpolator in this study). The maps generated by the KGP method seem to be consistent with the observed rainfall and reasonable with that produced by the best interpolator. A major advantage of the KGP method is that it eliminates the time-consuming and monotonous task of obtaining the variogram parameters and kriging weights as required with traditional kriging. Therefore, the KGP offers a viable alternative to traditional kriging for spatial interpolation and mapping of rainfall in hydrology.

The authors wish to thank DHI Water and Environment for providing GPKernel and the BoM for providing the necessary data.