Alternative approaches for estimating missing climate data:application to monthly precipitation records in South-Central Chile

Alonso Barrios,Guillermo Trincadoand René Garreaud

Abstract Background:Over the last decades interest has grown on how climate change impacts forest resources.However,one of the main constraints is that meteorological stations are riddled with missing climatic data.This study compared five approaches for estimating monthly precipitation records:inverse distance weighting(IDW),a modification of IDW that includes elevation differences between target and neighboring stations(IDWm),correlation coefficient weighting(CCW),multiple linear regression(MLR)and artificial neural networks(ANN).Methods:A complete series of monthly precipitation records(1995–2012)from twenty meteorological stations located in central Chile were used.Two target stations were selected and their neighboring stations,located within a radius of 25 km(3 stations)and 50 km(9 stations),were identified.Cross-validation was used for evaluating the accuracy of the estimation approaches.The performance and predictive capability of the approaches were evaluated using the ratio of the root mean square error to the standard deviation of measured data(RSR),the percent bias(PBIAS),and the Nash-Sutcliffe efficiency(NSE).For testing the main and interactive effects of the radius of influence and estimation approaches,a two-level factorial design considering the target station as the blocking factor was used.Results:ANN and MLR showed the best statistics for all the stations and radius of influence.However,these approaches were not significantly different with IDWm.Inclusion of elevation differences into IDW significantly improved IDWm estimates.In terms of precision,similar estimates were obtained when applying ANN,MLR or IDWm,and the radius of influence had a significant influence on their estimates,we conclude that estimates based on nine neighboring stations located within a radius of 50 km are needed for completing missing monthly precipitation data in regions with complex topography.Conclusions:It is concluded that approaches based on ANN,MLR and IDWmhad the best performance in two sectors located in south-central Chile with a complex topography.A radius of influence of 50 km(9 neighboring stations)is recommended for completing monthly precipitation data.

Keywords:Climatological data,Cross-validation,Artificial neural networks,Multiple linear regression

Background

The effects of climate on natural resources have become highly relevant(Cannell et al.1995).In forestry,there is an increasing interest to study the influence of climate on forest productivity(álvarez et al.2013),forest hydrology(Dai et al.2011),soil water availability(Ge et al.2013),and wood quality(Xu et al.2013).Nowadays,climate data are also required for parameterizing process-based simulators of tree growth(Sands and Landsberg 2002)and for studying forest water balance(Huber and Trecaman 2002),phenology processes(Codesido et al.2005)and to carry out pest and disease research(Ahumada et al.2013).To perform these studies,complete and homogenous climate data that covers a sufficiently long period of time is required(Teegavarapu 2012;Khosravi et al.2015).

Climate data often have missing information that limits their use(Alfaro and Pacheco 2000).Missing values in climate series affects parameter estimation when applying regression and multivariate analysis techniques(Ramos-Calzado et al.2008).In most cases,some techniques must be applied to estimate missing data.In forestry,there are few studies that have compared the accuracy of different approaches.Furthermore,factors that might affect their precision have not been studied in detail.

The simplest approach for imputing missing values involves the data being filled-in.The main limitation is that these approaches are suitable for small gaps and can only be applied to climate variables with a high degree of autocorrelation(Khosravi et al.2015),which is not the case for annual mean temperatures or precipitation values.A more common approach to complete missing data is to use information from neighboring meteorological stations(Vasiliev 1996),using techniques such as inverse distance weighting(IDW).Nonetheless,horizontal distance is not a measure of spatial autocorrelation(e.g.,Ahrens 2006;Ramos-Calzado et al.2008),especially when the region contains prominent topographic features or major water bodies.Indeed,two relatively close stations can feature substantial differences in their mean climate and climate variability if they are located at opposite sides of a mountain range.Spatial correlations could be quantified by calculating the correlation coefficient between time series obtained at different locations.Teegavarapu and Chandramouli(2005)found that replacing distances with correlation coefficients as weights improved estimation of missing precipitation data.The resulting method is known as a coefficient of correlation weighting(CCW),reported by Teegavarapu(2009).

