Random forest techniques for spatial interpolation of evapotranspiration data from Brazilian’s Northeast

Computers and Electronics in Agriculture 166 (2019) 105017

Contents lists available at ScienceDirect

Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag

Random forest techniques for spatial interpolation of evapotranspiration data from Brazilian’s Northeast

T

José Clodoalves da Silva Júniora, Victor Medeirosa, Cicero Garrozia, Abelardo Montenegrob, ⁎ Glauco E. Gonçalvesa, a b

Department of Statistics and Informatics, Federal Rural University of Pernambuco, Recife, PE, Brazil Department of Rural Technology, Federal Rural University of Pernambuco, Recife, PE, Brazil

ARTICLE INFO

ABSTRACT

Keywords: Evapotranspiration Machine learning Spatial interpolation Agrometeorology

Irrigated agriculture is the human activity responsible for the highest consumption of water from the environment. The evapotranspiration represents, in practice, the consumption of water by a culture and the quantitative information of this parameter assists in a large number of water management problems. By employing spatial interpolation methods, one can determine evapotranspiration in places where there is no information of this parameter. This work evaluates algorithms for spatial interpolation of evapotranspiration data in terms of precision and performance. It compares conventional strategies as the Inverse Distance Weighting (IDW) and Ordinary Kriging (OK), and machine learning strategies, represented by the Random Forest (RF) and a Random Forest variation for spatial predictions (RFsp). The evaluation uses data of evapotranspiration from automatic meteorological stations located in the northeast region of Brazil, in January 2017. The leave-one-out crossvalidation method was used to compare the precision of each interpolation algorithm. The results showed that RF obtained better results for the estimation of the reference evapotranspiration than conventional approaches, coming to reduce to roughly half the error obtained. Despite to resolve blocky artifacts present in RF spatial distribution, RFsp did not present better results than RF, obtaining a very similar performance to IDW and OK. Concerning computational performance, IDW obtained the shortest time to create the interpolation model. The IDW also obtained the smallest prediction times.

1. Introduction

whose data is an essential source for many academic studies (Xavier et al., 2016; Bautista et al., 2009; Xavier and Brochado, 2017). INMET delivers at its website hourly data from more than 500 automatic weather stations all over Brazil (INMET, 2019). This data can be used to produce evapotranspiration maps at multiple time scales employing spatial interpolation techniques. Inverse Distance Weighting (IDW) and Kriging are two algorithms for spatial interpolation widely known in literature due to its low estimating error. However, machine learning algorithms can be considered a promising alternative to conventional algorithms (Li and Heap, 2014). Machine learning-based methods were used, for example, to interpolate data from seabed mud content samples in southwest Australia (Li et al., 2011) and for interpolating air temperature at Mount Kilimanjaro (Appelhans et al., 2015). In both cases, the machinelearning-based interpolation algorithms reduced error measures when compared to the conventional ones. Despite being tested to interpolate these variables, to the best of our knowledge, these methods were not employed to interpolate

Knowledge about the evapotranspiration (ET) rate plays a vital role in agriculture because this data is a central component in the hydrologic cycle. The evapotranspiration computed for a day measures, in practice, the water demands of culture and the amount of water to be replenished through irrigation. This way, reference evapotranspiration (ETo ) is the core variable to promote optimized water management in a smart irrigation system. This information is also essential to other water management process such as flood and droughts forecasting (Krajewski et al., 2017). Methods for ETo estimation based on meteorological data require access to environmental variables as air temperature, net solar radiation, relative air humidity, and wind speed. Hence, meteorology agencies (at the national or regional level) support public services for providing these environmental data in a digital format. In Brazil, one of the national agencies that maintain an environmental online data service is the INMET (Instituto Nacional de Meteorologia - National Institute of Meteorology),

⁎

Corresponding author. E-mail address: [email protected] (G.E. Gonçalves).

https://doi.org/10.1016/j.compag.2019.105017 Received 6 February 2019; Received in revised form 12 July 2019; Accepted 19 September 2019 0168-1699/ © 2019 Elsevier B.V. All rights reserved.

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

evapotranspiration data in Brazilian’s Northeast. Therefore, this paper presents an essay about the evaluation of machine learning-based algorithms against conventional strategies for interpolating evapotranspiration data from the Brazilian’s Northeast recorded in January 2017. Beyond showing the overall quality in terms of error based metrics of the different interpolation algorithms, this paper intends to show results on the computational performance of each algorithm in terms of the time spent to compute the model and its interpolation. The principal motivation to show those results was the importance of computational performance when considering to automatically execute interpolation as part of an information system supporting real-time decisions. In this case, the time to compute the model influences the response time of the decision. Depending on the system’s requirements, developers can choose to trade off quality for performance. The remainder of this paper is organized as follows: Section 2 describes the parameters, data, and metrics employed in our evaluation; Section 3 presents the achieved results, which are summarized and discussed at Section 4. Finally, Section 5 presents our main conclusions, suggestions and future work.

