Easy-to-use spatial Random Forest-based downscaling-calibration method for 1 producing precipitation data with high resolution and high accuracy 2

. Precipitation data with high resolution and high accuracy is significantly important in 9 numerous hydrological applications. To enhance the spatial resolution and accuracy of satellite-based 10 precipitation products, an easy-to-use downscaling-calibration method based on spatial Random Forest 11 (SRF-DC) is proposed in this study, where the spatial correlation of precipitation measurements 12 between neighboring locations is considered. SRF-DC consists of two main stages. First, the 13 satellite-based precipitation is downscaled by SRF with the incorporation of high-resolution variables 14 including latitude, longitude, Normalized Difference Vegetation Index (NDVI), digital elevation model 15 (DEM), terrain slope, aspect, relief, and land surface temperatures. Then, the downscaled precipitation 16 is calibrated by SRF with rain gauge observations and the aforementioned high-resolution variables. 17 The monthly Integrated MultisatellitE Retrievals for Global Precipitation Measurement (IMERG) over 18 Sichuan province, China from 2015 to 2019 was processed using SRF-DC, and its results were 19 compared with those of classical methods including geographically weighted regression (GWR), 20 artificial neural network (ANN), random forest (RF), kriging interpolation only on gauge 21 measurements, bilinear interpolation-based downscaling and then SRF-based calibration (Bi-SRF), and 22 SRF-based downscaling and then geographical difference analysis (GDA)-based calibration 23 (SRF-GDA). Comparative analyses with respect to root mean square error (RMSE), mean absolute 24 error (MAE) and correlation coefficient (CC) demonstrate that: (1) precipitation importance of


SRF-based downscaling and then geographical difference analysis (GDA)-based calibration 23
(SRF-GDA). Comparative analyses with respect to root mean square error (RMSE), mean absolute 24 error (MAE) and correlation coefficient (CC) demonstrate that: (1) SRF-DC outperforms the classical 25 methods as well as the original IMERG; (2) the monthly-based SRF estimation is slightly more 26 accurate than the annual-based SRF fraction disaggregation method; (3) SRF-based downscaling and 27 calibration performs better than bilinear downscaling (Bi-SRF) and GDA-based calibration 28 (SRF-GDA); (4) kriging is more accurate than GWR and ANN, whereas its precipitation map loses 29 Topography could affect the regional atmospheric circulation and the spatial pattern of precipitation 165 through its thermal and dynamic forcing mechanisms (Jing et al., 2016;Jia et al., 2011). With the 166 increase of elevations, the relative humidity of the air masses increases through expansion and cooling 167 of the rising air masses, which brings precipitation (Jing et al., 2016). Thus, the precipitation-DEM 168 relationship has been widely employed to downscale satellite precipitation dataset. Here, the Shuttle with the spatial resolution of 90 m was downloaded from http://srtm.csi.cgiar.org/ and then resampled 171 to 1 km by the pixel averaging method. Since precipitation tends to be influenced by terrain variability 172 and terrain orientation, DEM derivatives including slope, aspect and terrain relief (Chen et al., 2020a) 173 were also used in the study. These derivatives were extracted from the SRTM DEM using ArcGIS 10.3. 174 The detailed information of all the datasets used in the study is shown in Table 1. 175

Methodology 177
The flowchart of SRF-DC is illustrated in Fig. 3, which includes three stages: data processing, 178 IMERG downscaling and downscaled IMERG calibration. It is noted that each IMERG pixel 179 represents the areal average precipitation within it, whereas rain gauge measurements are point-based. 180 Therefore, downscaling before calibration can decrease scale mismatch between pixel-based areal 181 precipitation and gauge-based point measurements.