Simple and multiple linear regressions have been successfully used to estimate precipitation(Pizarro et al.2009),and temperature(Xia et al.1999)in different topographical conditions.Alfaro and Pacheco(2000)compared different estimation approaches for missing precipitation data,including normal ratio and linear regression.They found that the best results were obtained when applying multiple linear regression;in agreement with the results reported by Xia et al.(1999)and Pizarro et al.(2009).

Recent studies used artificial neural networks for completing climate data(Kuligowski and Barros 1998;Khorsandi et al.2011;Ghuge and Regulwar 2013).Kuligowski and Barros(1998)compared the performance of artificial neural networks for completing six-hour precipitation data at six test stations from nearby stations,to four other approaches,such as the simple nearest-neighbor estimate,the arithmetic average,the inverse distance weighting and linear regression.They found that artificial neural networks and linear regression approaches produced the lowest overall errors.Khorsandi et al.(2011)compared four approaches including the artificial neural network,normal ratio,inverse distance weighting,and a geographical coordinate approach for completing missing monthly precipitation data.They found that artificial neural networks produced the best results compared to other approaches.Different artificial neural network designs have been developed and tested for missing data estimation.Coulibaly and Evora(2007)compared six different types of artificial neural networks and found that the multilayer perceptron(MLP)appears to be the most effective for completing missing daily precipitation values and missing daily maximum and minimum temperature values.

Several studies evaluated the predictive capability of different approaches for completing missing climate data,but few have evaluated the effects of the radius of influence when selecting neighboring stations(e.g.Chen and Liu 2012)in regions with complex topography.We tested the predictive capability of five reported approachesatcompletingmissingdataof monthly precipitationsfrom 1995 to 2012 from south-centralChile (around 37°S)along thewest slope of the Andes mountain range.This region features a climate transition between semiarid conditions in the north and more humid conditions in the south(e.g.Viale and Garreaud 2015).More importantly,the region exhibits a complex topography including a central valley flanked by the Andes mountain range,reaching over 2.000 m asl(above sea level).Our specific objectives are(i)to compare different approaches forestimating missing monthly precipitation data based on measures of precision and bias,and(ii)to evaluate the effect of the number of available neighboring stations within a radius of influence(25 and 50 km)on estimation precision.We selected monthly precipitation as the target climate variable because it is a limiting factor for fast-growing radiata pine plantations in Chile(Gerding and Schlatter 1995;álvarez et al.2013).

Methods

Data

Twenty meteorologicalstations from the DGA(Dirección General de Aguas)located in central Chile(36°–38°S,71°–72°W,Fig.1)with complete monthly precipitation records from January 1995 to December 2012 were selected(Table 1).Annual mean rainfall in this region ranges from 1000 to 2000 mm.

We selected stations Diguillin(number 1)and Mulchen(number 11)as target stations because they were surrounded by an equal number of neighboring stations with a radius of influence of 25 and 50 km(Table 1).Meteorological stations 1 and 11 were located in the Andean foothills at an elevation of 670 m asl and in the Central valley at an elevation of 130 m asl,respectively(Fig.1).In this part of the country there is marked seasonality,with winter(May to September)rainfall accounting for over 65%of annual accumulation and associated with widespread frontal systems crossing the region(e.g.Falvey and Garreaud 2007).Episodes of isolated convection are infrequent over this region and account for a very small fraction of the annual accumulation (Viale and Garreaud 2014).However,winter frontal rainfall is modified by the topography producing a marked enhancement over the western slope of the Andes relative to low-land values(Viale and Garreaud 2015).For instance,the horizontal distance between our target stations is less than 70 km but annual mean precipitation increases from about 1200 mm in the lower station(1)to 2100 mm in the higher station(11).On the other hand,annual precipitation across central Chile exhibits significant inter-annual variability where the standard deviation of annual accumulation is up to a third of the mean value due to the effects of the cold and warm phases of El Ni?o Southern Oscillation(ENOS;e.g.,Montecinos and Aceituno 2003;Garreaud 2009).