where N is the amount of points observed in the interpolated location’s neighborhood and F (xk) is the value measured at point xk . The weights W (xk) assigned to each point are given by

dk (x)

(2)

2.1.2. Ordinary Kriging Kriging adjusts a mathematical function to a specific number of points (or all points) within an established radius. From this function it is possible to determine the interpolated value for each desired point. Kriging considers that the direction and distance between the points of the sample reflect a spatial correlation that can be used to explain the variation of the measured variable surface. Several variations of the method can be found in the literature. In this work was applied Ordinary Kriging which is considered one of the most commonly used methods in spatial prediction (Li et al., 2011), in particular when spatial dependence is detected. The F (x) value interpolated at a point x by the Ordinary Kriging algorithm is given by

The main objective of our experiments is to measure performance in terms of quality of the generated models and the execution time of four different algorithms for interpolation of the reference evapotranspiration: Inverse Distance Weighting, Kriging, Random Forest and a Random Forest variation for spatial predictions. The first two represent the techniques conventionally used for spatial interpolation of environmental variables (Li and Heap, 2014), so the choice of IDW and Kriging is due to the fact that they are frequently used for interpolation of evapotranspiration (Xavier et al., 2016; Hodam et al., 2017; Malamos et al., 2017) and other environmental variables (Gerstmann et al., 2016; Shtiliyanova et al., 2017), producing small interpolation error. The last two algorithms evaluated in this study are from the machine learning area, and they have been used more recently as an alternative for the spatial interpolation of environmental variables. The Random Forest algorithm was selected due to its promising results obtained in recent studies, reaching low estimation error. In Li et al. (2011), this method was chosen from the existing machine learning techniques because of its superior performance over conventional methods to estimate other environmental variables, besides showing good combination with conventional methods. These hybrid models proved to be an excellent alternative to reduce the variability of the errors obtained. Commonly, the Random Forest is trained only with the coordinates of each spatial sample (latitude and longitude data), this way it tends to ignore in the modeling process the spatial correlation between the points, which can lead to biased predictions. Hengl et al. (2018) presents a recent proposal called Random Forest for spatial predictions (RFsp), that uses buffer distances of the observed points as explanatory variables, adding the effects of geographical proximity in the prediction process. This work also evaluates this variation. The next sections briefly describe all interpolation methods evaluated in this work.

N

F (x) =

kF

(x k )

(3)

k=1

where F (x k) is the value measured at the location, N is the number of measured values and k is the unknown weight for the value meaN sured at the k-th location, such that k = 1 k = 1. The process for generating the model includes two main steps: the adjustment of a function to the data of semivariance and the estimation of the weights k . Semivariance measures the degree of spatial dependence between two samples. The value of the semivariance between two points is directly proportional to the distance between them. Therefore, semivariance (h) is calculated using a set of pairs of observations of the form (F (x), F (x + h)) spaced by a separation vector h , using the following equation:

k th

(h) =

1 2 Nh

(F (x)

F (x + h))2

(F (x), F (x + h)) Nh

(4)

In Ordinary Kriging, the weight calculation is done in order to minimize the variance of the estimate, considering the constraint N = 1. Through the Lagrange multipliers method ( ), it is possible k=1 k to find a system of equations with n + 1 variables ( 1, 2 , …, n and ) and n + 1 equations by calculating the partial derivatives and equating them to zero, as shown in equation below: 1 C (x1,

x1) + x2) +

