Comparison of Spatial Interpolation Methods of Precipitation Data in Central Macedonia, Greece

Athanasios K. Margaritidis

Faculty of Agriculture, Forestry and Natural Environment, School of Agriculture, Department of Hydraulics, Soil Science & Agriculture Engineering, Thessaloniki, Greece.

DOI: 10.4236/cweee.2024.131002                

Abstract

The purpose of this paper is to investigate the spatial interpolation of rainfall variability with deterministic and geostatic inspections in the Prefecture of Kilkis (Greece). The precipitation data where recorded from 12 meteorological stations in the Prefecture of Kilkis for 36 hydrological years (1973-2008). The cumulative monthly values of rainfall were studied on an annual and seasonal basis as well as during the arid-dry season. In the deterministic tests, the I.D.W. and R.B.F. checks were inspected, while in the geostatic tests, Ordinary Kriging and Universal Kriging respectively. The selection of the optimum method was made based on the least Root Mean Square Error (R.M.S.E.), as well as on the Mean Error (M.E.), as assessed by the cross validation analysis. The geostatical Kriging also considered the impact of isotropy and anisotropy across all time periods of data collection. Moreover, for Universal Kriging, the study explored spherical, exponential and Gaussian models in various combinations. Geostatistical techniques consistently demonstrated greater reliability than deterministic techniques across all time periods of data collection. Specifically, during the annual period, anisotropy was the prevailing characteristic in geostatistical techniques. Moreover, the results for the irrigation and seasonal periods were generally comparable, with few exceptions where isotropic methods yielded lower (R.M.S.E.) in some seasonal observations.

Keywords

Interpolation, Kriging, I.D.W., Precipitation, Greece

1. Introduction

Global climate change is altering long-term precipitation patterns, ultimately leading to an increase in both droughts and flood events. As per the Intergovernmental Panel on Climate Change (IPCC) reports from 2022 [1] and earlier (e.g., IPCC 2007) [2] , many changes in the climate system are intensifying directly due to the rising global temperatures. Precipitation is a crucial factor in shaping the climate of a specific region, and having accurate information about precipitation patterns is essential for effective water resource management and is a prerequisite for conducting impact studies related to climate change concerns.

The rainfall regime in the Greek territory shows irregular behavior, both on a spatial and temporal scale regarding the volume and distribution of rainfall. More specifically, the irrigation water pattern shows a strong gradient between western (where the total precipitation is two or three times higher) and the rest of the regions. The transport of moisture from the western to the eastern Mediterranean corresponds to intensification of cyclones in the western Mediterranean and their following eastward motion [3] . Moreover, it is associated to cyclogenesis inside the Mediterranean region [4] . Surveys of rain data [5] [6] [7] [8] [9] for a large number of stations in Greece in time period 1951-1990, during which the data are relatively homogeneous, showed that the rainfall in the Greek area, during the aforementioned 40 years, in some areas shows a statistically significant trend of decrease. Especially in the last 20 years, they show a drier period. Also, a decrease in winter precipitation occurred, although the most statistically significant ones were in the northern and eastern parts as well as the western highlands [10] .

In the Kilkis Regional Unit (R.U.), an investigation into the annual rainfall trends revealed a decreasing pattern at the same monitoring stations over the identical period (1973-2008). These findings were statistically significant at the 5% confidence level, particularly at the Kilkis meteorological station [11] . The climate of the region is characterized as arid and hot, with limited annual precipitation. Doiran Lake in the year of 2000, leading to the conclusion that drought was one of the primary contributing factors to this alteration [6] . Precipitation values in the study area are one of the most basic spatially continuous input data to the various climate models which in turn are useful tools for management and planning of forests, agriculture as well as environmental assessments [12] . Nevertheless, there are not a few times that rainfall “gaps” appear [13] , in the range of the study time series which are due purely to economic reasons (lack of an observer) or other times again to the destruction of the measuring instruments at the observation station and in their untimely replacement. However, the missing values of precipitation in Kilkis R.U. of the geographical division of Central Macedonia (Greece) were fortunately few, while the method of filling the gaps of precipitation was the subject of a previous study [14] . The network of rain gauge stations in the study area appeared quite satisfactory given the high cost of installing a high-density array, combined with any unevenness of the topography.