The meteorological stations located in the Andes foothills show less variability in terms of mean annual precipitation than stations located in the Central valley(Fig.2).This can be partially explained by an increased amount of stations at higher elevations located in the Central valley(CV=62.8%)compared to the Andes foothills(CV=54.5%).

The Euclidean distance between target and neighboring stations were computed using the formulawhere xmand ymare the UTM coordinates of the target station and xiand yiare the UTM coordinates of the neighboring station.The radius of influence of 25 km included three neighboring stations and the radius of influence of 50 km included nine neighboring stations for each target station(Fig.1).Although the neighboring station 20 was 52 km away from target station 11 it was maintained in the analysis in order to have the same number of neighboring stations for each target station(Table 1).

Fig.1 Geographical location of the meteorological stations.Target stations are represented with a red dot()and neighboring stations with a black dot(●).Radius of influence of 25 km(dashed circle)and 50 km(continuous circle)

Table 1 Meteorological stations used for estimating monthly precipitation values

Minimum station density guidelines for different climatic and geographic zones have been established by the World Meteorological Organization(WMO 2008).In the study area,the corresponding network density of meteorological stations is～1.3 stations per 1000 km2,which is less than the minimum recommended network density for mountainous areas(4 stations per 1000 km2).Because existing network of climatological stations has a low density to explain the spatial variability of rainfall in mountainous regions at shorter time scales(e.g.hourly and daily)we used a monthly timescale for performing a comparison of approaches for estimating missing monthly precipitation data.Longer-timescales rainfall(e.g.monthly,seasonaland annual)tend to bemore spatially homogeneous than shorter-timescales rainfall(Cheng et al.2008;Girons-Lopez et al.2016).In addition,longer-timescale rainfall is of major importance for the evaluation of water availability for management of forest plantations(álvarez et al.2013).

Fig.2 Variation of the annual precipitation sum for meteorological stations at the Andean foothills(a)and Central valley(b).Meteorological stations are arranged by elevation,the target stations 1 and 11 are highlighted in dark gray.Red dotted lines indicate the inter-annual precipitation mean

Approaches for estimating missing data

We selected the following five reported approaches for estimating missing monthly precipitation data for the two target meteorological stations.All approaches were implemented and tested using the Statistical Analysis System-SAS(SAS Institute Inc.2009).

Inverse distance weighting(IDW)

Missing data from target station m are determined from the values observed in neighboring stations weighted by the inverse distance between the target and the neighboring stations.The missing data yj(m)at station m,based on the values observed in neighboring stations is given by,

where,n is the number of neighboring stations with information from the month to be estimated,dmiis the Euclidian distance between station i and m,and xj(i)is the observed value at station i,and k is the distance of friction ranging from 1 to 6(Vieux 2004).In this study,we used a value of k=2 suggested by Teegavarapu(2009).

Modified inverse distance weighting(IDWm)

Elevation has an important influence on precipitation(Golkhatmi et al.2012;Viale and Garreaud 2015),therefore we used the elevation differences between the target and neighboring stations to adjust IDW estimates.A revised version of the approach proposed by Chang et al.(2005)ensuring that the sum of the weights equals 1 was used.This approach considers not only the effect of Euclidian distances but also differences in elevation.Elevation differences were added to the base IDW formula as;

where hmiis the absolute elevation difference between the target and neighboring stations,and exponent a is a power parameter.Thus,hmimodifies the weights of IDW,prioritizing neighboring stations that are at the same or a close elevation of the target station giving them higher weights during the calculations.Values of the exponents a and k between 1 and 3 were tested,and a value of a=1 and k=1 were selected for computing the missing data.

Correlation coefficient weighting(CCW)

In this approach distance is replaced by Pearson’s correlation coefficients.The missing value j in a given month at the target station m is completed as,

where rmiis the Pearson’s correlation coefficient between the precipitation series of the neighboring station i and the incomplete series of the target station m,xj(i)is the monthly value observed at station i(Teegavarapu 2009).

Multiple linear regression(MLR)

