Acquisition of rainfall in ungauged basins: a study of rainfall distribution heterogeneity based on a new method

High-density rainfall data is always desired to capture the heterogeneity of precipitation to accurately describe the components of the hydrological cycle. However, equipping and maintaining a high-density rain gauge network involve high costs, and the existing rain gauges are often unable to meet the density requirements. The objective of this study is to provide a new method to analyze the spatiotemporal variability of the precipitation field and to solve the problem of insufficient site density. To this end, the proper orthogonal decomposition (POD) method is proposed, which can analyze the spatial distribution characteristics of rainfall fields to solve data shortages. To demonstrate the feasibility and advantages of the proposed methodology, four districts and counties (Hongshan District, Jianli County, Sui County, and Xuanen County) in Hubei province in China were selected as case studies. The principal results are as follows. (1) The proposed method is effective in analyzing the spatiotemporal variability of the rainfall field to reconstruct rainfall processes in ungauged basins. (2) Compared with the commonly used Thiessen polygon method, the inverse distance weighting method, and the Kriging method, POD is more accurate and convenient, and the root mean squared error is reduced from 3.22, 1.83, 2.19 to 2.09; the correlation coefficients are improved from 0.60, 0.85, 0.79 to 0.89, respectively. (3) The POD method performs particularly well in simulating the peak value and the peak time and can offer a meaningful reference for analyzing the spatial distribution of rainfall.


