Improvements in rainfall estimation over Bangkok, Thailand by merging satellite, radar, and gauge rainfall datasets with the geostatistical method

ABSTRACT Bangkok is located in a low land area, and floods frequently occur from rainfall, river discharge, and tides. High-accuracy rainfall data are needed to achieve high-accuracy flood predictions from hydrological models. The main objective of this study is to establish a method that improves the accuracy of precipitation estimates by merging rainfall from three sources: an infrared channel from the Himawari-8 satellite, rain gauges, and ground-based radar observations. This study applied cloud classification and bias correction using rain gauges to discriminate these errors. The bias factors were interpolated using the ordinary kriging (OK) method to fill in the areas of estimated rainfall where no rain gauge was available. The results show that bias correction improved the accuracy of radar and Himawari-8 rainfall estimates before their combination. The merged algorithm was then adopted to produce hourly merged rainfall products (GSR). Compared to the initial estimated product, the GSR is significantly more accurate. The merging algorithm increases the spatial resolution and quality of rainfall estimates and is simple to use. Furthermore, these findings not only reveal the potential and limitations of the merged algorithm but also provide useful information for future retrieval algorithm enhancement.


Introduction
Bangkok is the capital of Thailand and the commercial center of the country. Urbanization has rapidly developed in Bangkok over the past 30 years (Edelman, 2022). Land cover has changed to a high density of tall buildings without vegetation cover and large areas of concrete, which are the most significant factors related to urban flooding after medium and heavy rainfall (Adulkongkaew et al., 2020;Visessri & Ekkawatpanit, 2020). Urban floods can cause transportation delays, power shortages, and significant economic losses, as well as spread disease (Wongsa et al., 2019). Generally, numerical weather prediction (NWP) and hydrological models are used to predict rainfall and floods in advance, respectively. However, to achieve high-precision predictions, the models require high spatiotemporal resolution for the precipitation input (Shen et al., 2018).
Rain gauges measure rainfall, but this method is unreliable due to the influence of wind, rain gauge density, wetting losses, and systemic errors (Ren & Li, 2007). At present, remote sensing techniques, such as ground-based weather radar or weather satellites, are broadly utilized due to the higher spatial and temporal resolutions of their estimates (Woldemeskel et al., 2013). Although ground-based weather radar can provide high spatial (1 km) and temporal (5 min) resolution, performing radar calibration is difficult (Abro, 2018;Villarini & Krajewski, 2010) because it is often limited by beam blockages due to complex terrain, overshooting beam altitude at key distances from the radar, or attenuation of the radar signal for radars with short wavelengths such as C-band and Xband radars (Rico-Ramirez & Cluckie, 2008;Friedrich et al., 2007;Giangrande & Ryzhkov, 2005).
Where a weather satellite can provide rainfall estimates on a global scale, satellite rainfall retrievals remain highly inaccurate due to their uncertainty compared to rain gauges (Maggioni et al., 2016). The physical constraints of the sensor, such as resolution, sampling, accuracy, precision, and noise, also impact the uncertainties in satellite precipitation observations. Furthermore, the unknown variance of the physical and statistical relationships in space and time can lead to a significant amount of uncertainty when algorithms are employed to convert the brightness of clouds or raindrops detected by satellites into precipitation (Ji et al., 2022;Stephens & Kummerow, 2007).
To obtain highly accurate rainfall data in both space and time, much effort has been expended into obtaining quantitative precipitation estimations (QPE) by merging radar, satellite, and gauge rainfall data (Sun et al., 2018;Nunez et al., 2008). For example, gaugeradar merged hourly rainfall products have a very high spatial and temporal resolution but have some gaps within the radar coverage area. The performance results have shown large differences in gauge rainfall and the gauge-radar merged rainfall products (Zhang et al., 2021;Kalinga & Gan, 2007Li & Shao, 2010).
To overcome the rainfall gap caused by combining point rain gauge measurements and Tropical Rainfall Measuring Mission (TRMM) gridded rainfall, a nonparametric kernel smoothing process was used (Nerini et al., 2015). The results showed an improvement in rainfall visual efficiency, especially in areas with low rain gauge density, but the accuracy is lower than the gauge-radar merged rainfall product (Chao et al., 2018;Chen et al., 2020;Grimes et al., 1999;Hu et al., 2015;Li & Shao, 2010;Ji et al., 2022;Woldemeskel et al., 2013). In Australia, Woldemeskel et al. (2013) have used thin plate smoothing splines (TPSS) and a modified inverse distance weight (MIDW) algorithm to merge daily rain gauge and TRMM 3B42 data. The study concluded that using interpolation techniques prior to integrating satellite rainfall estimates with a rain gauge can boost rainfall estimation, particularly in areas with few rain gauges (Hu et al., 2015;Pang et al., 2021;Yu et al., 2020;Zhang et al., 2021).
In the Ziwuhe Basin of China, researchers have used geographically weighted regression and wind speed to combine satellite and gauge rainfall resulting in improved correlation (R) and root mean square error (RMSE) values (Chao et al., 2018). Their merging process improved the spatial resolution and quality of the satellite product and was also simple to use. These results revealed that, for sparse gauge density, geographically weighted regression could skillfully estimate daily rainfall.
Furthermore, Wetchayont et al. (2013) have used a combination of radar, a geostationary satellite, and rain gauge data in northeastern Thailand. Cloud classification was performed to mitigate non-rainfall data before combining the radar and satellite data, which showed significant improvement (Wetchayont et al., 2013). Over northern and western China, a new gauge-radar-satellite combined rainfall dataset has been generated by applying a local gauge correction (LGC) and optimal interpolation (OI) merging strategies (Shen et al., 2018). The merged rainfall dataset was more consistent in all subregions and across all seasons; however, the cold season had more variability than the warm season. Furthermore, the combined rainfall dataset measured heavy rainfall more accurately than the original technique, demonstrating the benefit of full spatial coverage.
The city of Bangkok has numerous tall buildings that could block the weather radar signal. Additionally, because the radar is operated at a base scan with a 1° elevation angle, part of the horizontal data from the radar may be missed. This can cause the radar data to provide an underestimation of estimated rainfall. Satellite data is the best choice to fill in that missing data. However, the previously mentioned merged techniques were developed for merging gauge precipitation with radar or satellite precipitation, and extensive work has been performed all over the world on the combined gauge-radar-satellite precipitation measurements to enhance QPE quality temporally and spatially (Chao et al., 2018;He et al., 2018;Hu et al., 2015;Pang et al., 2021;Berndt et al., 2014;Sideris et al., 2014;Wetchayont et al., 2013;Woldemeskel et al., 2013;Yu et al., 2020;Zhang et al., 2021).
In particular, the recently launched next-generation geostationary meteorological satellites with high spatial and temporal resolution, such as the Himawari-8 satellite, have not yet been added to these research efforts. Therefore, this paper aims to improve rainfall estimates with high spatiotemporal resolution by merging gauge rainfall data with radar and Himawari-8 satellite data from Bangkok using a straightforward method. The findings not only reveal the potential and limitations of a merged algorithm but also provide useful information for future retrieval algorithm enhancement.