The ultimate goal of this study is the spatial interpolation of rainfall variability with various deterministic and geostatic techniques [15] - [26] , and the selection of the optimal interpolation based on the smallest root of root mean square error of prediction (R.M.S.E.). Deterministic techniques include inverse distance checks (I.D.W.) as well as radial basis function (R.B.F.) checks, which project spatial distribution surfaces according to the degree of similarity (I.D.W.) or smoothing (R.B.F.) exploiting the known precipitation values of stations as input data. Geostatic techniques on their part (Kriging algorithms) estimate with various assumptions the value of precipitation at an unknown point in the area taking into account the spatial correlation and statistical properties of neighboring known points [27] . Each method has advantages and disadvantages, and the best method is the one that appears to the best fit the characteristics of the input data. The study area of the present research consists of 12 meteorological stations distributed as evenly as possible in the geographical area of Kilkis R.U.

2. Materials and Methods

2.1. Study Area

Kilkis R.U. belongs to the region of central Macedonia, bordering to the north with F.Y.R.O.M., to the south with Thessaloniki R.U., east with Serres R.U., and west with Pella R.U. Based on the available data from the meteorological station of Kilkis, the climate of the region is characterized as dry and hot (semi-arid) with limited annual rainfall.

From a geomorphological point of view, the relief of Kilkis Prefecture (Figure 2), is determined by the morphological characteristics of its individual hydrological basins according to Figure 1. The basins into which the Prefecture is divided are the Doiran lake basin (only a part of the basin from the side of Greece), the Gallikos River basin, the basin of the Axios River and finally the Strimonas River basin (a part of Prefecture of Kilkis) (Figure 1). The data on atmospheric precipitation received by the R.U. of Kilkis have been obtained from 12 specific stations of the study area.

Figure 1. Study area, water basins and the distribution of selected stations (Basemap: Esri & OpenStreetMap).

2.2. Data

The provision of the rain gauge data, which were used for the preparation of this thesis, was made following a relevant request to Kilkis Sanitary Improvement Department, of Y.P.E.K.A. as well as from E.K.V.Y. These data come from the file of the above services, from 12 stations within the study area. Table 1 below shows the meteorological stations that were used in this research, the watersheds in which it is located, its area, the average annual rainfall of each station (time range of 36 years), as well as the altitude of the location of each station from the mean sea level.

The range of observation data as well as station elevations are listed in Table 1. Figure 2 illustrates the location of the study area as well as the locations of rain gauge stations.

The stations have a single chronological beginning of data recording as well as the same chronological end of recording. The number of stations used is the maximum possible for the study area. A preliminary examination of the monthly values showed that the station records contain few measurement gaps. Specifically, in the time series Nov. Doiran 2 missing values appeared in the year 2001 during the months of January and February, one value in the year 2007 in the same time series in the month of July and one more in 2008 in the month of August [14] , which are almost zero percentage of the total range of values (0.08%).

Figure 2. Study area and locations of meteorological stations in Kilkis R.U. (Greece) station Nov. Doiran (F.Y.R.O.M.).

Table 1. Characteristics of weather stations in the study area.

The rain gauge data of the work stations are cumulative monthly values of rain, as recorded at each meteorological station. The data were studied on an annual basis (hydrological year), on a seasonal basis, as well as during the non-rainy-dry-season. A hydrological year is defined as a continuous 12-month period, chosen in such a way that changes, generally in available water supplies, are minimal and the remaining water resources at the end of the period are reduced to a minimum. This period, for the regions of the Northern Hemisphere with a temperate climate, to which Greece also belongs, begins in October of each year and ends in September of the following year. In seasonal values, winter refers to the months of December, January, February, spring to the months of March, April, May, summer to the months of June, July, August and autumn, to the months of September, October, and November. In the hydrological year, two periods can be distinguished: the wet period (October-April), during which the largest amounts of rain are recorded in the Greek area, with the consequence that there is an increasing trend of water reserves, and the dry-rainless period (May-September), by in which rainfall is generally reduced.

2.3. Spatial Interpolation Methods

2.3.1. Deterministic Methods

Inverse Distance Weighted (I.D.W.)

It is applied on the condition that points adjacent to each other show similar precipitation values compared to those located at a far distance [16] . The I.D.W. method in predicting an unknown value of a variable (e.g. precipitation) at a location, it uses the already known values from points around it locally. So any local effect on the value of the unknown point decreases with distance (hence its name, distance weighted inversely). Valid according to the following formula:

$\hat{Z}\left(s_0\right)={\displaystyle {\sum }_{i=1}^N{\lambda }_{i}Z\left(s_i\right)}$ (1)

where: $\hat{Z}\left(s_0\right)$ , the unknown prediction value for a point locally in the region $s_0$ .

N, is the set of points around the unknown prediction point.

λ, the weights for each point N, which decrease with distance.

$\hat{Z}\left(s_i\right)$ , the observation value at locations $s_i$ .

The weights for each known point in space are calculated by the following formula:

${\lambda }_i=\frac{d_{i0}^{-p}}{{\displaystyle {\sum }_{i=1}^N{d}_{i0}^{-p}}}$ ${\displaystyle {\sum }_{i=1}^N{\lambda }_i}=1$ (2)

As the distance increases, the weight decreases from the value p, while as $d_{i0}$ the distance between the prediction points $s_0$ and its neighboring points $s_i$ . The optimal p-value is determined by minimizing the mean square error of prediction. The R.M.S.E. is obtained by calculating the cross validation. Thus, the cross validation value is subtracted from each measurement point and compared to the predicted value at that point. The value of the mean square error of prediction changes as the value of R.M.S.E. also determines the critical value of the variable p. As the distance between the prediction point and the known neighboring points increases, the value of the variable p decreases and the spatial interpolation effect is local [16] in the present study, a value of 2 was used for the p variable, while values greater than 1 are usually used.

Radial Basis Functions (R.B.F.)

Radial Basis Function (R.B.F.) deterministic control is a series of interpolation techniques where the surface of the study area should pass through each measured value of the sample. There are five (5) different basis functions, thin-plate spline, spline, multiquadric function and inverse multiquadric spline. Each of the above basis functions has a different interference surface. Radial Basis Function (R.B.F.) control illustrates one form of artificial neural networks.

R.B.F. method tries through the measured values of the sample on a surface to minimize its total curvature. The way in which the specific adjustment is changed also achieves the best match between the (rainfall) values. The ultimate goal is to minimize the curvature between points which depends on which of the different basis functions will be chosen. The R.B.F. methods are accurate interpolations while differing from global and local interpolations, which appear more imprecise as they do not require the surface to pass through the point values. The difference of the R.B.F. with the I.D.W. method is that the latter will never be able to predict values greater than the maximum measurement value or less than the minimum value. In contrast, the R.B.F. can predict values greater than the maximum sample value and less than the minimum. The optimal value is obtained with the criterion of the cross validation value just like in the I.D.W. method.

Polynomial simulation functions of points in space are characterized by a degree m, while a term r symbolizes the various constraints on the choice of the type of radial basis function. When the value of the parameter r = 0 then there are no restrictions on the function, while when the value of the parameter r = 1 the only restriction is that the function is continuous. When the value of the parameter r = m + 1 then the constraints depend on the degree of the variable m. For value m = 1, r = 2, the radial basis function is called linear spline, for value m = 2, r = 3, quadratic spline and for value m = 3, r = 4 cubic spline. The radial basis functions create smooth projection surfaces, while the results they provide are satisfactory. Finally, because the predictions are exact, the control of radial basis functions is locally sensitive to outliers.

2.3.2. Geostatic Methods

Ordinary Kriging

It follows the following formula:

$Z\left(s\right)=\mu +\epsilon \left(s\right),$ (3)

where μ is an unknown constant mean value and ε is the value of the error between point Z(s) and the constant mean value (Figure 3). It is one of the simplest forecasting methods that have remarkable flexibility. It can also work for trending data while to express autocorrelation between projection points it uses either semivariograms or covariances (they are the mathematical forms that express autocorrelation).

In addition, it can use various transformations, as well as remove the voltage from the focus area with a criterion of the smallest measurement error. It is used for prediction maps as well as error display maps.

Universal Kriging

The geostatic Universal Kriging method follows the following formula:

$Z\left(s\right)=\mu \left(s\right)+\epsilon \left(s\right),$ (4)

where μ(s) some deterministic function, ε the value of the error between the point Z(s) and the deterministic function μ(s).

The autocorrelation in Universal Kriging is modeled by the random errors ε(s). The Universal Kriging method uses either semivariograms or convolutions which are the mathematical forms of expressing autocorrelation and additionally uses various transformation models, which in turn will be mentioned below, from which any trends should be removed, leaving some room for measurement errors (Figure 4).