Introduction
Rainfall process is a highly heterogeneous process covering an extensive range of scales in time and space (Marani 2005;Nicótina et al. 2008). As the main input term of hydrological models, the variability of rainfall constitutes a significant source of uncertainty in hydrological simulations (Alexis and Guy 2004); therefore, accurate and reliable estimates of the location and intensity of precipitation are crucial for water resources management and hydrological studies (Van de Beek et al. 2010).
Applying rain gauges is one way to obtain accurate rainfall data (Sohn et al. 2010; Van de Beek et al. 2010); however, the limitation of this method is that it only provides point measurements and therefore lacks information on the spatial variability ( Van de Beek et al. 2010) unless used in a rain gauge network with sufficient density. Lebel et al. (1987) found that if rainfall measurements are based on ground measurements only, their accuracy depends on the spatial variability of the rainfall process and the density of the rain gauge network. Faurès et al. (1995) concluded that the use of a single rain gauge can lead to large uncertainties in runoff estimations for a small-scale (4.4 ha) catchment. This literature suggests that high-density rain gauge networks are needed to capture rainfall heterogeneity to accurately describe the components of the hydrologic cycle. However, equipping and maintaining high-density rain gauge networks involve high costs (Terink et al. 2018). In most areas, the sparse and nonuniform distribution of rain gauges does not meet the requirements of hydrologic simulations (Terink et al. 2018).
With the development of remote sensing technology, satellite rainfall monitoring has become a viable option to complement rain gauge observations. Several high-resolution satellite rainfall estimates are increasingly available at a quasi-global scale (Km et al. 2021). These products have the potential to address the temporal and spatial sampling limitations of rainfall stations. However, remotely sensed rainfall products are subject to systematic bias and random error (Kidd 2010), and their performance is highly affected by the topography and climatic conditions of the investigated regions (Romilly and Gebremichael 2010;Worqlul et al. 2014;Baez-Villanueva et al. 2018;Mekonnen et al. 2020). In addition, to represent the full spatial distribution of rainfall, massive satellite data are required; however, due to satellite orbit limitations, many satellite measurements are not continuous, and long records of satellite data are extremely limited (Segond et al. 2006).
The fusion of measured and satellite rainfall data is a new trend in quantitative rainfall analysis, and the consistent spatial and temporal resolution of multi-source data is a prerequisite for data fusion. In addition, different engineering projects have different demands for rainfall data with different spatial and temporal resolutions. Therefore, the interpolation method of precipitation, i.e., the spatial variation analysis of precipitation, has been widely studied (Dey and Mujumdar 2019). There are various spatial interpolation methods for rainfall data, such as the Thiessen polygon method, the inverse distance weighting method (IDW) method, and the Kriging method, which are commonly used (Jin and Heap 2011). The Thiessen polygon method proceeds as follows (Croley and Hartmann 1985): The interpolated points are connected two by two and the vertical bisectors of the connecting lines are made. The vertical bisectors intersect to form several polygons, thus dividing the large watershed into several subregions, each with one rainfall station, and the rainfall value in each region is represented by the measured value of the rainfall station within the region. The Thiessen polygon method assumes that rainfall varies abruptly across the boundary of the polygon and rainfall is uniformly distributed within the subregion, which is unphysical. The IDW method assumes that the interpolated points are more influenced by the closer stations than the far ones, and thus uses the reciprocal of the squared distance as the weighting factor (Chen et al. 2015). The disadvantage of this method is that the interpolation error is large when the data is significantly different from the neighboring points. The Kriging method (Armstrong and Matheron 1986) considers not only geographical features but also spatial correlation and variability of rainfall data and uses a covariance function to make the variables unbiased and optimal. However, when there are few rainfall stations, it is quite difficult to obtain the covariance function necessary for the Kriging method, which will also affect the interpolation effect to some extent (Gilbert and Simpson 1985). In addition, these methods are interpolated for a single time point, and the acquisition of rainfall processes for long time series requires iterative operations, which is very time-consuming.
To accurately and efficiently simulate the rainfall processes with high temporal resolution in large basins, it is crucial to fully consider the heterogeneity of rainfall fields. Many factors may contribute to the variability of precipitation, including tropical disturbances, El-Niño southern oscillation, Madden Julian oscillation, and so on, while the influence of spatial factors is the most significant (Pariyar et al. 2020). Variations in rainfall processes at different locations are often directly reflected in the temporal distribution of rainfall; therefore, they can be considered as "spatiotemporal variations". For example, owing to the different locations of sites B and A, rainfall at site B might occur later than at site A, and the rainfall pattern at site B might be gentler (steeper) owing to weakening (enhancement) of the air mass (Palynchuk and Guo 2011;Lin and Jhong 2015). Several scholars have tried to represent rainfall in continuous space and time mathematically. For instance, Northrop (1998) developed a spatial-temporal model based on the Bartlett-Lewis process for the temporal evolution of storm cells and the Neyman-Scott process for spatial arrivals. Wilks (1998) proposed an extension from a single station to multiple stations by driving singlestation models at various locations in New York with temporally independent but spatially correlated random numbers. Cowpertwait et al. (2002) proposed a spatial-temporal model based on the Neyman-Scott mechanism and tested the model on nine stations in Italy. Yang et al. (2005) used generalized linear models for the generation of multisite rainfall series in Southern England. Abas et al. (2014) presented a stochastic rainfall model for the generation of hourly rainfall data in an urban area in Malaysia. However, due to the influence of topography and climate, there is a large spatial gradient in rainfall intensity, such models currently assume spatial stationarity is not credible. (Alam and Elshorbagy 2015;Kwon et al. 2016). Improvements to include nonstationary properties of precipitation, to represent (1) the temporal characteristics of rainfall and (2) the spatial features of the data such as the location of high values of rainfall, topographic effects, or spatial continuity, are still under development .
Proper orthogonal decomposition (POD) is an effective dimensionality reduction method that uses basis functions extracted from a large number of known data to describe the original data in an optimal approximation (Aubry 1991). Because of its ease of interpretation and broad applicability to data from both simulations and experiments (Kutz et al. 2016), the POD technique is now used widely in fields such as fluid mechanics, aerodynamics, economics, statistics, and psychology to identify and analyze spatial distribution characteristics of some variables. For example, Meric et al. (2012) applied the POD method to extract information on the dynamic spatiotemporal evolution of global stock markets during the 2008 economic crisis, Cheveigné and Simon (2007) applied the POD method to neuroscience to extract the evolutionary structure of neurons in the human brain, and Varraso et al. (2012) used POD to extract and analyze the main spatial transmission of epidemic characteristics. However, this method has rarely been used in rainfall analysis. Considering that the structure of the vibration, flow, and aerodynamic fields is similar to that of the rainfall field, these fields also change 1 3 constantly with time and space. Besides, the velocity in the flow field and the pressure in the aerodynamic field are similar to the rainfall intensity in the precipitation field. Therefore, the POD method also has potential for applicability in the field of rainfall simulation, where it can be used to analyze the spatial distribution characteristics of rainfall fields and as an interpolation-like method to obtain rainfall data in ungauged basins.
This study is conducted with the following two objectives: (1) to develop a POD-based method to extract and analyze the spatial distribution characteristics of rainfall fields. This method can be used as an interpolation-like method to reconstruct rainfall processes in ungauged basins to obtain accurate simulations of rainfall processes with high temporal resolution; and (2) to test the applicability of the proposed method through comparative analysis with the interpolated results obtained using the Thiessen polygon method, the IDW method, and the Kriging method.