Study area and rain gauge network
Bangkok is located in the central region of Thailand, between 13.45° and 14.10° N and 100.30° and 101.00° E (Figure 1). The climatic characteristics of the region are related to the onset of the Northeast Monsoon in the winter season and the Southwest Monsoon in the summer season. The Department of Drainage and Sewerage (DDS), Bangkok Metropolitan Administration (BMA), has installed a network of rain gauges, including 130 automated tipping rain gauges. These weather gauges record rainfall at 5-minute intervals with 0.5-mm sensitivity and operate within areas with radar coverage. The monthly precipitation average varied from 14 to 323 mm between 1991 and 2020, according to the rain gauges, which cover 1,569 km 2 (Drainage and Sewerage Department, 2021).
The quality of the data from the rain gauges used in this study was examined based on their spatial and temporal consistency using the double mass curve approach. The inconsistent rain gauges were eliminated from the analysis. The rain gauge networks were used as a reference to measure the direct rainfall, were considered a "true" measurement, and were used to validate and correct the radar-and satellite-derived rainfall data. The gauge rainfalls accumulated from 5 to 60 minutes before the next step of the process. The 130 rain gauges were randomly split into two sets. Half of the gauges (65), set A, were appointed bias correction gauges (indicated by blue dots in Figure 2), and the other half of the gauges (65), set B, were used as reference gauges for comparison with the estimated rainfall in the cross-validation process.
Data from the rain gauges, radar, and Himawari-8 satellite for the storm events were collected and retrieved separately during the rainy season from August to October 2017. Furthermore, the events where the radar, gauge, and Himawari-8 data were not all available for the whole event were not included. Finally, ten of the precipitation events were selected for analysis (Table 1). Figure 1. Location of Nongchok radar (red star) and 130 automatic rain gauges in Bangkok, Thailand: A set (blue dots) for the bias correction process and B set (green dots) for the performance assessment process. The red circle shows radar range coverage (120 km).