1) Types of Models

When the experimental semivariogram shows variations (e.g. extreme value phenomena in precipitation values), an adjustment is made to theoretical models [28] [29] which are briefly described below.

1) Spherical model:

$\gamma \left(h\right)=\{\begin{array}{ll}{C}_0+A\ast \left(\frac{3}{2}\left(\frac{h}{r}\right)-\frac{1}{2}{\left(\frac{h}{r}\right)}^3\right)\hfill & \gamma \iota \alpha h\le r\hfill \\ C_0+A\hfill & \gamma \iota \alpha h>r\hfill \end{array}$ (5)

Figure 3. Illustration of the Ordinary Kriging spatial prediction method.

Figure 4. Illustration of the Universal Kriging spatial prediction.

With r the range of influence within which the differences between points are spatially dependent, h the spatial lag interval, C₀, the non-spatial noise or nugget variance, and A the structural variance.

2) Exponential model:

$\gamma \left(h\right)=C_0+A\ast \left(1-{\text{e}}^{\frac{-3h}{r}}\right)\gamma \iota \alpha h>0$ (6)

3) Gaussian model:

$\gamma \left(h\right)=C_0+A\ast \left(1-{\text{e}}^{{\left(\frac{-3h}{r}\right)}^2}\right)\gamma \iota \alpha h>0$ (7)

It is observed that as the values of the spatial lag h increase the value γ(h) increases asymptotically tending to the upper limit called the threshold (partial sill) (Figure 5). In the Gaussian or exponential model the threshold does not meet the asymptote at any point.

Figure 5. Example semivariogram with range, nugget and partial sill.

2) Anisotropy on Models

Anisotropy is the appearance of a change in the properties of a variable in terms of its position but also its orientation in space. The only difference between anisotropic and isotropic models is that the former provides additional information concerning the direction within the parameter boundaries. Anisotropy is a property of a spatial process in which the spatial dependence (autocorrelation) changes with both distance and direction between two locations [30] . This means that the anisotropic model reaches the sill faster in some directions than others. That is, similar things are more similar over long distances in some directions than in others [30] . In the present study and purely for the sake of completeness, the anisotropy was used only in the case of the spatial distribution of precipitation with the Kriging method, while the optimal value of the p parameter was chosen which minimizes the mean squared prediction error of the model. In the same method, the corresponding equations describing the isotropic models were checked.

The parameter values in the anisotropic models in the present study were values for spatial noise (nugget), threshold (partial sill) and range of influence values (Figure 5). Table 2 below shows the various combinations tested in the spatial investigation of precipitation for Ordinary Kriging and Universal Kriging interpolation respectively.

3. Explore Data

Prior to calculating the spatial interpolation methods, the time series data were collected and qualitatively examined. Both the various histograms (Figure 6) and the semivariograms as well as the Q-Q plots, as well as those of the covariances, were a criterion in the separation of spatial interferences with the various spatial interpolation methods. All the above elements had as their ultimate purpose the illustration of the best graphic technique.

The trend was included in some methods (Table 3), while it was not taken into account in some others (Ordinary Kriging).

Figure 6. Histogram of the distribution of precipitation (top) and Q-Q plot curve (bottom) in the annual time series of the study area before transforming data into a Log-transformation.

Table 2. Combinations of models with or without the effect of trend in area for the Ordinary and Universal Kriging spatial controls.

Table 3. Statistical characteristics of precipitation data.

Cross Validation Interpolation

In this work, the results of precipitation in 13 different combinations of spatial interpolations were investigated while comparing them with the criterion of choosing the optimal spatial interpolation based on the smallest Root Mean Square Error of prediction as well as the Mean Error. In all sub-cases, 4 deterministic and 9 geostatic through the software package ArcGis 10, the mean values of the Log time series were selected for each time period separately in the area. The time periods tested with the spatial methods were the hydrological year, the dry season (May-September) or irrigation season, the four seasons as well as the months for 36 hydrological years (1973-2008). There are various validation criteria for spatial interpolation methods, which exhibit sensitivities depending on the input data. Root mean square error as well as mean error were used in this particular study. These two specific criteria were chosen as they are the only ones in the I.D.W. methods and R.B.F. which are included in the spatial interpolation methods of the specific study. They are expressed by the following formulas:

$\text{R}\text{.M}\text{.S}\text{.E}\text{.}=\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(\hat{z}\left(s_i\right)-z\left(s_i\right)\right)}^2/n}}$ (8)

$\text{M}\text{.E}\text{.}={\displaystyle {\sum }_{i=1}^n\left(\hat{z}\left(s_i\right)-z\left(s_i\right)\right)/n}$ (9)

where $\hat{z}\left(s_i\right)$ the normalized Log forecast values and $z\left(s_i\right)$ the recorded normalized values of precipitation for each study time period separately and n the number of intervention points (meteorological stations). The method of spatial interpolation with the smallest values of R.M.S.E. will be the optimal combination.

4. Results and Discussion

In Kilkis R.U., the spatial interpolation of precipitation during the period of the hydrological year was studied with both deterministic and geostatic methods [31] . The average annual rainfall in the 12 meteorological division of Central Macedonia (Greece), as well as in the sole representative of the Doiran basin, station Nov. Doiran of F.Y.R.O.M., (Figure 7) shows in a 3D projection the annual trend at the 12 stations in two directions (north and west) [30] .

An upward trend is observed from south to north and a strong U-shaped trend from west to east. The curves in Figure 7 depict the trend with a manipulation criterion of the average value of the points (stations) of the precipitation in the normalized Log annual time series of the area during the time range 1973-2008.

Table 3 above shows the most important statistics such as the minimum, mean and maximum value of precipitation, the median, the skewness, the kurtosis as well as the standard deviation.

By observing, it is evident that the Log-transformed rainfall series better simulates the Gaussian distribution as the mean value of rainfall is very close to the median, the skewness is close to 0 and the convexity is close to 3 (moderate distribution). In the present study, Log transformed was selected compared to Box-Cox, yielding optimal rainfall prediction values (lower R.M.S.E. values). The parameter λ in the Box-Cox transformations was set equal to 0.5, while in the Log transformations it was set to the value 0.

In general, the anisotropy in combinations of spatial interpolation, Universal and Ordinary Kriging gave the most accurate results and appeared optimal in most study time periods. Figure 10 presents the outcomes related to isotropy during the hydrological year, showcasing the diverse combinations of Universal Kriging. In contrast, Figure 10 (bottom), illustrates anisotropy patterns observed during the irrigation period using the same combinations of Universal Kriging.

Figure 7. Spatial distribution of the trend in the average values of the precipitation of the hydrological year in the 12 meteorological scales of the area under study during the time period 1973-2008.

More in detail, in the annual observations, the deterministic techniques presented a smaller approximation compared to the geostatic ones (Table 4). The average error values showed that only in the I.D.W. method the predicted precipitation value was slightly overestimated compared to the recorded value, while the opposite was found for all other combinations of spatial interpolation techniques. The most reliable according to all the combinations tested was the anisotropic Universal Kriging method [2.3] (Table 4) with constant trend subtraction and the use of the Gaussian transformation. This particular combination yielded the lowest R.M.S.E. of the order of 0.1444 (Table 4), while the isotropic method yielded an R.M.S.E. value of 0.1590 (Figure 8). In the irrigation season (May- September), the anisotropic methods of isotropic isotopes also came to be advantageous.

Table 4. Pooled R.M.S.E. results and M.E. of all combinations of spatial interpolations for the time periods of the year, irrigation period and four seasons.

The anisotropic Universal Kriging method by exponentially removing the stress and by testing the Exponential model emerged as the best with the lowest R.M.S.E. equal to 0.2191 (Figure 9). Same price R.M.S.E. gave the previous combination with a constant trend removal but with a greater deviation in the value of the mean error which was the main reason for its rejection as optimal.

Moreover, in the seasonal observations (winter, summer and autumn) the Ordinary Kriging method and the combination [1.3] of the Universal Kriging method showed the same M.E. values, and R.M.S.E values in the anisotropic method, (Table 4). Finally, in the spring time the anisotropic combination of Ordinary Kriging produced a value of R.M.S.E. equal to 0.2508. Collectively, the results are presented in detail in Table 4, while the most important of them (annual and irrigation season) are shown in Figure 10 below. Annually, the isotropic combinations [1.1], [1.2], [1.3], [2.1], [2.2], [2.3], Universal Kriging and for the irrigation season the same anisotropic combinations.