Methodology
The POD method, also known as the Karhunen-Loève procedure (Karhunen 1946;Loève 1955), is an effective dimensionality reduction method that can extract modes based on optimizing the mean square of the field variables being examined (Taira et al. 2017). It provides an objective algorithm to decompose a set of data into a minimal number of modes, each reflecting the influence of one or some physical factors on the field being decomposed, and these modes can be used to make the optimal approximate description of the original data. In this study, POD is implemented to decompose the precipitation field ( ∈ m×n ) into several orthogonal modes that reflect the influence of several physical factors on the rainfall field, and the modes reflecting the influence of spatiotemporal factors can be found through the subsequent comparative analysis. By reconstructing the rainfall field with the remaining modes after removing the modes reflecting spatiotemporal characteristics, the rainfall field structure independent of spatiotemporal factors can be obtained, and the rainfall processes in the ungauged basins can be reconstructed indirectly. The detailed steps and physical meanings of the variables involved in using the POD method to obtain the spatiotemporal structure of the rainfall field and to obtain data for ungauged basins are described below (Fig. 1).

Sampling and preprocessing of data
To study the influence of spatiotemporal factors on the rainfall field structure, two data sampling sessions are required. The first sample matrix X A consists of a certain rainfall process of sufficient duration at multiple stations in the study area, and the other sample matrix X B is manually spatiotemporally normalized based on the previous sample. By manual alignment (Roja et al. 2021), the rainfall processes at different stations start and end simultaneously in X B . With the above data collation, the original determinant X A has spatiotemporal variations, and the processed determinant X B can be considered to be uniformly distributed in time and space. X A and X B can be denoted as follows: where n is the length of the rainfall duration and m reflects the number of rain gauges in the study area. In fluid dynamics, the measured state u x, t n is related to the physics, i.e., the velocity, vorticity, and stream functions at each discretized spatial grid point. For precipitation data collected from a rainfall field, the rainfall intensity at different rain gauges at time t n can be denoted as u x, t n , which contains the mean of u (u(x) ) and the fluctuating rainfall intensity v x, t n . It is customary to subtract u(x) before computing the POD modes. The fluctuating rainfall intensity after removing u(x) at time t n can be denoted as follows: and the fluctuating rainfall field after removing u at time t n can be denoted as follows: where v x, t i stands for the ith fluctuating rainfall intensity at different stations in the study area.

POD calculation on the sample matrices
Various methods can be used to calculate POD, and the method selected for use in this study is the popular SVD (singular value decomposition)-based approach (Eckart and Young 1936). SVD can decompose X A and X B to get the modes and time factors. The mode represents the distribution structure of the rainfall parameter (intensity) in a certain physical sense, and the corresponding time factor reflects the trend of the mode over time. The purpose of applying POD in this study is to find the mode(s) that reflect spatiotemporal variations of the rainfall field, and to achieve this purpose, ̃ and ̃ are input to POD, respectively.
(2) v x, t n = u x, t n − u(x), In matrix form, the data matrix ̃ and ̃ can be decomposed with SVD as follows ( Fig. 2): are composed of the modes A and B . A reflects the spatial distribution information of rainfall intensity containing rainfall spatial variability, and B reflects the structure of the rainfall field after removing rainfall spatial variability. The n × n matrix = [a A1 , a A2 , ⋯ , a An ] and = [a B1 , a B2 , ⋯ , a Bn ] contain coefficients representing the temporal evolution of the modes and , respectively. The superscript T represents a conjugate transpose matrix operation (Higham et al. 2017). Matrix Σ holds singular values ( 1 , 2 , ⋯ , n ) along its diagonal. λ j ( j = 2 j ) is the energy contribution of the jth mode to the total variance, indicating the influence of the jth mode on the overall rainfall field. λ is listed in descending order, i.e., 1 ≥ 2 ≥ ... ≥ n ≥ 0 , the larger the value of λ, the greater the mode energy contribution. Furthermore, the percentage rainfall intensity energy contribution (E) of each POD mode can be obtained as follows:

Reconstruction of rainfall processes in ungauged basins
To find the modes that characterize the spatiotemporal factors of the rainfall field, the original data and the manually normalized data are subjected to POD in the previous section. It is assumed that the energy of the first r-order modes is close to the overall energy of the precipitation field, i.e., the modes reflecting the spatiotemporal influence factors of the rainfall field are included in the first r-order modes, and therefore, the first r-order modes are singled out for further comparison and analysis.
By comparing the mode values, numerical distributions, and temporal evolution coefficients of the predominant modes obtained from two POD calculations with different inputs, those modes with significant differences are excluded because they are considered to be related to spatiotemporal properties. Then, a homogeneous precipitation field is reconstructed with the remaining predominant modes. The fluctuating rainfall intensity v(x, t) can be expressed as a linear superposition of all orthogonal modes ( (x) ) multiplied by their corresponding time factors ( a(t)): where j is the subscript of the remaining modes after removing those modes related to spatiotemporal properties. Then, the precipitation field can be reconstructed as follows: Because the temporal and spatial variations of the precipitation field are eliminated, the errors between precipitations at close locations are greatly reduced, and the precipitation processes in the ungauged basin can be represented by precipitation processes at an adjacent station (Fig. 3).

Verification methods
To assess the performance of the methods used in this study, we used three performance statistics widely used to compare the model performance such as root mean square error (RMSE), relative error (RE), and correlation coefficient (Corr). Based on the results, performance of these methods is presented in detail.

Root mean squared error (RMSE)
The RMSE (Eq. (9)) estimates the standard deviation between the observed value and the simulated value. A small RMSE value indicates better performance, whereas a higher value indicates a poor performance.

Relative error (RE)
The RE represents the magnitude of the deviation of the simulated value from the true value and is calculated as follows:

Correlation coefficient (Corr)
The Corr is a statistical indicator to reflect the closeness of the correlation between variables and its value ranges between −1 and 1 representing a perfect negative and positive correlation, respectively. The correlation coefficient is calculated by the product-difference method and is based on the deviation of the two variables from their respective means. The formula for calculating correlation can be represented as: where Cov(X,Y) is the covariance of X and Y, Var[X] is the variance of X, and Var[Y] is the variance of Y.

Study area and data
Hubei Province is located in central China, between latitude 29° 01′ 53′′ -33°6′ 47′′ N and longitude 108°21′ 42′′-116° 07′ 50′′ E with a total area of 185900km 2 . Mountains, hills, and plain lakes account for 56%, 24%, and 20% of the total area of Hubei province, respectively. Except for the alpine climate in the high-mountain areas, most parts of the province belong to the subtropical monsoon humid climate, with sufficient light energy, abundant rainfall, and long frost-free periods. The average annual precipitation in Hubei province is between 800 and 1600 mm. The regional distribution of precipitation shows a decreasing trend from south to north (Fig. 4), with a maximum of 1400-1600 mm in the southwest and a minimum of 800-1000 mm in the northwest. The seasonality of precipitation distribution is obvious, with the most in summer (300-700 mm) and the least in winter (30-190 mm), respectively. The rainy season is from mid-June to mid-July, with the largest rainfall and intensity. To study the applicability of the method recommended in this paper, we selected several extreme rainstorms in Hubei Province for analysis. Hourly rainfall data at 56 stations in Hubei province are downloaded from European Centre for Medium-Range Weather Forecasts.
, 4 Results and discussion