Random Forest (RF) 185
RF is an ensemble of several tree predictors such that each tree relies on a random and independent 186 selection of some samples and features but with the same distribution (Breiman, 2001). The general 187 framework of RF is shown in Fig. 4. Specifically, each decision tree is constructed by randomly 188 collecting some training data with replacement, while the others are used to assess the tree performance 189 (sample bagging). When constructing each tree, only a random subset of features is selected at each 190 decision node (feature bagging). In the end, the majority vote for classification or the average 191 prediction of all trees for regression is used to obtain the final output. Overall, RF includes three 192 parameters to set: number of trees, depth of the tree, and number of features. In this study, the RF regression model was performed with the freely available codes, downloaded 203 from the website (https://code.google.com/archive/p/randomforest-matlab/downloads). 204

Spatial Random Forest (SRF) 205
In essence, the classical RF is a non-spatial statistical technique for spatial prediction, as it neglects 206 sampling locations and general sampling pattern (Hengl et al., 2018). This can potentially cause 207 sub-optimal estimations, especially when the spatial autocorrelation between dependent variables is 208 high. To this end, a spatial RF (SRF) is proposed in this study. The general formulation of SRF is as 209 follows: 210 where p(s 0 ) is the estimated precipitation at location s 0 , e is the fitting residual, f(▪) is the function 212 constructed by SRF, and X s and X ns are the spatial and non-spatial covariates, respectively. 213 In addition to spatial coordinates, one spatial covariate (X s ) is computed to account for the spatial 214 autocorrelation of precipitation measurements between neighboring locations: 215 where s i is the ith neighbor of s 0 , z(s i ) is the precipitation data of s i , w i is its weight, and n is the number 217 of neighbors. 218 In previous studies (Zhang et al., 2021;Li et al., 2017), the inverse distance weights (IDW) were 219 widely used. However, the IDW method only resorts to the spatial distance between the estimated 220 location and its neighbor locations, and does not consider the spatial autocorrelation between the 221 neighboring locations. To overcome this limitation, ordinary kriging (OK)-based variogram is adopted 222 to estimate the interpolation weights in this study by solving the following linear system: 223 where  is Lagrange parameter and     is the semivariogram.

225
It can be concluded that the variogram-based weights consider the spatial autocorrelation not only 226 between the known locations, but also between the known locations and the interpolated location 227 (Berndt and Haberlandt, 2018). In practice, the experimental semivariogram can be estimated from 228 sample data as follows (Goovaerts, 2000): 229 where n is the number of data pairs with the attribute z separated by distance h. 231 To obtain the semivariogram at any h, a theoretical semivariogram model should be fitted to the 232 experimental values. There are four commonly used theoretical semivariogram models: the spherical, 233 Gaussian, exponential, and power models. The best one with the highest fitting R 2 was used in the 234 study. 235

Working procedure of the proposed method 236
The detailed steps of SRF-DC are as follows (Fig. 3): 237 (1) Each pixel value of the 10 km IMERG was re-estimated by OK interpolation with its k nearest 238 neighbors (e.g. k=8) to obtain the interpolated IMERG (termed as 10km s I ), the 10 km IMERG 239 was interpolated by OK to obtain the interpolated 1 km IMERG ( where e is the fitting residual. 258 The 1 km precipitation data ( 1km C ) was produced based on the constructed relationship in step In this study, residual correction was ignored during downscaling and calibration, as many previous

Comparative methods 271
In the study, the performance of SRF-DC was comparatively assessed under three manners.
where n is the number of testing points, and E i and O i are the estimated and observed precipitation at 300 location i, respectively. 301  Table 2. Results show that all methods have poor 325 results for these observations. A possible reason is that high precipitation is often caused by 326 complicated environmental factors, which cannot be sufficiently explained by the constructed 327 predictors-precipitation relationships. In terms of ME, SRF-GDA ranks the first, which is followed by 328 kriging and SRF-DC. However, their ME values are less than -70 mm. With respect to RMSE and 329 MAE, kriging performs the best, which is closely followed by SRF-DC, while with respect to CC, 330 SRF-DC with the value of 0.64 outperforms the others. 331  all gauge stations. Overall, the RMSEs tend to be larger in the middle area, since it has higher 343 precipitation than the other areas (Fig. 1). BPNN (Fig. 7d) yields the poorest result, where many 344 stations have the RMSE values greater than 60 mm. It is followed by GWR (Fig. 7f). RF (Fig. 7c) and 345 kriging (Fig. 7e) are better than GWR and BPNN at most stations. SRF-DC (Fig. 7a) and SRFdis (Fig.  346 7b) are more accurate than the classical methods, especially at the stations in the middle area. Since the wettest month is July 2018 (Fig. 2), it is taken as an example to show the precipitation 351 maps of SRF-DC and some classical methods. Moreover, the semivariogram of kriging derived from 352 the original IMMERG and its prediction error map are shown, since they play an important role in the 353 performance of kriging and SRF-based methods. Results (Fig. 8) indicate that precipitation produced 354 by all the methods have spatial distribution patterns similar to IMERG, with much high precipitation in 355 the middle and low precipitation in the east. The ML-based methods have more spatial details of 356 precipitation than IMERG due to the inclusion of high-resolution predictors for precipitation estimation. 357

Results and analysis 302
The kriging map is so smooth that many details and variations of precipitation pattern are lost. This is 358 expected as it only uses ground measurements for the interpolation. RF shows obvious unnatural 359 discontinuities at the bottom. GWR suffers from systematic anomalies, with the values clearly greater 360 than their neighbors. In comparison, SRF-DC produces a good precipitation map. The semivariogram and prediction error map of OK are shown in Fig. 9. Obviously, OK has a 365 spherical model with the nugget variance (C 0 ) of 10.0 m 2 , sill (C 0 +C) of 10,560 m 2 , residual sum of 366 squares (Rss) of 8,800,611 m 2 , range (A 0 ) of 321,000 m, and fitting R 2 of 0.962, respectively (Fig. 9a). 367 The prediction error map (Fig. 9b) illustrates that the errors in the west are larger than in the east, and 368 in the boundary are larger than in the inner. It can be inferred that large errors are mainly located in the 369 areas with the sparse distribution of rain gauges. Moreover, the error magnitudes are not related to 370 RMSE distribution (Fig. 7) and precipitation pattern (Fig. 8). 2020b). Thus, they should be carefully selected to produce accurate precipitation data. 381

Model 382
In previous studies, the most commonly adopted model was GWR ( results than the original IMERG in the study region. RF and kriging outperformed GWR. Nevertheless, 386 the two methods have some shortcomings. For example, the precipitation map of kriging was so 387 smooth that many details were lost, and RF did not consider the spatial autocorrelation of precipitation 388 measurements. In comparison, SRF-based methods with the consideration of spatial autocorrelation 389 information demonstrated higher accuracy than the classical methods. Moreover, SRF-DC yielded 390 slightly better results than Bi-SRF, SRF-GDA and SRFdis. 391

Environmental predictors 392
NDVI, latitude, longitude and DEM-based parameters were commonly adopted as predictors to  Based on RF, the relative importance of each predictor (i.e. predictor importance estimate) is shown 400 in Fig. 10. Obviously, precipitation from kriging interpolation has the most importance. This is 401 because the interpolated value is directly related to precipitation. Kriging estimation is followed by 402 the downscaled precipitation. Longitude is the third most important variable, which is followed by The three LSTs also have a great impact on the precipitation estimation, where LST D seems slightly 408 more important than LST N and LST D-N . NDVI has a slight effect on the precipitation, which ranks last 409 but one. This might be due to the fact that NDVI is influenced by both precipitation and temperature 410 in the study site, and the low temperature above certain elevations hinders the vegetation growth. It is 411 less likely that the response of vegetation to precipitation has the delay, since SRF-DC on the monthly 412 scale is more accurate than SRFdis on the annual scale. 413 Among the 12 predictors, aspect has the least importance. This conclusion was also obtained by Ma

Temporal scale 418
Temporal scale has a great effect on the selection of predictors for precipitation estimation. There is a 419 debate on whether NDVI should be taken as a predictor for downscaling and calibration of monthly 420 precipitation. Some (Duan and Bastiaanssen, 2013;Immerzeel et al., 2009)  found that SRF-DC on the monthly scale was slightly more accurate than that on the annual scale (i.e. 427 SRFdis). This indicates that the response of vegetation to precipitation has no obvious time delay, and 428 NDVI can be used for monthly precipitation estimates. 429

Easy-to-use feature 430
Since the classical RF did not consider the spatial information in the modeling process, Hengl et al. layers-based RF, SRF is highly effective. Moreover, with the variogram-based kriging interpolation, the 441 spatial autocorrelation of precipitation not only between the gauge locations, but also between the 442 estimated location and gauge locations is taken into account. Thus, SRF has the merits of accuracy, 443 effectivity and ease of use. 444

Limitations and further researches 445
Although SRF-DC shows promising results than the classical methods, it still suffers from some 446 Berndt, C. and Haberlandt, U.: Spatial interpolation of climate variables in Northern 524 Germany-Influence of temporal resolution and network density, Journal of Hydrology-Regional 525 Studies