The ordinary least squares method is used to fit a line between the observed data from the target station and several neighboring stations.We used a stepwise selection process to ensure that each station in the final linear model contributes to the accuracy of the estimate without compromising the goodness of fit.The linear model has the following form,

where yj(m)is the observed monthly value from the target station m,xj(i)is the observed value in the neighboring station i and βiare the parameters to be estimated(Freund et al.2006).

Artificial neural networks(ANN)

An artificial neural network is a computational model inspired structurally and functionally in biological neural networks(Coulibaly and Evora 2007).The architecture of the designed artificial neural network corresponds to a feed forward multilayer perceptron with one hidden layer with ten neurons(see e.g.Dreyfus 2005;Teegavarapu and Chandramouli 2005).The observed values in the neighboring stations are used for the input layer and the estimated values for the target station are obtained for the output layer.To model the transformation of values through the layers a sigmoid function was used for the hidden layer and linear activation was used for the outer layer.Training of the artificial neural network was performed by using the standard error as criterion,applying the Levenberg-Marquardt training algorithm(Khorsandi et al.2011;Ghuge and Regulwar 2013).The artificial neural network was built,trained and simulated using the SAS NEURAL procedure(SAS Institute Inc.2009).

Cross-validation and statistical evaluation

Because complete monthly precipitation records were available for all meteorological stations,we simulated missing values using cross-validation for evaluating the accuracy of the estimation approaches.Cross-validation is a technique used for assessing how generalized the results of a statistical analysis are compared to an independent dataset(Chen and Liu 2012).For each target station,data were randomly partitioned into 10 nearly equally sized folds containing 21 or 22 monthly precipitation records(about 10%of total data).Subsequently,10 estimation and validation iterations were performed,where 9 folds were used to estimate model parameters and the remaining fold was used to validate the method.Refaeilzadeh et al.(2009)reported that 10 folds are the most common because it allows estimations to be made with 90%of the data,producing representative data.

The performance and predictive capability of the approaches for completing missing monthly precipitation records were evaluated using the ratio of the root mean square error to the standard deviation of measured data(RSR).

where yj(m)and ?j(m)are the observed and estimated expected monthly precipitations at station m during the month j,respectively,?mis the observed mean and n is the number of missing values.

The RSR standardizes the root mean square error(RMSE)using the observed standard deviation.RSR varies from the optimal value of 0,which indicates zero RMSE or residual variation and therefore a perfect estimation,to a large positive value(Moriasi et al.2007).Percent bias(PBIAS)measures the average tendency of the estimated data to be larger or smaller than their observed counterparts(Moriasi et al.2007).On the contrary,Nash-Sutcliffe efficiency(NSE)is a normalized statistic that determines the relative magnitude of residual variance compared to measured data variance(Nash and Sutcliffe 1970).NSE indicates how well the plot of observed versus estimated data fits the 1:1 line(Moriasi et al.2007).

For testing the main and interactive effects of the radius of influence(e.g.number of neighboring stations)and estimation approaches,we applied a two-level factorial design considering the target station as the blocking factor(Quinn and Keough 2002),

where yijklis RSR calculated in the lthcross-validation iteration within the kthestimation approach within the jthradius of influence within the ithtarget station,Siis the target station(block),Rjis the radius of influence,Akis the estimation approach,(R×A)jkis the interaction between radius of influence and estimation approach and eijklis the error term.To confirm significant differences between factors(radius of influence or estimation approach)the Student–Newman–Keuls(SNK)test was used(Quinn and Keough 2002).A p-value of 0.05 was considered significant.

Results

Predictive capability of estimation approaches

The ANN and MLR approaches produced the best results for nearly all statistical criteria at both target stations 1 and 11,presenting a lower bias and higher precision compared to the other approaches(Table 2).On the contrary,the CCW approach showed the worst performance in terms of bias and precision for all target stations and radius of influence combinations.The variant IDWmproduced better results than IDW for all targetstationsand radiusofinfluence combinations,indicating that the inclusion of elevation differences improved the predictive capability.This result was somewhat expected given the existence ofa vertical precipitation gradient in this mountainous region.