Figure 2.
A merging technique was employed to combine the gauge, Himawari-8 satellite, and radar in this work.

Ground-based radar datasets
A ground-based C-band radar station is 231 meters above mean sea level in Bangkok, Thailand (13.83° N, 100.85° E). The DDS operates this station, which has 681 azimuths and one elevation (min = 1° and max = 3°). Radiation with a wavelength of 5 cm and a beam diameter of 1° is emitted, with a maximum range of 120 kilometers. This research was conducted based on radar reflectivity data with 5-min temporal resolution and 42-m horizontal spatial resolution from August to October 2017. The raw radar volume scans were processed, and quality control was used to flag errors, such as ground clutter, anomalous propagation, beam blocking, and attenuation due to intense rainfall conditions, prior to the calculation of QPEs as follows.
• First, removing clutter: a static clutter map created from the long-term reflectivity statistics was subtracted from the radar scan to eliminate clutter and other nonmeteorological echoes (Chang et al., 2009). • Second, filling the gaps: the radar dataset may contain gaps resulting from the clutter removal. Interpolation methods, such as the ordinary kriging (OK) approach, were employed to fill the gaps. In order to create hourly resolution, the original resolution data were transformed into a Cartesian grid with a spatial resolution of 1 km 2 horizontally and 1 km vertically (Chen et al., 2020;Grimes et al., 1999;Sinclair & Pegram, 2005). • Third, classifying the clouds: convective and stratiform clouds were extracted from the radar reflectivity using a method proposed by Steiner et al. (1995), and the threshold was set at 30 dBZ (Wetchayont et al., 2013). • Fourth, calculating the rainfall: the radar reflectivity (dBZ) was converted to radar reflectivity factor (Z) in mm 6 m −3 , as shown in Equation (1). Then, the radar reflectivity factor was converted to the rainfall rate (R) in mm h −1 using the Z-R relationship (Z = a R b ), in which a and b are coefficients that indicate different precipitation types, which are retrieved from the distribution of raindrop sizes (Sauvageot, 1994). In this study, Equation (2) (Burgess, 1992) and Equation (3) (Marshall & Palmer, 1948) were applied to the radar reflectivity for convective and stratiform clouds, respectively, as follows.
dBz ¼ 10 log 10 Z (1) where dBZ and Z are the radar reflectivity (decibel) and the radar reflectivity factor (mm 6 m −3 ), while Rr is the radar rainfall estimate (mm h −1 ).