Reconstruction of rainfall processes at ungauged basins
Because of the large spatially varying gradients of precipitation in Hongshan District, Jianli County, Sui County, and Xuanen County in Hubei province, these locations are assumed to be the ungauged basins for rainfall reconstruction (Table 1). Several extreme rainfall events in Hubei Province in 2020 were used as the original input (X A ) of POD to obtain the modes with spatiotemporal variability, and the manually normalized data (X B ) were input to obtain the modes with spatiotemporal variability removed. Figure 5 shows the energy contribution and the cumulative energy contribution of the modes obtained by decomposing X A and X B as POD inputs, respectively. The energy contribution of the modes derived from the decomposition of X A shows a greater descending trend, where the cumulative energy contribution of the first two modes is approximately 10% higher than that of the modes from the X B decomposition. This indicates that the X A -decomposed modes are more representative, which is consistent with the conjecture that the unprocessed data contain more spatiotemporal information. From the fifth mode onward, the energy contribution  of the X B -decomposed modes becomes larger than that of the X A -decomposed ones. And after the ninth mode, the cumulative energy contribution of the modes obtained from the POD of these two data sets is almost the same. In both cases, the cumulative energy contribution of the first six modes exceeds 80%, so these modes are selected as the predominant modes for subsequent analysis and reconstruction of the precipitation field. Figure 6 shows the predominant modes and their corresponding time factors before spatiotemporal normalization, i.e., raw data X A are used as input to the POD. As the energy contribution of the mode decreases, a higher proportion of null values is observed for modes with values close to zero, which is in line with the physical meaning of the energy contribution. In addition, as the mode's energy contribution decreases, the time factor corresponding to the mode fluctuates more frequently, which indicates that the mode becomes more unstable.
Because the energy contribution of the first mode exceeds 27.44%, it can be assumed that the reconstruction based on the first mode reflects the main structure of the precipitation field. The values of the first mode are mostly negative and its corresponding time factor is in the inverted triangular shape; therefore, a positive triangular rainfall intensity process can be obtained by multiplying the first mode and its corresponding time factor (Eq. (7)). The majority of the values of the second mode are positive, with the corresponding time factor appearing in the shape of a triangle. A similar rainfall process that increases firstly and then decreases can be obtained by multiplying the two. This can be interpreted as the structure of this precipitation field being a single-peak rainfall process. The values of the third to sixth modes are mostly around zero, and their corresponding time factors fluctuate frequently. This indicates that they play little role in the reconstructed rainfall field, but there may still be modes reflecting spatiotemporal factors in them.
To find the modes characterizing the variations of time and space, the processed data X B are used as input to POD to get the modes after spatiotemporal normalization. Similar to the modes obtained from the original data as input, the proportion of null values and the fluctuation frequency of the time factors increase with the decrease of the energy contribution of the modes. However, some variations are observed in the predominant modes and their corresponding time factors due to different inputs (Fig. 7). Compared with the first mode before spatiotemporal normalization, the values of the first mode after spatiotemporal normalization are also mostly negative, decreasing from northwest to northeast and slightly higher in the southeast. The main pattern of precipitation after spatiotemporal normalization is also a single-peak rainfall with a backward peak, which indicates that the main pattern of precipitation does not change significantly after the normalization of time and space. In addition, the distribution of the values of the fourth mode is almost the same before and after normalization; the values of the fifth mode show an increasing trend from west to east and reach the maximum at northeast, and the trends of time factors of the fifth mode before and after the normalization are the same. Therefore, it can be assumed that the first, fourth, and fifth modes do not reflect the influence of spatiotemporal factors on the structure of the precipitation field. However, the second, third, and sixth modes change significantly due to the normalization of time and space. Therefore, it can be assumed that the superimposed effect of the second, third, and sixth modes reflects the spatiotemporal variations of the precipitation field. Therefore, the reconstructed precipitation field after removing these three modes can be regarded as the result of removing the spatiotemporal variations. Figure 8 shows the reconstruction of precipitation fields before and after spatiotemporal normalization. The comparison shows that the rainfall processes at different stations are more consistent and the peaks are more pronounced after the removal of spatiotemporal variations. This confirms that modes 2, 3, and 6 contain information on the spatial and temporal distribution of rainfall, and modes 1, 4, and 5 better reflect the peak characteristics of rainfall.