Estimation approaches showed a decrease in RSR and PBIAS,as well as an increase in NSE,when they were applied to the higher elevation target station 1 compared to the lower target station 11(Table 2).In comparison to other approaches,IDW and CCW increase RSR and PBIAS and decrease NSE when the radius of influence increased from 25 to 50 km,that is,when the number of neighboring stations increased from 3 to 9.

Comparison of estimation approaches

The ANOVA showed significant differences(p＜0.0001)between estimation approaches(Table 3).Even though ANN and MLR have the lower RSR values(Table 1),the SNK multiple comparison test showed no significant differences with IDWm(Fig.3a).Additionally,IDWmhad a more significant difference than IDW and CCW(Fig.3a).This indicates that including elevation differences into the IDW significantly contributed to the improvement of its performance.The worst RSR values were obtainedwhen applying the CCW approach(Fig.3a)and its RSR values increased when the radius of influence was increased(Fig.3b).

Table 2 Predictive capability of the estimation approaches by target station and radius of influence(number of neighboring stations)

The ANOVA showed no significant effect of the radius of influence on RSR values,indicating that similar estimates of missing data can be obtained when considering 3 or 9 neighboring stations.However,as shown in Table 3 a significant interaction between the radius of influence(R)and the estimation approach(A)was detected(Table 3).This seemingly contradictory result is due to the opposite impact of the radius of influence on the method’s performance:RSR increases when increasing the radius of influence(Fig.3b)in IDW and CCW.In contrast,in the other approaches the RSR decreases when the radius of influence increases(Fig.3b).

Discussion

In this study,we compared five alternative approaches for estimating missing monthly precipitation records in two sectors in south-central Chile with complex terrain.The ANN and MLR showed higher precision and in most cases a lower bias compared to the other approaches.However,the precision(as per RSR)of IDWmwas not significantly different from ANN and MLR,according to the SNK test(p＜ 0.05).The ANOVAindicated that the radius of influence in terms of RSR did not significantly affect their predictive capability.However,this result can be explained by the significant interaction between the radius of influence(R)and the estimation approach(A).Therefore,an additional ANOVA was performed to evaluate the effects of the radius of influence on the predictive capability considering only the best three approaches:ANN,MLR and IDWm.For these approaches the radius of influence had a significant effect(p=0.036).Therefore,weconcludethatestimatesbased on nine neighboring stations located within a radius of 50 km are recommended for completing missing monthly precipitation datain theseregionswith complex topography.

Table 3 Analysis of variance for estimation approaches

Fig.3 SNK multiple comparisons for average RSR between estimation approaches(a)and interaction between radius of influence and estimation approach(b).In the upper panel,error bars represent the standard deviation and the letters(a,b,c)represent group methods that are not significantly different at α=0.05

Is there a “best”method?

Past studies have reported that the artificial neural network approach(ANN)was the best at estimating missing monthly precipitation records compared to otherapproaches(Teegavarapu and Chandramouli 2005;Khorsandi et al.2011).Coulibaly and Evora(2007)tested different neural networks architectures for completing daily precipitation records and found that the best method was the multilayer perceptron used in our study.In contrast,Alfaro and Pacheco(2000)in Costa Rica and Pizarro et al.(2009)in central Chile found that multiple linear regression(MLR)was the best method for filling in gaps in annual and monthly precipitation series,respectively.Thus,past research showed that ANN and MLR have emerged as robust methods for completing missing data in different geographical and climate settings(Kuligowski and Barros 1998).

Impact of elevation

The inclusion of elevation differences between the target and neighboring stations as a weight modifier to the IDWmsignificantly improved its performance.This is in agreement with studies that showed that including elevation differences in IDW had a positive impact on its predictive capability(Chang et al.2005;Golkhatmi et al.2012).Recently,Khosravi et al.(2015)used an altitude ratio(elevation of the target station divided by elevation of the neighboring station)to enhance the efficiency of the geographical coordinate method for completing gaps in annual precipitation series.