Himawari-8: Advanced Himawari Imager datasets
Himawari-8 is a geostationary weather satellite operated by the Japan Meteorological Agency (JMA) since October 7, 2014. It is situated 35,793 km above the equator at 140.7° E.
The Advanced Himawari Imager instrument on the Himawari-8 satellite can capture data on 16 channels from visible to infrared, including three visible wavelength bands (red, green, and blue), to observe the atmosphere of the Earth. In this study, full disk images with a spatial resolution of 2 km × 2 km were collected every 10 min. We downloaded the Himawari-8 data with geo-correction accuracy Version 01 Yamamoto et al., 2020) at the same time as the ten selected rainfall events in Thailand from the Center for Environmental Remote Sensing (CEReS), the data server of Chiba University (ftp://hmwr829gr.cr.chiba-u.ac.jp/). Then the Himawari-8 data were averaged to hourly resolutions. Two different infrared brightness temperature (TBB) datasets, at 10.4 μm (TBB13) and 12.4 μm (TBB15), were used in this study. The geometric correction was used to correct cloud positional errors by using the top height of clouds provided by Hamada and Nishi (2010). Then, a split-window method (SWM) was applied to the Himawari-8 data to discriminate among different cloud types by comparing the water vapor absorption in two nearby infrared channels. The threshold strategy used by SWM is based on the relationship between TBB13 and its deviation (BTD) from TBB15 as a function of TBB13 (Inoue, 1987). Using a cloud classification threshold based on SWM for the summer season (Purbantoro et al., 2018), nine cloud types were identified over the Nongchok radar station: high cumulonimbus, middle cumulonimbus, cumulus, dense cirrus, ice cloud, water cloud, thick cirrus, cirrus, and thin cirrus. Then only the TBB13 data of the cumulonimbus clouds were used to estimate rainfall rates (Rs) using the non-linear regression empirical formula (Vicente et al., 1998) given in Equation (4). The estimated satellite rainfall was interpolated by the OK method into the same grid size as the radar data (1 km × 1 km).

Bias correction
The bias factor was obtained from the ratio of the true value to the corresponding estimated value. The estimated rainfalls were multiplied by the bias factors to eliminate bias values. The area-mean bias factor was retrieved from the ratio of the areal mean of the true rainfall to the areal mean of the estimated rainfall (Anagnostou Emmanouil et al., 1998). Thus, following the same explanation, the bias factor reaches a maximum value when the areal mean of both the true and the estimated rainfall is at maximum.
The assumption was that the true rainfall was the rain gauge observation while the estimated rainfall was the radar-and satellite-based measurements. Regarding the different spatial scales of the rain gauge, radar-based, and satellite measurements, a matchup technique (Chu et al., 2002) has been proposed as follows: while the rain gauge is a point measurement, rain estimates from the radar and satellite were obtained by averaging the values within a space window of 3 × 3 pixels of the radar-based and satellite data (9 km 2 ) with the center of the window on the rain gauge sites.
To remove rainfall estimation bias, a multiplicative correction was used to adjust the bias through a bias factor (B) (Anagnostou Emmanouil et al., 1998;Jongjin et al., 2016;Martinaitis et al., 2018;Nelson et al., 2015;Wang et al., 2012;Li et al., 2019), as shown in Equation (5). The B values can be obtained at the Rr and Rs pixels, where only rain gauges are available. Therefore, to derive the B value at the remaining pixels, the B values were interpolated using the OK method. By multiplying the B value with Rr and Rs at the corresponding position, the corrected bias of the rainfall estimations (Rr′ and Rs′) was achieved.
where Rg and Re are the rainfall from the gauge and estimate (radar or satellite rainfall), and i denotes the i-th gauge site.

Merging rain gauge, Himawari-8 satellite, and radar rainfall (GSR)
To merge rainfall data from different sources, the precipitating cloud regions were first collected from the radar and satellite raw data, and then the rainfall for each cloud type was measured. To create an optimized algorithm to combine the measured precipitation, an error relationship of multi-source precipitation measurements was applied. We used weighted averages of error variance to combine all the rainfall datasets. Finally, cross-validation was employed to evaluate the accuracy of the merged rainfall data ( Figure 2). Statistical objective analysis with a weighted average of the hourly rain gauge, satellite, and ground-based radar data over each pixel was applied to merge the rain gauge, radar, and Himawari-8 rainfall (hereafter, GSR) data, expressed in Equation (6) GSR ¼ Wr i Rr where Wr and Ws are the weights of the rain radar and the Himawari-8 estimate at the i-th grid, respectively, derived from Equation (7) and (8). Rr′ and Rs′ are the bias-corrected rainfall data from radar and Himawari-8 at the i-th grid as shown in Equation (9) .
The GSR is obtained from the weighted averages of two error variances (Grimes et al., 1999;Hu et al., 2015;Shen et al., 2018;Wetchayont et al., 2013;Woldemeskel et al., 2013). The spatial distribution of the B values highly influences the accuracy of GSR rainfall. The weights (Wri and Wsi) are computed based on the error variances of each rainfall estimate at a grid location, as follows: where ε r and ε s are the radar and Himawari-8 error variances at the i-th grid, which are estimated through the B value at the i-th grid, where ε e is the radar or satellite error variance. B ei is the bias factor at the i-th grid and B e is the averaged bias factor of the radar or satellite data.

Performance assessment
To assess the performance of the hourly GSR rainfall data, cross-validation was used.
The cross-validation is carried out by comparing the hourly rainfall data from the GSR and the rain gauges (set B). The performance scores of the GSR rainfall, which measure the degree of improvement over the gauge rainfall, were constructed (Ebert, 2007), and some error metrics, such as correlation coefficient (r), mean error (ME), mean absolute error (MAE), and RMSE (Abro et al., 2019(Abro et al., , 2020Zhu et al., 2021), were computed according to the following Equation (10)-(13), respectively.
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi P n i¼1 R gi À R gi RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 n where n is the total number of observations in the analysis and R ei and R gi are the estimated (radar, Himawari-8, and GSR rainfall) and observed rainfall at an individual rain gauge, respectively.

Cloud classification
The analysis of rain events was performed by selecting case studies that exhibit the average precipitation characteristics over Bangkok in the rainy season (from June to October 2017) (Promchote et al., 2016). On average, there are around ten rain events due to the influences of the summer monsoon, which can provide favorable dynamic and thermodynamic conditions for heavy precipitation (Day et al., 2015;Pillai & Annamalai, 2012). Using the data from the rain gauge network, episodes of rain events were selected. To be used for the case studies, at least ten rain gauges must have detected the rain and the radar image must contain a cluster of storms.
In the ten rain cases summarized in Table 1, the maximum and minimum rainfall averages were found on October 2, 2017 (4.02 mm h −1 ) and October 17, 2017 (2.01 mm h −1 ).

Radar cloud classification
Because of significant vertical convection, the horizontal structure of the reflectivity field was used as a marker for vertical motion characteristics. The stratiform regions are rather uniform, in contrast to the convective regions, which have an elevated reflectivity compared to their surroundings. This study used a 30-decibel (dBZ) threshold on radar reflectivity to distinguish between the convective and stratiform types (Wetchayont et al., 2013). Figure 3 is an example of the radar reflectivity (a, c, e, and g) and classified results (b, d, f, and h) in October 2017. On these days, radar detected a large area of precipitation caused by a mesoscale convective system (MCS). The results showed that the radar reflectivity and associated cloud type had agreed according to the reflectivity threshold utilized in this investigation. Additionally, this seasonal feature is observed in all precipitation structures dominated by the stratiform field.

Himawari-8 cloud classification
The results of cloud classification from the radar and Himawari-8 data were consistent with corresponding cloud types. Figures 3 and 5 show examples of radar reflectivity at 3 km, TBB13 (left panel), and the cloud classification results (right panel). Figure 4a shows the TBB13 of the rain event on October 2, 2017. The cloud classification threshold by Purbantoro et al. (2018) was applied to the Himawari-8 data; initially, there was no precipitating cloud (Figure 4b). Dense cirrus clouds dominated the results. Then, comparing the TBB13 and the gauge rainfall data for each cloud type revealed that the original threshold values in Purbantoro et al. (2018) developed for Himawari-8 with data from around Japan are inappropriate for use with cloud systems in Thailand (Figure 4e). Accordingly, the TBB13 threshold value for cumulonimbus cloud types was reduced to 200 K, and the BTD threshold value was reduced to 2 K. The threshold values of BT and BTD listed in Figure 4d were implemented by investigating a number of images (not limited to the ten datasets in Table 1). The original thresholds had been empirically adjusted and verified by object analysis and comparison with gauge rainfall data. The classified clouds from the SWM with adjusted threshold results in nine cloud categories, as shown in Figure 4c. Comparison between the TBB13 and the gauge rainfalls showed reasonable agreement with the corresponding cloud type ( Figure 4F). Almost all cumulonimbus pixels (rain clouds) were successfully identified. Cumulonimbus is the most common precipitating cloud in Thailand (Aumjira & Trivej, 2018). The cloud top data often exhibited heights over 15 km. Moreover, the precipitating clouds have TBB13 nearly below 200 K. Furthermore, the BTD of the cumulonimbus clouds (TBB13-TBB15) has a range of values, −4 to 1, Cumulonimbus which has been observed by multifunctional transport satellites (MTSAT) (Wetchayont et al., 2013).

Comparison between radar and Himawari-8 cloud classification
An example of cloud classification results can be seen for October 2, 2017; the radar (Figures 3(a-d)) and the Himawari-8 (Figures 5(a-d)) data show a squall line convective system that produced a large area of rainfall. Likewise, on October 3 and 23, 2017 (Figures 3(e-h) and 5(e-h), cumulus clouds, which often develop into cumulonimbus clouds due to the conventional process, formed over Bangkok. Figure 5 shows the distribution of classified clouds collected by the SWM for nine cloud categories at the same time as the radar data (Figure 3). The classification system identified almost all of the cumulonimbus cells (rain clouds), and the SWM with adjusted thresholds from this study is suitable for the cloud system over Bangkok.

Bias correction
Half of the rain gauges (set A) were chosen to adjust the rainfall bias on ten representative rainy days in 2017 (Table 1) for both radar and Himawari-8 rainfall data. The multiplicative bias correction must be computed separately for each grid cell. The OK interpolation technique was used to obtain the B values at each grid cell from 65 rain gauge sites. For radar and Himawari-8 estimates, the mean estimated rainfall ranged from 0.0012 to 3.80 mm h −1 and 0.02 to 2.53 mm h −1 , respectively, and the mean B values ranged from 0.0 to 17,916.6 mm h −1 and 0.0 to 725.0 mm h −1 , respectively. Figure 6 shows scatter plots of the rainfall measured by radar and rain gauges: (a) without bias correction and (b) with bias correction.
The numbers of matching pairs for each case in a set of radar-gauge pairs and Himawari-8-gauge pairs are shown in Table 1. The number of matched pairs differs due to differences in rainfall areas and sensing methods. The number of matchups is the largest on October 2 for radar (291) and , while the smallest number of matchups are 10 and 17 pairs on October 17 for radar and Himawari-8, respectively. The correlation coefficient (r), RMSE, and bias were used as metrics to assess the overall efficiency of the multiplicative bias correction process, as seen in each figure. The bias change improved the accuracy of rainfall forecasts substantially, shown by an increase in the correlation coefficient between hourly rainfall calculated from rain gauges and radar estimates, from 0.08 to 0.62 mm h −1 . After applying the bias correction, the RMSE, ME, and MAE all decreased. The scatter plots also show that the multiplicative process improves radar rainfall estimates when they are below 30 mm h −1 but underestimates high precipitation values. Figures 6c and d show a scatter plot of the rain gauge versus Himawari-8 rainfall values without and with bias correction. The multiplicative approach can also be used to increase Himawari-8 rainfall estimates. From 0.04 to 0.75 mm h −1 , the correlation coefficient Figure 6. Scatter plots of rain rates between the estimated versus the gauge, with and without bias adjustment of rainfall obtained from the radar (b and a), and rainfall obtained from the Himawari-8 (d and c), respectively. Data were derived from the chosen ten rainy cases. Blue, red, green, purple, yellow, orange, brown, black, grey, and maroon colors present rainy cases of numbers 1 to 10, respectively.
between the hourly rainfall measured by rain gauges and Himawari-8 estimates improved significantly. Unadjusted rainfall estimates had an RMSE of 4.35 mm h −1 , while the adjusted rainfall estimates had an RMSE of 2.95 mm h −1 . Additionally, the adjusted Himawari-8 rainfall tended to overestimate rainfall rates below ~20 mm h −1 and largely underestimated rainfall rates above that. However, some data showed no response to the multiplicative method in radar and Himawari-8 rainfall, according to those rain gauges located at the radar and Himawari-8 sites, as the rain estimates indicated them to be areas without rainfall. Figure 6b presents the significantly high performance of radar rainfall on October 2-3, 2017 (blue dots) due to high numbers of matching pairs compared to other cases. Meanwhile, Figure 6d reveals that there is no variation in the results in each case as long as the matching pair size is the same.
The results of the bias-corrected rainfall for the radar and Himawari-8 estimates for the sample rainy days, October 3 (upper panel) and 16 (lower panel), 2017, are shown in Figures 7 and 8, respectively. The bias corrections tended to significantly amplify the estimated rainfall of both radar and Himawari-8. The spatial distribution of the precipitation zones, as well as rainfall magnitudes, increased. The greatest rainfall amplification is found at rain gauges with high rainfall rates, except for the rain gauges located within a space window with zero estimated rain, as shown in Figures 7(e-f) and 8(e-f). Box areas in those figures represent rainfall location and quantity inconsistent between the radar and Himawari-8 rainfall fields and the gauge rainfall, which will continuously add uncertainty to the adjusted rainfall. Many data are missing from the radar reflectivity and the TBB data from Himawari-8 within the boxes, while the surface rainfall was observed by rain gauges.
Bias correction using a geostatistical method will not impact the corrected rainfall over that region. However, the adjusted rainfall for both radar and Himawari-8 estimates reveals that the bias correction with the geostatistical method is highly effective for low to medium rather than high rainfall rates. B is strongly linked to spatial variation in the Each column shows (a) radar reflectivity in decibels (dB); hourly rainfall distributions of (b) the unadjusted and (c) the adjusted rain of radar estimates. Dots represent rain gauge observations and the color scale indicates rainfall rate (mm h −1 ). rainfall field; therefore, the bias correction process can improve the magnitude and spatial distribution of the areas of precipitation.

Merging rain gauge, Himawari-8 satellite, and radar rainfall (GSR) results and performance assessment
The rain gauge, radar, and Himawari-8 rainfall data were merged using Equation (5), whereas the merging weights (Wri and Wsi) were calculated using Equation (6) and Equation (7). The weights are stationary over time but vary spatially from 0 to 1. The radar and the Himawari-8 rainfalls were calibrated from the rain gauge data; these adjusted rainfalls were used. Two data sources were merged to provide more accurate rainfall estimates. Figure 9 shows the rainfall fields on October 3, 2017 (07:00 UTC) and October 16, 2017 (15:00 UTC) based on ground-based radar (a and d), Himawari-8 (b and e), and GSR (c and f), respectively. The Himawari-8 rainfall fields have greater coverage than the radar fields, while the radar rainfall fields provide higher accuracy and a finer spatial scale.
Nevertheless, by merging these data into the GSR, the rainfall fields cover a greater area than before, with improved accuracy. As specified in the bias correction section, certain rainfall values from radar and Himawari-8 cannot be enhanced due to zero rain coverage within the space window. These areas with differences between radar and Himawari-8 rainfall are not in the same location. To fill the gaps not provided by radar or Himawari-8 rainfall, the merging method ( Figure 9) is utilized.
Cross-validation was performed using data from ten rainy events to determine the GSR rainfall estimation performance. Figure 10 shows a scatterplot of GSR versus gauge rainfall for sets A and B of rain gauges. The results are for ten rainy events over the study area at the rain gauge locations for each set. Figure 10a shows perfect results on GSR rainfall versus set A of the gauge rainfall with an r-value of 0.99, and a decrease in RMSE from 2.95 Each column shows (a) brightness temperature in K; hourly rainfall distributions of (b) the unadjusted and (c) the adjusted rain of Himawari-8 estimates. Dots represent the rain gauge location and the color scale indicates the rainfall rate (mm h −1 ). and 3.47 mm h −1 to 0.78 mm h −1 ; ME from −0.65 and −0.18 to −0.08; MAE from 0.74 and 0.75 to 0.16. Set A of rainfall gauges was used in the bias correction procedure; thus, the cross-validation performance after merging should improve.
In contrast, Figure 10b compares the GSR and gauge rainfall from set B, which was not included in any processes. The results show that the GSR rainfall maintains good performance but had more scatter, with an r-value of 0.90, RMSE of 2.03, ME of −0.12, and MAE of 0.41, although there was more dispersion due to underestimation. These underestimations were probably caused by the bias factors used in the bias correction and because the merged process was obtained from set A of rainfall gauges by the OK method. These results suggest  Figure 10. Scatter plot of rain rate between the estimated versus the gauge, the GSR against the a set rain gauges (a), and the GSR against the B set rain gauges (b), respectively, data obtained from the chosen ten rainy cases. that the GSR rainfall estimates will perform better when all the rainfall gauges in the network are used.

Discussion
This study formulated a framework that adopted a cloud classification and a bias correction combination approach to merge gauge, radar, and Himawari-8 rainfall data in the tropical climate. The cloud classification results showed that the Steiner et al. (1995) technique works well on the radar data relative to Wetchayont et al. (2013). The precipitating clouds in Bangkok were not well represented by the cloud classification based on the initial SWM threshold (Purbantoro et al., 2018), which was developed using data from Japan. The thresholds of TBB13 and BTD criteria are too high and low, respectively, for the Bangkok cloud system. Additionally, such criteria are modified according to the climate systems. Thus, we adjusted the SWM threshold to suit the cloud systems over Bangkok using the same method as Wetchayont et al. (2013). Then, we obtained rainfall estimates using appropriate equations for each cloud type. The SWM threshold values are slightly different from those of previous works (Purbantoro et al., 2018) as the BTD has no negative values; this may vary because of the difference in the climate system.
The comparison between the estimated (from radar and Himawari-8) and the measurements (from rain gauges) showed only small improvements for higher rain rates than using only one equation for estimates from radar and satellite data. This may be due to the limitations in radar and satellite measurements, as previously discussed. Missing data is common; the current system misses many cloud features observable on radar  because the Nongchok radar station has only one elevation angle operation, causing the radar assessment of rainfall to be underestimated. As shown in Figure 6b, the matched pairs size effect had higher performance on October 2-3, 2017, due to a larger number of matching pairs, which resulted in higher adjustment performance. Thus, a small dataset size is insufficient for adjusting the rainfall among different rainy cases. The radar should be observed from various elevation angles. Overall, the radar and Himawari-8 rainfall tend to underestimate the rainfall.
The rainfall fields of estimation and the rain gauges have different geographical variability. This variation in position is primarily caused by wind drift, changes in raindrop size due to evaporation (Dai et al., 2019), lack of detection in local deep convection by remote sensing systems (Li & Shao, 2010), or missed detection of small-scale precipitation systems due to low rain gauge density (Dai et al., 2013;Ren & Li, 2007). An increased number of rain gauges within the network is needed to further overcome this problem.
However, bias correction was used in this study to mitigate these errors. The bias correction performs well for all rainfall rates except for an area with no estimated rainfall within the space window and its surrounds. The adjusted rainfall shows a greater improvement in the radar rainfall than the Himawari-8 rainfall. These findings agree with other studies of Thailand (e.g. Wetchayont et al. (2013)).
By merging with the bias factor obtained from the geostatistical method, the GSR rainfall dataset improves not only the magnitude but also the spatial distribution. The bias factor is closely related to spatial variation in the rainfall field. Usually, the estimated rainfall bias correction is limited by the low density of traditional rain gauge networks. However, the bias correction combined with the geostatistical method has the advantage of the low-density rain gauge network and high spatiotemporal resolution compared to the common bias correction method. This agrees with the general conclusions drawn by Wetchayont et al. (2013), Chao et al. (2018), and Chen et al. (2020).
The performance of the GSR rainfall product for overall events showed greater reliability with smaller numbers of missing data. An overestimation (underestimation) of low (high) rainfall rates was found in this study and is consistent with past literature (Chen et al., 2020;Hasan et al., 2016;Shen et al., 2018;Steiner & Smith, 2002;Wetchayont et al., 2013). In this regard, the insufficient reproduction of the extremes of stronger rainfall rates at a single point of the GSR rainfall product is due to the upper limit of the sensitivity of the radar and satellite sensors to rainfall remote sensing.
Overall, the cross-validation results performed in good agreement with the rain gauges and reduced the statistical error indicators for the GSR rainfall product in Bangkok ( Figure 10). The findings of this study are an essential data source for hydrological models to ultimately improve runoff simulations, even though the rain gauge data used for bias correction and validation do not have the same quantities and geographical locations. Furthermore, the merging procedure can be reproduced in other areas, especially in regions with low rain gauge network density, complex topography, or data loss in radar scan patterns.

Conclusions
This research presented a scheme for combining gauge, radar, and satellite rainfall observations, as well as using cloud classification and bias correction as additional factors to reduce the uncertainties of rainfall estimation. The cloud classification is adopted to radar reflectivity and TBB data from the Himawari-8 satellite to retrieve the rain rates for each cloud type. Moreover, the bias correction has been applied to both the radar and Himawari-8 rainfall data. With the surface rainfall measurements from a rain gauge network as a reference, the bias correction with the geostatistical method can significantly improve the correlation between radar and Himawari-8 data with the rain gauge dataset and reduce the difference for most events. The adjusted rainfall improved the gauge-radar correlation from 0.08 to 0.62 and reduced the RMSE from 4.41 mm h −1 to 3.47 mm h −1 for ten selected events. Meanwhile, the adjusted rainfall improved the Himawari-8-gauge correlation from 0.04 to 0.75 and reduced the RMSE from 4.39 mm h −1 to 2.95 mm h −1 . These results revealed that the bias correction using rain gauges as a reference is an important process of all estimated rainfall quality control.
The merged gauge-radar-satellite scheme showed better performance and is recommended for future high-accuracy rainfall estimation. An improvement of two-thirds in the root mean squared error was obtained using a linear weighted combination method compared to the original radar and satellite rainfall rates. Almost all gauge locations where the GSR rainfall was evaluated significantly improved their rainfall estimations. The interpolation of the bias factor as an additional process in the combination algorithm further improved the precipitation estimates compared to the algorithm without interpolation of the bias factor. The consistency between the surface and the GSR rainfall demonstrates the effectiveness of this method.
Some problems remain with the proposed scheme, such as the limitations of the sensitivity of the radar and satellite sensors and different rain field locations. However, this study is a preliminary attempt at merging the estimated rainfall from different sources, which can be easily applied to other study areas. In future work, the long-term continuous data from various climate and geographical conditions will be explored to improve the performance of the proposed method. Further, the GSR rainfall has high spatiotemporal rainfall data on the catchment scale, and the resulting data is available in remote areas because of extended radar and satellite data coverage. Using the GSR data with hydrological models to simulate river flow and model extreme events would be an interesting next step for this study.