Comparative analysis of different methods
Based on the results obtained in the previous section, the rainfall processes at the ungauged sites are reconstructed according to the method introduced in Sect. 2.3. To illustrate the  implementability and superiority of the proposed method, three most commonly used interpolation methods in current hydrological engineering: the Thiessen polygon method, the IDW method, and the Kriging method, were selected for comparative study. The rainfall data obtained by different methods as well as the measured rainfall data in the ungauged stations are shown in Fig. 9. It can be seen that a 24-h rainfall varies in magnitude and time distribution in different regions of Hubei Province, and the interpolation results of different methods applied in different regions also have great differences.
In comparison, the Thiessen polygon method performs poorly in the case study of this paper, especially in Jianli County where the rainfall gauge distribution is relatively sparse, because for district-and county-scale regions, the simplification of assuming an abrupt change in the spatial distribution of rainfall at a specific boundary is very rough. The accuracy of the interpolation results for rain intensity above 5 mm/h is low, with a Corr of only 0.6 and an RMSE of 3.22 mm/h (Fig. 10, Table 2). In addition, the interpolation results of peak time in Hongshan District and Sui County are late, which is unfavorable to the prevention and control of rain floods. Therefore, the Thiessen polygon method is only an alternate interpolation option for small areas where other methods cannot be employed. The interpolation performance of the Kriging and IDW methods is slightly better overall, but the performance of the peak value and the peak time is not satisfactory.
The peak value of the rainfall obtained by the IDW method shows a significant improvement compared to the Thiessen polygon method, with an RMSE of only 1.83 mm/h and a Corr of 0.85 (Fig. 10, Table 2). However, the interpolation of the IDW method for rainfall exceeding 15 m/h tends to be small and the performance in terms of peak time is not satisfactory. For example, there are large errors in the simulation of peak time in Jianli and Xuanen counties. This may be attributed to the fact that the IDW method mainly focuses on the magnitude of rainfall at the adjacent stations and does not give enough attention to the peak time. The interpolation results of the Kriging method are mostly a little smaller than the measured values, and the peak time is not accurate. This is unfavorable from the perspective of practical engineering safety. The Corr of the Kriging method is 0.79 and the RMSE is 2.19 mm/h, outperforming the Thiessen polygon method but worse than the IDW method (Fig. 10, Table 2). This large error of the Kriging method may be triggered by establishing the covariance matrix with insufficient dense sites in this study area.
The POD method exhibits a more accurate simulation of the peak value and peak time. The median and arithmetic mean of the relative errors of the POD simulations for the peak portion of the rainfall are minimal (Fig. 11), which is very meaningful for hydrological design. But perhaps because the extracted modes contain not only the spatiotemporal characteristics of the rainfall field but also part of the initial and final rainfall, the POD method is inaccurate in this case in simulating rainfall less than 5 mm/h at the beginning and end of the rainfall period, resulting in an RMSE of 2.09 mm/h. However, the POD method is more accurate for rainfall over 10 mm/h and the Corr exceeds 0.89 (Fig. 10, Table 2).

Conclusions
Spatial interpolation of rainfall has been regarded as the most common method to obtain rainfall data in ungauged basins, and the nonconsistency of rainfall spatial distribution is the source of error in the interpolation method. In this study, POD is introduced to analyze the spatiotemporal variations of rainfall fields and is used as an interpolation-like method to obtain rainfall data in ungauged basins. POD is a dimensionality reduction method that can be used to analyze the spatial distribution characteristics of different variable fields. Similar to its application to extract the flow structure in fluid mechanics (Redha et al. 2018), information on dynamic market transactions in finance (Meric et al. 2012), neuronal evolutionary structures in neuroscience (Cheveigné and Simon 2007), and virus transmission characteristics in epidemiology (Varraso et al. 2012), our study demonstrated that POD can efficiently identify and extract the structures of the rainfall field. By comparing the results obtained from two different sets of data as input, it is possible to find the modes that reflect the characteristics of the spatial distribution of the rainfall field, and the rainfall field reconstructed by removing these modes can be regarded as spatially uniformly distributed. Based on this, the rainfall processes in the ungauged stations can be reconstructed using the measured rainfall data from adjacent stations. This is particularly advantageous for urban areas with large spatially varying rainfall gradients.
In the application to the case study area in Hubei Province, POD is proved an effective and reliable approach for obtaining rainfall data in ungauged basins. Compared with the Thiessen polygon method, the IDW method (Goovaerts 2000), and the Kriging method (Haberlandt 2007), the first advantage of the proposed method is that it has a significant improvement in the simulation of the peak value and peak time. Second, the use of the POD method to remove spatiotemporal correlation modes to reconstruct rainfall in ungauged basins is once and for all, and it is more efficient to obtain long series of rainfall data with the POD method than other methods that require repeated interpolations. Despite the above advantages, this method has the potential for further improvement. For example, how to separate the parts of the same mode that may contain different characteristics of the rainfall field and how to quantitatively identify the modes that characterize spatial and temporal properties remain to be discussed. It is hoped that this study will serve as a stepping stone for the readers to become familiar with POD, analyze complex issues associated with spatial heterogeneity of rainfall fields, and further advance the developments of this Fig. 11 Relative error box plots of the peak portion of simulated rainfall obtained by different methods method in analyzing the effects of other factors (e.g., distance between stations and the ocean, slope, slope direction, wind speed) on the spatial distribution of rainfall data.