2 C (x2,

1 C (x1, 1 C (x1,

xn) +

2 C (x2, x n) +…+ =1

…

1

2.1.1. Inverse distance weighting IDW assumes that each measured point has influence (weight) on the prediction point, which will be inversely proportional to the distance to the forecasting location. Therefore, the algorithm gives greater weights to the points closer to the prediction point. Formally, the value F (x) interpolated at a point x is given by

+

2 +…+ n

2 C (x2,

x1) +…+ n C (xn, x1) + x2) +…+ n C (xn, x2) + n C (x n ,

xn) +

= C (x1, x) = C (x2, x) = C (xn, x) (5)

where C (x i, xj) is the semivariance calculated using the function adjusted in the previous step and considering the separation vector h between the points x i and xj . Given the weights, it is possible to estimate the value of the measure at point x i using the Eq. 3. 2.1.3. Random forest Random Forest is a supervised machine learning classifier based on the construction of a set of decision or regression trees. For making a tree, it is chosen as a random subset of variables that will be used to

N

W (x k ) F (x k )

p

where dk (x) is raised to the power value p and represents the distance between the interpolated coordinate x and one of the observed points x k , which can be calculated by Euclidean distance, Hamming distance, among others. In most works, including this paper, the Euclidean distance is used. Note that the weights in Eq. (2) are calculated such that N W (x k) = 1. i=1

2.1. Methods

k=1

p

N i=1

2. Materials and methods

F (x) =

dk (x)

W (x k) =

(1) 2

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

determine the outcome of the forecast. The number of trees (ntree ) and the number of variables available for selection in each split (mtry ) are two important parameters in the Random Forest training process (Houborg and McCabe, 2018). The Random Forest training algorithm (for both classification and regression) is as follows (Canty, 2002): 1. Draw ntree bootstrap samples from the original data. 2. For each of the bootstrap samples, grow an unpruned classification or regression tree, with the following modification: at each node of the created tree, rather than choosing the best split among all predictors, randomly sample mtry predictors and choose the best split from among those. Prediction of new data is made by aggregation of the predictions of the ntree trees. The final forecast will be the most common value returned (for classification problems), or the average value returned (for regression problems) by the decision trees that compose the forest. More details about Random Forest can be found in Breiman (2001). Although already used for spatial predictions, coordinate-based Random Forest ignores the spatial correlation between the points and the overall pattern of the samples during the modeling process, resulting in a non-spatial approach for spatial predictions. This characteristic may bias the results, leading to overestimated or underestimated predictions, especially where the spatial correlation in the estimated variable is high, and the sample point patterns show a clear trend. To solve this problem, Hengl et al. (2018) presents a framework that includes the effects of geographic proximity in the prediction process, where buffer distances from observation points are used as explanatory variables instead of the coordinates.

Fig. 1. Location of the automatic weather stations in the Brazilian’s northeast.

The INMET provides automatic weather station data updated hour by hour, so each analyzed station has up to 24 records in each day of the period used in the experiment. It was got from INMET the following environmental variables: hourly maximum and minimum air temperature, hourly maximum and minimum relative humidity of the air, hourly average wind speed at 2 meters, and hourly solar radiation. The time records were aggregated so that they could generate a single daily record using the daily maximum (Tmax ) and minimum (Tmin ) air temperature, the daily maximum (RHmax ) and minimum (RHmin ) air relative humidity, the average wind speed at 2 meters (u2 ) and the sum of daily solar radiation (Rs ). 2.3. Data quality

2.2. Study area and data

Before the calculation of the reference evapotranspiration, a validation was made for data quality control, in order to exclude stations that have discrepant data at each day. A station is removed from the data set for interpolation on a given day, if at least one of its meteorological variables (Tmax , Tmin, RHmax , RHmin, u2 and Rs ) shows values out its respective range as shown in Table 2. The value of Rs is validated according to the extraterrestrial radiation Ra . These rules were based on Xavier et al. (2016), which states these as the typical values for Brazil. Each day analyzed in the experiment period is considered as an independent basis for validation. This way, a station removed at a day can be considered at the other days, if its data is consistent. As can be seen in Table 2, no station was removed based on temperature, relative humidity, and wind speed variables. However, data from 112 records were removed due to solar radiation. While it seems a high value, Fig. 2 details the distribution of these removals at each day under study. As demonstrated, in the worst case, January 20, 117 stations were used in the evaluation.

The weather data for the calculus of the reference evapotranspiration were collected from the automatic weather stations of the INMET. These weather stations measure and record a set of environmental parameters (air temperature, air pressure, air humidity, rain, and others) according to a predefined schedule. These raw data are then sent to a database and made available for queries on the INMET’s site (INMET, 2019). The proposed experiment uses data from the automatic weather stations located in the Brazilian’s Northeast, an amount of 136 stations, with data from January 1 to 31, 2017. Within this period, only 6 of these stations had no available meteorological data, as shown in Table 1. Therefore these automatic weather stations were not included in the experiment, thus totaling 130 stations used. Fig. 1 shows the location of each weather station considered in this work. This Figure also shows the evapotranspiration values observed on January 1, 2017, at each station. The high values occur at the center of this region as expected, due to the semi-arid climate (B zone according to Koppen’s climate classification). The semi-arid climate is hinterland climate, typically found in Brazilian’s Northeast. This climate appears basically in landscapes where annual rainfall drops on average to less than 800 mm (Alvares et al., 2013). The high evapotranspiration values are explained due to the high evaporation rates and high average temperatures.

2.4. FAO 56 Penman-Monteith method For the calculus of the daily reference evapotranspiration (ETo ) at the meteorological stations it was used the FAO 56 Penman-Monteith method (Allen et al., 1998), as it is a standard way to estimate reference evapotranspiration from meteorological data. The Penman-Monteith

Table 1 Automatic weather stations that did not have data in the period of the experiment Weather station

State

Conde Estreito Itabaiana Nossa Senhora da Glória Senhor do Bonfim Uruçuí

Bahia Maranhão Sergipe Sergipe Bahia Piauí

Table 2 Validation rules applied to weather data.

Coordinates 11°48′S 06°39′S 10°67′S 10°20′S 10°26′S 07°28′S

Variable

37°35′W 47°25′W 37°47′W 37°43′W 40°08′W 44°20′W

Tmax , Tmin RHmax , RHmin Rs u2

3

Validation

0

30 Tmax , Tmin < 50 RHmax , RHmin < 100 0.03Ra Rs < R a 0

u2 < 100

No. of removals 0 0 112 0

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

Fig. 2. Stations removed from data at each day due to validation rules.

method uses standard environmental data that can be measured or derived from measured data and is defined by the following equation:

0.408*

Rn

ETo =

G +

Y * 900 * u2 * (es T + 273

+ Y (1 + 0.34*U2 )

Rnl =

=

(6)

(7) and the psychrometric parameter Y is calculated using the equation below (8)

in which Patm is the local atmospheric pressure (in kPa) that can be approximated by the elevation of the given weather station. The difference between es and ea is the saturation deficit and can be calculated using the expressions below:

ea =

es *RH 100

17.27*T T + 237.3

0.35

(11)

The data used to create the interpolation model of each algorithm at each day was formed by the longitude and latitude of the meteorological station (independent variables in the models), also obtained on INMET’s database, and the ETo observed, which will be the dependent variable in the generated models. It should be noted that the models generated in this work are only spatial models, not taking the time into account. This way, it is observed that the experiments at each day are independent of each other. For the IDW, p = 2 was used because this weighting parameter produced the best results in the revised works that deal with the use of IDW (Andes and Cox, 2017; Ly et al., 2011). In the interpolation of the chosen point, the five nearest meteorological stations were used to compose the interpolated point, as suggested by Xavier et al. (2016). For Ordinary Kriging, the exponential model was used to model the semivariance because this model offered the best results in our preliminary tests. In the Random Forest and its variation for spatial predictions, the number of trees was selected as 500 for all cases, as used in Genuer et al. (2010), and for the mtry was used the proposed default of p/3 for regression forests (Liaw and Wiener, 2002), where p is the number of predictors. The leave-one-out cross-validation method was used to compare the accuracy of each interpolation method used. For each day under analysis the following method was applied: one of the meteorological stations is removed from the training data set and, with the remaining data, spatial models are created based on each one of the algorithms used (IDW, Kriging, RF, and RFsp), then, for each algorithm, the evapotranspiration at the coordinates of the removed station is predicted and recorded. This

(T + 237.3)2

es = 0.6108*exp

Rs Rso

2.5. Quality of the interpolation models

17.27 * T T + 237.3

Y = 0.665*10 3*Patm

0.14 ea ) 1.35

where is the constant of Stefan-Boltzmann (MJ m 2 day 1); Tmaxk is the maximum absolute temperature during the 24-h period (in K); Tmink is the minimum absolute temperature during the 24-h period (K); and Rso is the clear-sky solar radiation (MJ m 2 day 1). The clear-sky solar radiation Rso is obtained by Rso = (as + bs ) * Ra , where as + bs is the fraction of extraterrestrial radiation reaching the Earth on days of clear-sky (n = N) and Ra is the extraterrestrial radiation constant (0.0820 MJ m 2 ). Considering these formula, the daily ETo value was calculated for each weather station. The Fig. 3 presents some spatial statistics of ETo for each analyzed day. It can be observed that there is a low spatial variability of ETo at each day in the studied region, with the coefficient of variation ranging between 0.2 and 0.35.

ea)

