The methodology employed is as follows.  (i) First, Regional frequency analysis (RFA) with the L-moments method. RFA combines data from nearby stations with similar characteristics (e.g., geographical proximity, elevation, and mean precipitation, etc.). L-moments are defined as certain linear combinations of order statistics of a random sample. In order statistics, random samples are ordered from the smallest to the largest. L-moments-based methods are less sensitive to outliers (unusually high) values in the data. (ii) Second, Kriging is applied to interpolate station data to the grid.  Kriging is an interpolation method considered better than conventional interpolation methods (e.g., distance weighting) because kriging accounts for proximity in data, spatial correlation, and also estimates interpolation uncertainty. (iii) Third, we quantify uncertainty in PF estimates. The methodology is shown to be applicable to two widely different regions, the Kissimmee-Southern Florida watershed (flat topography) and the Sacramento-San Joaquin watershed (complex topography). The PF estimates obtained from our methodology are shown to be statistically indistinguishable from the NOAA-Atlas 14 estimates. Moreover, our estimates are shown to be reasonably calibrated using probability integral transform (PIT)- a method popular in probabilistic forecasting used to estimate the accuracy of forecasts.