Figure 11 (down) shows the results of deterministic spatial interpolations (I.D.W. Average, R.B.F. Complete Regularized Spline, R.B.F. Multiquadric) during the annual, irrigation and seasonal periods. In the context of deterministic methods and across all temporal study periods, the Inverse Distance Weighting (I.D.W.) method exhibited the highest Root Mean Square Error (Table 4), whereas the rest deterministic spatial distribution techniques demonstrated R.M.S.E. values approaching zero, with minimal deviations among them (Figure 11, Table 4).

Figure 8. Semivariogram of isotropic combination of spatial interpolation during the hydrological year (U.K. [2.3]).

Figure 9. Semivariogram of anisotropic combination of spatial interpolation during the irrigation period (U.K. [1.2]).

Figure 10. R.M.S.E. in combinations of spatial interpolation of precipitation in the annual (top) and for the irrigated season (bottom) observations.

Figure 11. R.M.S.E. in combinations of spatial interpolation of precipitation in the annual, irrigation and seasonal periods with deterministic methods (top) and geostationary (bottom).

The I.D.W. Average method appears marginally better than the R.B.F. (Complete), except for the autumn season, while the R.B.F., (Multiquadric) shows to deviate from the other two methods, except for the winter and autumn. Furthermore, in the same Figure 11 the isotropy and anisotropy in the spatial interpolation method Ordinary Kriging (constant remove trend) are presented. In most time periods the anisotropy shows optimal values compared to the isotropic method (Ordinary Kriging) (Table 4).

In the below Figure 12 are depicted the most important combinations of spatial interferences in Kilkis R.U. Figure 12 focuses mainly on the geostatic spatial interpolation methods which achieved optimal results compared to the deterministic ones. The stations are depicted, while the distribution of precipitation in the area is show in different colors. Areas with cool colors are characterized as areas with low annual precipitation, while with warm colors (red shades) the stations with the highest.

Figure 13 below shows the most important combinations of spatial interferences in Kilkis R.U., during the irrigation period (May-September) the time range 1973-2008. The validation of each method separately was carried out on the four methods (I.D.W., R.B.F., Ordinary and Universal Kriging) for three types of models (Spherical, Exponential, Gaussian) only for the Universal Kriging method. The inclusion of trend was absent in the Ordinary Kriging method, while constant and exponential trend were respectively incorporated within the Universal Kriging.

Finally, isotropy as well as anisotropy in methods was involved for the 2 aforementioned results based on the M.E. criteria and R.M.S.E., while Figure 14 shows some maps of deterministic methods (I.D.W., I.D.W. Average. and R.B.F. Complete Regularized Spline) in the annual and irrigation precipitation observations in the region. One of the characteristic features in the deterministic spatial interpolation methods are the so-called “bull eyes” around the locations of the rainfall stations as shown in Figure 14.

Figure 12. Spatial interpolation combinations of precipitation over the period of the hydrological year.

Figure 13. Spatial interpolation combinations of precipitation during the irrigation season.

Hydrological Year (Deterministic methods)

I.D.W.

I.D.W. Average Regularized Spline Irrigation Period (Deterministic methods) I.D.W. I.D.W. Average Regularized Spline

Figure 14. Deterministic combinations of annual and irrigation period of spatial interpolation of precipitation.

5. Conclusion

The aim of the present study was to visualize the rainfall spatially at the 12 meteorological stations of the Kilkis R.U. of the geographical division of Central Macedonia (Greece), as well as in one and only station Nov. Doiran on the part of F.Y.R.O.M. The methods used were deterministic and geostatic with the main purpose of comparing the forecast values and presenting the best method. The initial cumulative monthly precipitation values deviated from the model of the standard normal distribution and thus the input data were normalized to the corresponding Log cumulative monthly precipitation values. Log time series of 36 hydrological years beginning on 1973 and ending in 2008 were again input data to the ArcGis 10 software package and evaluated according to some validation criteria. The control criteria for all combinations of the spatial interpolation were the Mean Error (M.E.) and the Root Mean Square Error (R.M.S.E.). The results verified according to the literature the dominance of geostatic methods generally yielded better validation results as it performed the most optimal spatial interpolation methods, e.g. during the irrigation season the combination [1.2] of Universal Kriging. Seasonally, the combination of Ordinary Kriging presented exactly the same results as those of Universal [1.3] for winter, summer autumn, while in spring the optimal combination was that of Ordinary Kriging. Finally, the anisotropic combination [2.3] of Universal Kriging was the best spatial interpolation of the annual observations (Figure 11).