is the slope of the vapor pressure curve in relation to the where temperature (kPa °C 1); Rn is the daily net radiation (MJ m 2 day 1); G is the total daily flow of heat in the ground (MJ m 2 day 1); Y is the psychrometric coefficient (kPa °C 1); u2 is the wind speed at 2 meters high (ms 1); es is the vapor saturation pressure (kPa); ea is the current vapor pressure (kPa); T is the average air temperature (°C). While T and U2 are obtained directly from INMET (T is given by the mean between Tmax and Tmin ), and G 0 for daily estimates, the value of can be calculated by the expression

4098* 0.6108*exp

Tmaxk + Tmink (0.34 2

(9) (10)

in which RH represents the average relative humidity of the air (%), which can be estimated from the mean between RHmin and RHmax . The daily net radiation (Rn ) is defined by the difference Rn = Rns Rnl , in which Rns (net solar radiation) is employed using the equation Rns = (1 a)* Rs , where a represents the albedo for grass (typically a value of 0.23) and the solar radiation Rs (MJ/m2) that comes from the INMET data. The following formula defines the net longwave radiation Rnl . 4

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

Fig. 3. Spatial evapotranspiration statistics observed by day.

process repeats E times, where E is the number of weather stations with valid data at the day. The stations are removed circularly so that each station is removed from the training set only once. Following the described method, for each day under analysis and for each interpolation algorithm, the following statistical metrics are computed to evaluate methods: correlation coefficient (R), Bias, Mean Absolute Error (MAE), and Root Mean Squared Error (RMSE). Eqs. (12)–(15) define each one of those metrics, respectively:

corresponds to the estimation of the reference evapotranspiration using the latitudes and longitudes from the Brazilian’s northeast. 2.6. Algorithms’ computing performance Besides the metrics explained in Section 2.5, which indicates the quality of the interpolation models, some indicators of the computational performance of each algorithm are also measured. Such a study is essential to provide data that can drive the choice of a better algorithm for each case. For example, if the time to compute the model is a critical requirement, one can trade off the quality of the solution by the time to solve it. The time to generate the model is defined as the period spent to create the interpolation model. The performed activities for the elaboration of each model are dependent on the interpolation algorithm. For example, in OK, it is the time interval needed for the solution of the linear system, while in the Random Forest, it is the time interval for building the set of trees. The prediction time of the model was also measured, indicating the period spent to estimate evapotranspiration for a given coordinate. It depends on the internal operations made at each interpolation algorithm. For example, for IDW it equals the time interval taken to select the neighbors close to the given location and calculating the weighted mean, while in RF it is the time interval to go through the forest, from the root node of each tree to one of its leaf nodes, according to the criteria defined on each node. The values obtained for the generation time of the model and the prediction of the model along the cross-validation of each station were averaged. The evaluation was implemented in RStudio, using the R language. For the interpolation algorithms were used well-known libraries in the R language: for IDW and OK was used the gstat library,1 for the RF was used the randomForest library2 and for RFsp it was used the ranger library.3 These libraries were chosen because they were validated and widely used in the area (Li and Heap, 2014). The experiments were performed on a computer with the following configuration: Intel Core I3-3217U 1.8 GHz processor and 4 GB RAM.

E

EToi R=

ETo

EToi

ETo

i=1 E

ETo) 2