The Pearson’s correlation coefficients between the target and surrounding stations are presented in Fig.4.The values are moderately high(typically larger than 0.8)which is somewhat contradictory with the poor performance of the CCW method(e.g.,Fig.3a).We speculate that such a high correlation coefficient is due to the marked annual rainfall cycle,which is common among the seasons in this region and therefore this coefficient has little impact on the estimate of monthly precipitation values.Also there is a negative relationship between Pearson’s correlation coefficients and the elevation differences between the target station and its neighboring stations(Fig.4).This allows us to conclude that neighboring stations located at similar altitudes to the target station have a close relationship.

Impact of the radius of influence

Fig.4 Relationship between the Pearson’s correlation coefficients and absolute elevation differences between each target station and their neighboring stations:(a)target station 1(Andean foothills)and(b)target station 2(Central valley)

Even though the ANOVA showed that the radius of influence has a non-significant effect on precision(RSR),this factor interacted significantly with the evaluated approaches(Fig.3b).An increase of the radius of influence around the target station improved the predictive capability of only three of the evaluated approaches:ANN,MLR and IDWm.However,CCW and IDW showed a decreased performance when the radius of influence increased from 25 to 50 km,probably due to the association between decreased precipitations at the target and neighboring stations when distance from the target station increased(Johansson and Chen 2003;Mair and Fares 2011).Chen and Liu(2012)evaluated the IDW for interpolating rainfall data and found that the optimal radius of influence was in most cases up to 10–30 km.They also reported that the interpolation accuracy of this approach could become inferior when the number of considered rainfall stations exceeds the optimal value.

Conclusions

This study found that approaches based on artificial neuralnetworks(ANN),multiple linearregression(MLR)and IDWmhad the best performance in two sectors located in central-south Chile with a complex topography.Inclusion of elevation differences and Euclidian distances between targets and neighboring stations as weight modifier in the IDWmsignificantly improved overall estimates.Because the predictive capability of the three best approaches was significantly affected by the number of neighboring stations(radius of influence),we conclude that estimates based on nine neighboring stations located within a radius of 50 km are needed for completing missing monthly precipitation data.

Abbreviations

A×R:Interaction between estimation approach and radius of influence;A:Estimation approach;ANN:Artificial neural networks;ANOVA:Analysis of variance;CCW:Coefficient of correlation weighting;CV:Coefficient of variation;DGA:Dirección General de Aguas;ENOS:El Ni?o Southern Oscillation;IDW:Inverse distance weighting;IDWm:Modified inverse distance weighting;;m asl:Meters above sea level;MLP:Multilayer perceptron;MLR:Multiple linear regression;NSE:Nash-Sutcliffe efficiency;PBIAS:Percent bias;R:Radius of influence;RMSE:Root mean square error;RSR:Root mean square error to the standard deviation;SAS:Statistical analysis system;SNK:Student–Newman–Keuls test;UTM:Universal Transverse Mercator

Funding

This work was supported by the National Fund for Scientific and Technological Development(FONDECYT)[Project 1151050]and the first author gratefully acknowledges funding from Chile’s Education Ministry through the program MECESUP2[Project UCO0702].

Availability of data and materials

The data used in this study are available in public repositories of the Dirección General de Aguas(DGA;available at http://snia.dga.cl/BNAConsultas/reportes).

Authors’contributions

AB collected the data and performed the statistical analysis.AB and GT drafted the manuscript.GT and RG revised it critically for important intellectual content.AB,GT and RG gave final approval of the version to be published.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Escuela de Graduados,Facultad de Ciencias Forestales y Recursos Naturales,Universidad Austral de Chile,Box 567 Valdivia,Chile.2Departamento de Ciencias Forestales,Facultad de Ingeniería Forestal,Universidad del Tolima,Box 6299 Ibagué,Colombia.3Instituto de Bosques y Sociedad,Facultad de Ciencias Forestales y Recursos Naturales,Universidad Austral de Chile,Box 567 Valdivia,Chile.4Departamento de Geofísica,Facultad de Ciencias Físicas y Matemáticas,Universidad de Chile,Santiago,Chile.

Received:6 February 2018 Accepted:3 July 2018