Acknowledgements

I would like to thank the local Department of Land Reclamation of Kilkis, (Ministry of Environment, Energy and Climate Change) and the Greek Biotope/ Wetland Centre for the provision of precipitation data.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Thackeray, C.W., Hall, A., Norris, J. and Chen, D. (2022) Constraining the Increased Frequency of Global Precipitation Extremes under Warming. Nature Climate Change, 12, 441-448. https://doi.org/10.1038/s41558-022-01329-1
[2]	Parry, M.L., Canziani, O., Palutikof, J., Van Der Linden, P. and Hanson, C. (2007) Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, 976.
[3]	Repapis, C., Metaxas, D. and Zerefos, C.S. (1978) Spatial and Seasonal Climatology of Sensible Heat Flux over the Mediterranean Sea. Univ. of Ioannina Tec. Report No. 138.
[4]	Luterbacher, J., et al. (2006) Chapter 1. Mediterranean Climate Variability over the Last Centuries: A Review. In: Lionello, P., Malanotte-Rizzoli, P. and Boscolo, R., Eds., The Mediterranean Climate: An Overview of the Main Characteristics and Issues, Elsevier, Amsterdam, 27-148. https://doi.org/10.1016/S1571-9197(06)80003-0
[5]	Karpouzos, D., Kavalieratou, S. and Babajimopoulos, C. (2010) Trend Analysis of Precipitation Data in Pieria Region (Greece). European Water, 30, 31-40.
[6]	Myronidis, D., Stathis, D., Ioannou, K. and Fotakis, D. (2012) An Integration of Statistics Temporal Methods to Track the Effect of Drought in a Shallow Mediterranean Lake. Water Resources Management, 26, 4587-4605. https://doi.org/10.1007/s11269-012-0169-z
[7]	Dalezios, N. and Bartzokas, A. (1995) Daily Precipitation Variability in Semiarid Agricultural Regions in Macedonia. Greece. Hydrological Sciences Journal, 40, 569-585. https://doi.org/10.1080/02626669509491445
[8]	Νastos, P. and Zerefos, C. (2010) Climate Change and Precipitation in Greece. Hellenic Journal of Geosciences, 45, 185-192. https://doi.org/10.5194/adgeo-25-45-2010
[9]	Κambezidis, H., Larissi, I., Nastos, P. and Paliatsos, A. (2010) Spatial Variability and Trends of the Rain Intensity over Greece. Copernicus Publications on Behalf of the European Geosciences Union. https://doi.org/10.5194/adgeo-26-65-2010
[10]	Feidas, H., Noulopoulou, C., Makrogiannis, T. and Bora-Senta, E. (2007) Trend Analysis of Precipitation Time Series in Greece and Their Relationship with Circulation Using Surface and Satellite Data: 1955-2001. Theoritical and Applied Climatology, 87, 155-177. https://doi.org/10.1007/s00704-006-0200-5
[11]	Margaritidis, A. (2021) Site and Regional Trend Analysis of Precipitation in Central Macedonia, Greece. Computational Water, Energy and Environmental Engineering, 10, 49-70. https://doi.org/10.4236/cweee.2021.102004
[12]	Ranhao, S., Zhang, B. and Jing, T. (2008) A Multivariate Regression Model for Predicting Precipitation in the Daqing Mountains. Mountain Research and Development, 28, 318-325. https://doi.org/10.1659/mrd.0944
[13]	Teegavarapu, R.S.V. (2009) Estimation of Missing Precipitation Records Integrating Surface Interpolation Techniques and Spatio-Temporal Association Rules. Journal of Hydroinformatics, 11, 133-146. https://doi.org/10.2166/hydro.2009.009
[14]	Μargaritidis, A. and Karpouzos, D. (2015) Homogeneity of Rainfall Time Series Analysis in the Wider Area of Doiran Lake. Proceedings of the 9th Panhellenic Conference of the Society of Agricultural Engineers of Greece, Thessaloniki, 8-9 October 2015, 123-132. https://egme.gr/EGME_PRAKTIKA/PDFS/9_conference/EGME%209o%20proceedings.pdf
[15]	Kyriakidis, P., Jinwon, K. and Miller, N. (2001) Geostatistical Mapping of Precipitation from Rain Gauge Data Using Atmospheric and Terrain Characteristics. Journal of Applied Meteorology and Climatology, 40, 1855-1877. https://doi.org/10.1175/1520-0450(2001)040<1855:GMOPFR>2.0.CO;2
[16]	Keblouti, M., Lahbassi, Q. and Hamouda, B. (2012) Spatial Interpolation of Annual Precipitation in Annaba-Algeria—Comparison and Evaluation of Methods. Energy Procedia, 18, 468-475. https://doi.org/10.1016/j.egypro.2012.05.058
[17]	Goovaerts, P. (2000) Geostatistical Approaches for Incorporating Elevation into the Spatial Interpolation of Rainfall. Journal of Hydrology, 228, 113-129. https://doi.org/10.1016/S0022-1694(00)00144-X
[18]	Ikechukwu, M.N., Ebinne, E., Ufot, I. and Ndukwu, I.R. (2017) Accuracy Assessment and Comparative Analysis of IDW, Spline and Kriging in Spatial Interpolation of Landform (Topography): An Experimental Study. Computational Water, Energy and Environmental Engineering, 9, 354-371. https://doi.org/10.4236/jgis.2017.93022
[19]	Stefanidis, S. and Stathis, D. (2018) Spatial and Temporal Rainfall Variability over the Mountainous Central Pindus (Greece). Climate, 6, Article No. 75. https://doi.org/10.3390/cli6030075
[20]	Gofa, F., Mamara, A., Anadranistakis, M. and Flocas, E. (2019) Developing Gridded Climate Data Sets of Precipitation for Greece Based on Homogenized Time Series. Climate, 7, Article No. 68. https://doi.org/10.3390/cli7050068
[21]	Prasad, M.S.G. and Susma, N. (2016) Spatial Prediction of Rainfall Using Universal Kriging Method: A Case Study of Mysuru District. International Journal of Engineering Research & Technology, 4, 1-4. https://www.ijert.org/spatial-prediction-of-rainfall-using-universal-kriging-method-a-case-study-of-mysuru-district
[22]	Philips, D., Dolph, J. and Marks, D. (1992) A Comparison of Geostatistical Procedures for Spatial Analysis of Precipitation in Mountainous Terrain. Agricultural and Forest Meteorology, 58, 119-141. https://doi.org/10.1016/0168-1923(92)90114-J
[23]	Guan, H., Wilson, J. and Makhin, O. (2005) Geostatistical Mapping of Mountain Precipitation Incorporating Autosearched Effects of Terrain and Climatic Characteristics. Journal of Hydrometeorology, 6, 1018-1031. https://doi.org/10.1175/JHM448.1
[24]	Ijaz, H., Gunter, S., Pilz, J. and Yu, H.L. (2010) Spatio-Temporal Interpolation of Precipitation during Monsoon Periods in Pakistan. Advances in Water Resources, 33, 880-886. https://doi.org/10.1016/j.advwatres.2010.04.018
[25]	Javari, M. (2017) Geostatistical Modeling to Simulate Daily Rainfall Variability in Iran. Cogent Geoscience, 3, Article 1416877. https://doi.org/10.1080/23312041.2017.1416877
[26]	Bernardi, M.S., Carey, M., Ramsay, J.O. and Sangalli, L.M. (2018) Modeling Spatial Anisotropy via Regression with Partial Differential Regularization. Journal of Multivariate Analysis, 167, 15-30. https://doi.org/10.1016/j.jmva.2018.03.014
[27]	Wu, H. and Xu, C. (2013) Comparison of Spatial Interpolation Methods for Precipitation in Ningxia, China. International Journal of Science and Research (IJSR), 2, 181-184. https://www.ijsr.net/archive/v2i8/MDIwMTMyMTQ=.pdf
[28]	Cressie, N.A.C. (1993) Statistics for Spatial Data. Wiley Series in Probability and Statistics. Wiley, Hoboken. https://doi.org/10.1002/9781119115151
[29]	Burrough, P.A. and McDonnell, R. (1998) Principles of Geographical Information Systems.
[30]	ESRI Inc., “ArcGIS Help 10,” ESRI Inc., Redlands.
[31]	Hartkamp, A.D., Kirsten, D.B., Stein, A. and White, J.W. (1999) Interpolation Techniques for Climate Variables. Natural Resources Group. https://original-ufdc.uflib.ufl.edu/UF00077518/00001/1j