(EToi

(EToi

i=1

Bias = ETo MAE =

RMSE =

1 E

ETo)

2

(12) (13)

ETo E

EToi

EToi

(14)

i=1

E i=1

(EToi E

EToi )

2

(15)

In these equations, EToi and EToi are, respectively, the reference evapotranspiration value observated and the reference evapotranspiration interpolated on the day in a weather station i , EToi is the mean spatial reference evapotranspiration observed in one day and ETo is the mean interpolated reference evapotranspiration. The correlation coefficient R measures the linear dependence of the variables, ranging from −1 to 1. Thus, the closer R is to 1, more perfect is the positive correlation between the two variables, so if one increases, the other increases too. On the other hand, the closer R is to −1, the more perfect is the negative correlation between two variables. If R is 0, the two variables do not depend linearly on each other. Bias indicates whether the interpolated estimate tends to be smaller or larger than the observed data (0 is the ideal value). RMSE and MAE measure the accuracy, so that when the observed data and estimated data are similar, RMSE and MAE are close to 0, indicating a more accurate interpolation. RMSE calculates the square of the deviation between the estimated values and observed values, being more sensitive to significant errors. In order to provide a visual measure of the quality of each interpolation algorithm, it is also showed the spatial distribution of the interpolated reference evapotranspiration. The spatial distribution

1 2

12. 3

5

https://www.rdocumentation.org/packages/gstat/versions/1.1-5. https://www.rdocumentation.org/packages/randomForest/versions/4.6https://www.rdocumentation.org/packages/ranger/versions/0.11.2.

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

3. Results

interpolation of IDW and OK presents a smooth transition between the different ETo regions. The RF algorithm presented blocky artifacts which are due to its inherent characteristic. It works building different decision trees each one over a random subset of the training data. Each node at a decision tree corresponds to a decision question in the form “Is the latitude less than −8.03?” or “Is the longitude greater than −34.93?”, grouping data into blocks of similar values. In other words, RF focuses on data similarity instead of distance, which can indirectly capture aspects as altitude or vegetation cover that can change ETo . Each decision node in a tree of a random forest divides data according to a straight line parallel to longitude or latitude, depending on the decision question. Such characteristic is the principal cause of the blocky artifacts present in the RF spatial distribution. Increasing the number of trees in the forest can help to mitigate those artifacts – because each tree specializes in a subset of the data –, but even using 500 trees, as is the case, the spatial distribution of the random forest tends to present a blocky structure. Similar to the conventional algorithms, RFsp also presented smooth transition between the different ETo regions. The blocky artifacts presented by RF smooth and even vanish in RFsp. It occurs because it considers the distance between observation points as explanatory variables instead of the coordinates. This way, the divisions made at each node imposes lines at many different angles, which allows the RF to provide better borders, smoothing the surface.

This section presents and discusses the results achieved in the experiments. Section 3.1 presents the results about the accuracy of the interpolation models with each of the algorithms under analysis, and Section 3.2 shows the performance of the algorithms in terms of the time to generate the model and the time to interpolate evapotranspiration for a specific coordinate. 3.1. Accuracy analysis As a general result, it can be observed (from Figs. 4–8) that the Random Forest (RF) provided more accurate results than the other evaluated algorithms for all metrics, which is an impressive result since IDW and OK are well known in the literature by their low interpolation errors. Among the machine learning algorithms, RF was superior to the Random Forest for spatial predictions (RFsp) in all metrics too. Concerning the correlation between calculated and interpolated values (R), as shown in Fig. 4, a high positive correlation was observed for the RF, with values ranging from 0.93 to 0.96. RFsp obtained results much closer to IDW and OK algorithms, with slightly better results than these conventional algorithms at a few days. This result can be better visualized at Fig. 5 which shows four scatter plots considering data from all the valid stations at the first day of our dataset (January 1st, 2017). These scatter plots relate the 122 observed ETo with the ones estimated by each interpolation algorithm during the cross-validation. In these graphs, a perfect interpolation estimate should lie on the line with slope 1. As one can see the RF algorithm approximated better the line than RFsp, OK, and IDW. The RMSE, at Fig. 6, shows that RF provides the smallest error. RF achieved a much better result than the other algorithms with values between 0.5013 mm/day (at day 26) and 0.6982 mm/day (at day 2), which equates to an error of 8.48% and 11.64% relative to the ETo calculated at these two days. RFsp, in turn, presented a performance very similar to conventional algorithms in general, being only slightly better than IDW and OK in less than half of the analyzed days. The best results of RF can be confirmed by observing that, on average, the RMSE of RF is about half the RMSE obtained by IDW, OK, and RFsp at each day. The MAE results, shown in Fig. 7, compliments the RMSE results. An important aspect that should be highlighted is that the value obtained with this metric is slightly lower than those obtained with the RMSE, indicating that there is a small variability in the measured errors. In the case of RF, for example, a minimum MAE of 0.3990 mm/day (6.75% error on the average calculated ETo ) was reached at January 26th and a maximum of 0.5497 mm/day (9.04% error on the mean observed ETo ) at January 7th. RF achieved the best result among all analyzed algorithms. This superior result of RF can be verified again, noting that MAE for RF is again about half the MAE obtained by IDW, OK, and RFsp at each day analyzed. RFsp obtained very similar performance to IDW and OK in this metric also. From the Bias, shown in Fig. 8, one can see that RF and OK tend to underestimate evapotranspiration. In the case of the IDW and RFsp it is observed that they present an erratic behavior, sometimes underestimating, sometimes overestimating. In addition, the IDW has the highest bias on most days. The spatial distribution for the estimated ETo for each algorithm is shown at Fig. 9(a), (b), (c) and (d). The observations correspond to the interpolated ETo in the latitudes and longitudes that embrace the northeast region of Brazil on January 1, 2017. This day had 122 valid weather stations for estimation of reference evapotranspiration (the black dots represents these stations). The conventional algorithms for spatial interpolation obtained more significant variation in the estimates. For IDW, the ETo estimates ranged from 2.8181 to 10.7058, for OK estimates ranged from 3.2237 to 9.4588. For RF and RFsp algorithms the variations were from 3.4591 to 8.6243 and from 3.4602 to 9.1810, respectively. Moreover, the spatial

3.2. Computational performance In addition to the metrics that seek to indicate the quality of the models generated by the interpolation algorithms, the time to create the interpolation model (modeling time) and the time to predict ETo (prediction time) were also analyzed. Table 3 shows that the algorithm IDW obtained the shortest modeling time, ranging from 0 to 0.2 ms (the zero case occurs because the library used to measure time in R has its precision in milliseconds). This result occurs due to the relative simplicity of the IDW algorithm compared to the other algorithms analyzed. RFsp was the slowest algorithm with modeling time ranging between 5.05 to 7.08 s, which is due to the longer time necessary for the calculation of buffer distance of each observed point for all other points, this way for n points, n2/2 distances must be calculated. The modeling time for RFsp also is impacted by the construction of the set of decision trees used in this algorithm. Due to this reason, RF obtained the secondworst modeling time (from 51.1 to 61.5 ms). The inferior performance of RFsp can be proven by observing that, on average, the modeling time is more than 11000% greater than for RF. Regarding the prediction time for ETo , IDW reached the lowest value again, reaching a maximum of 1.7 ms and a minimum of 0.6 ms, as shown in Table 4. RF stands at the second lace, and it takes between 1.8 and 4.5 ms to predict ETo for one given coordinate. In the case of OK, the values ranged from 5.6 to 8.4. Finally, RFsp was the slowest algorithm with prediction time between 7.5 and 10.6 ms. The low coefficient of variation indicates that, although small, the mean values are representative of each interpolation algorithm and were little influenced by unexpected events generated by the operating system in which they were executed. 4. Discussion The results for RF based on coordinates confirm superior performance, in terms of accuracy, of the machine learning algorithms for spatial interpolation of other environmental variables seen in previous works (Li et al., 2011; Appelhans et al., 2015). The RF obtained more accurate results than other analyzed algorithms, reducing mean error metrics by approximately half. It was also observed that the estimates offered by the Random Forest tend to slightly underestimate the 6

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

Fig. 4. Correlation coefficient R for the algorithms under analysis. Fig. 5. Scatter plot of estimated ETo for 1 January 2017 using, respectively, (a) IDW, (b) OK, (c) RF, and (d) RFsp.

observed values of ETo . RFsp, in its turn, despite to resolve blocky artifacts present in the spatial distribution of coordinate-based RF, did not present better results than RF in terms of accuracy. In other words, comparing both RF and RFsp, there is an evident trade-off between a smooth spatial distribution and the measurement

accuracy. When compared to IDW and OK, which are known to provide excellent results for ETo spatial interpolation, RFsp obtained a very similar accuracy. Regarding the computational performance of the algorithms, the time to generate the model, was smaller than 70 ms for all algorithms, 7

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

Fig. 6. Root mean square error between observed and interpolated evapotranspiration for the four algorithms under analysis.

Fig. 7. Mean absolute error between observed and interpolated evapotranspiration for the four algorithms under analysis.

except for RFsp, which obtained the worst results. This result indicates that, for the number of stations used in this work, only RFsp, would not be able to be used for the daily estimation of ETo in an agrometeorological information system, for example. It should be noted that in such systems the interpolation model must be generated every day (or every hour, if more precision is needed) after the acquisition of the data from the previous day (or previous hour). Except for RFsp, the generation of the model for all other algorithms will be impacted much more by the acquisition of data, which involves database queries and data transformation, than by the generation of the model correctly. It should also be noted that the measurement of the execution time of an algorithm depends mainly on the size of the data input. In the experiment of this work, the data entry corresponds to the meteorological stations of the Northeast Region of Brazil which were about a hundred stations, but the increase in the number of stations will increase the time to generate the model. From the point of view of the time to predict, it is observed that RF could attend more than 310 interpolation requests per second, sequentially. This can be considered a high number, since the experiments were conducted in low-end hardware and an information system of the type that is glimpsed in this work would be executed in high-end

servers. Besides, since there are no dependencies between competing executions, the generated daily (or hourly) spatial model can be used in a distributed manner, being copied to multiple servers in order to satisfy interpolation requests in information systems. This work consisted of a prospective analysis about the use of a machine learning algorithm for the interpolation of ETo , which is, to the best of the authors’ knowledge, pioneering work on this topic. As a general result, this paper confirmed that the Random Forest algorithm is an attractive substitute for conventional interpolation algorithms widely used for estimating ETo . This result indicates that machine learning algorithms are viable alternatives to interpolate ETo and to compose agrometeorological information systems that can assist the farmer in his property at a low cost. It should also be noted that this viability is based on both precision and computational performance. 5. Conclusion This work evaluates, in terms of accuracy and computing performance, algorithms for the spatial interpolation of evapotranspiration data, comparing conventional strategies and machine learning strategies. The evaluation focused on the interpolation of evapotranspiration 8

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

Fig. 8. Bias between calculated and estimated ETo for the four algorithms under analysis.

Fig. 9. Spatial distribution of the reference evapotranspiration.

overcame, in terms of accuracy, the RFsp, a variation of the RF algorithm recently proposed to cope with spatial predictions. RFsp, in turn, achieved results very similar to IDW and OK. The algorithms were also evaluated about their computation performance in terms of execution time. The time obtained by all algorithms, but RFsp, is satisfactory for its usage in information systems to estimate daily (or hourly) evapotranspiration. The RFsp showed as an unfeasible option since the time to generate a model based on this algorithm increases at a fast pace with the number of spatial points in the model. The superior performance that the machine learning algorithms demonstrated for spatial interpolation of other environmental variables (e.g., percentage of sea mud (Li et al., 2011) and air temperature (Appelhans et al., 2015)) was also demonstrated in the estimation of the reference evapotranspiration in the experiment carried out. Due to the promising results achieved through this experiment in the selected period and region, a new experiment using data corresponding to the whole country in a more extended period will be carried out in the future. The performance of other machine learning algorithms will also be evaluated for comparison with the results of the algorithms already used in this work, as well as the performance of hybrid solutions combining conventionally used techniques and machine learning strategies (Li et al., 2011). Since such an approach merges the advantages of the two areas, this combination of machine learning techniques is expected to provide even more accurate results in the interpolation of evapotranspiration.

Table 3 Modeling time statistics (in milliseconds). Algorithm RF IDW OK RFsp

Minimum

Maximum

Mean

Median

Std. dev.

51.1 0 5.9 5053.1

61.5 0.2 8.9 7078.8

55.4 0.1 7.2 6107.1

54.9 0.1 6.9 6146.2

2.8 0 0.8 375.1

Table 4 Evapotranspiration prediction time statistics (in milliseconds). Algorithm RF IDW OK RFsp

Minimum

Maximum

Mean

Median

Std. dev.

1.8 0.6 5.6 7.5

4.5 1.7 8.4 10.6

3 1 6.8 8.8

3 1 6.7 8.9

0.7 0.2 0.7 0.6

data in the Brazilian’s Northeast, with data from the period of January 2017. The results showed that the coordinate-based Random Forest, a machine learning algorithm, obtained superior results for the estimation of ETo in relation to IDW and OK, which are algorithms widely used in the literature for the spatial interpolation of several different environmental variables, and are known for offering low error levels in the estimates of evapotranspiration. The coordinate-based RF also 9

Computers and Electronics in Agriculture 166 (2019) 105017

J.C. da Silva Júnior, et al.

Acknowledgements

Hodam, S., Sarkar, S., Marak, A.G., Bandyopadhyay, A., Bhadra, A., 2017. Spatial interpolation of reference evapotranspiration in india: comparison of idw and kriging methods. J. Inst. Eng. (India): Ser. A 98, 511–524. Houborg, R., McCabe, M.F., 2018. A hybrid training approach for leaf area index estimation via cubist and random forests machine-learning. ISPRS J. Photogram. Remote Sens. 135, 173–188. INMET, INMET’s automatic weather stations data webpage, , 2019. (Access in: 03/12/ 2019). Krajewski, W.F., Ceynar, D., Demir, I., Goska, R., Kruger, A., Langel, C., Mantilla, R., Niemeier, J., Quintero, F., Seo, B.-C., et al., 2017. Real-time flood forecasting and information system for the state of Iowa. Bull. Am. Meteorol. Soc. 98, 539–554. Li, J., Heap, A.D., 2014. Spatial interpolation methods applied in the environmental sciences: a review. Environ. Model. Software 53, 173–189. Liaw, A., Wiener, M., et al., 2002. Classification and regression by randomforest. R News 2, 18–22. Li, J., Heap, A.D., Potter, A., Daniell, J.J., 2011. Application of machine learning methods to spatial interpolation of environmental variables. Environ. Model. Software 26, 1647–1659. Ly, S., Charles, C., DEGRé, A., 2011. Geostatistical interpolation of daily rainfall at catchment scale: the use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrol. Earth Syst. Sci. 15, 2259–2274. Malamos, N., Tsirogiannis, I.L., Tegos, A., Efstratiadis, A., Koutsoyiannis, D., 2017. Spatial interpolation of potential evapotranspiration for precision irrigation purposes. In: Proceedings of the 10th World Congress on Water Resources and Environment, Athens, Greece, pp. 5–9. Shtiliyanova, A., Bellocchi, G., Borras, D., Eza, U., Martin, R., Carrére, P., 2017. Krigingbased approach to predict missing air temperature data. Comput. Electron. Agric. 142, 440–449. Xavier, F., Brochado, M.L.C., 2017. Use of spatial visualization for pattern discovery in evapotranspiration estimation. Proc. XVIII Brazilian Symp. Geoinform. 346–356. Xavier, A.C., King, C.W., Scanlon, B.R., 2016. Daily gridded meteorological variables in Brazil (1980–2013): daily gridded meteorological variables in Brazil (1980–2013). Int. J. Climatol. 36, 2644–2659.

The authors would like to thank CNPq (Grant No. 424979/2016-0) for its financial support. References Allen, R.G., Pereira, L.S., Raes, D., Smith, M., et al., 1998. Crop evapotranspirationguidelines for computing crop water requirements-fao irrigation and drainage paper 56. Fao, Rome 300, D05109. Alvares, C.A., Stape, J.L., Sentelhas, P.C., de Moraes, G., Leonardo, J., Sparovek, G., 2013. Köppen’s climate classification map for brazil. Meteorol. Z. 22, 711–728. Andes, L.C., Cox, A.L., 2017. Rectilinear inverse distance weighting methodology for bathymetric cross-section interpolation along the mississippi river. J. Hydrol. Eng. 22, 04017014. Appelhans, T., Mwangomo, E., Hardy, D.R., Hemp, A., Nauss, T., 2015. Evaluating machine learning approaches for the interpolation of monthly air temperature at Mt. Kilimanjaro, Tanzania. Spatial Stat. 14, 91–113. Bautista, F., Bautista, D., Delgado-Carranza, C., 2009. Calibration of the equations of hargreaves and thornthwaite to estimate the potential evapotranspiration in semiarid and subhumid tropical climates for regional applications. Atmósfera 22, 331–348. Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. Canty, A.J., 2002. Resampling methods in r: the boot package. The Newsletter of the R Project 2, 3. Genuer, R., Poggi, J.-M., Tuleau-Malot, C., 2010. Variable selection using random forests. Pattern Recogn. Lett. 31, 2225–2236. Gerstmann, H., Doktor, D., Gläßer, C., Möller, M., 2016. PHASE: A geostatistical model for the Kriging-based spatial prediction of crop phenology using public phenological and climatological observations. Comput. Electron. Agric. 127, 726–738. Hengl, T., Nussbaum, M., Wright, M.N., Heuvelink, G.B., Gräler, B., 2018. Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables. PeerJ 6, 1–49.

10