This vignette first provides an overview of the rainfall frequency estimates from the National Weather Service (NWS) followed by some examples on how to obtain and plot the data using the rainfreq package.
Rainfall frequency estimates for the USA from the NOAA National Weather Service's (NWS) division of Hydrometeorological Design Studies Center (HDSC) are often used in the the design of dams and other hydraulic structures and also in environmental planning and management. Data from NOAA NWS is available in various formats, including a user interface to extract the desired information. However, there is a lot of data and it is available in raw format as a large number of 1-km resolution GIS files.
The rainfreq package provides functionality to access 1-km rainfall frequency estimates in GIS format provided by the NWS' PF Data Server. The goal of the rainfreq package is to make the retrieval and analysis of this GIS data easier. Moreover, rainfreq also comes with NWS datasets on record point rainfall measurements.
After installing the package, load the package along with RCurl (for data extraction) and SDMTools, raster and maps for GIS analysis and graphics.
require(rainfreq)
require(RCurl)
require(SDMTools)
require(raster)
require(maps)
The main function provided by rainfreq is extract_freq. This could be used to extract data for any desired region. The default invocation of extract_freq of gets the 100-year 24-hour rainfall for the Southeast USA.
x_se <- extract_freq()
#> downloading data. might take a few moments ...
#> reading contents of the zipped file ...
print(x_se)
#> class : RasterLayer
#> dimensions : 1480, 1796, 2658080 (nrow, ncol, ncell)
#> resolution : 0.008333, 0.008333 (x, y)
#> extent : -94.92, -79.96, 24.46, 36.79 (xmin, xmax, ymin, ymax)
#> coord. ref. : NA
#> data source : in memory
#> names : layer
#> values : 6596, 16976 (min, max)
In order to obtain the 1000-year 48-hour rainfall for the Midwest, change region_name, storm_RP and storm_duration arguments accordingly.
x_mw <- extract_freq(region_name = "mw", storm_RP = 1000, storm_duration = "48h")
#> downloading data. might take a few moments ...
#> reading contents of the zipped file ...
print(x_mw)
#> class : RasterLayer
#> dimensions : 1934, 3239, 6264226 (nrow, ncol, ncell)
#> resolution : 0.008333, 0.008333 (x, y)
#> extent : -109.4, -82.38, 33.31, 49.42 (xmin, xmax, ymin, ymax)
#> coord. ref. : NA
#> data source : in memory
#> names : layer
#> values : 2804, 19478 (min, max)
Similarly, in order to obtain the 10-yr 6-hour rainfall for Hawaii, change the region_name, storm_RP and storm_duration arguments accordingly.
x_hi <- extract_freq(region_name = "hi", storm_RP = 10, storm_duration = "6h")
#> downloading data. might take a few moments ...
#> reading contents of the zipped file ...
print(x_hi)
#> class : RasterLayer
#> dimensions : 800, 1310, 1048000 (nrow, ncol, ncell)
#> resolution : 0.004167, 0.004167 (x, y)
#> extent : -160.3, -154.8, 18.91, 22.24 (xmin, xmax, ymin, ymax)
#> coord. ref. : NA
#> data source : in memory
#> names : layer
#> values : 1557, 11917 (min, max)
One could also obtain the record storm measurements provided by NWS using rainfreq.
data(rain_max_usa)
head(rain_max_usa)
#> Duration Amount_in Amount_mm Location Lat
#> 1 1-min 1.23 31 Unionville, Maryland 38.80
#> 2 5-min 2.03 52 Alamogordo Creek, New Mexico 34.66
#> 3 15-min 3.95 100 Galveston, Texas 29.29
#> 4 30-min 7.00 178 Cambridge, Ohio 40.00
#> 5 42-min 12.00 305 Holt, Missouri 39.45
#> 6 1-hr 13.80 351 Burnsville 6 WNW, West Virginia 38.88
#> Lon Start_Date Estimate
#> 1 -76.13 1956-07-04
#> 2 -104.39 1960-06-05
#> 3 -94.79 1871-06-04
#> 4 -81.58 1914-07-16
#> 5 -94.33 1947-06-22
#> 6 -80.77 1943-08-04 Yes
data(rain_max_world)
head(rain_max_world)
#> Duration Amount_in Amount_mm Location Lat Lon
#> 1 1-min 1.23 31 Unionville, Maryland, USA 38.80 -76.13
#> 2 3-min 1.75 44 Haughton Grove, Jamaica 18.33 -77.98
#> 3 5-min 2.48 63 Porto Bello, Panama 9.55 -79.65
#> 4 8-min 4.96 126 Fussen, Bavaria, Germany 47.87 12.17
#> 5 15-min 7.80 198 Plumb Point, Jamaica 17.93 -76.78
#> 6 20-min 8.10 206 Curtea-de-Arges, Romania 45.12 -24.42
#> Start_Date Estimate
#> 1 1956-07-04
#> 2 1925-09-30
#> 3 1911-11-29
#> 4 1920-05-25
#> 5 1916-05-12
#> 6 1889-07-07
The output from extract_freq is designed to be consistent with the RasterLayer class of the SDMTools package. This consistency enables the use of GIS functions for analysis and graphics provided by SDMTools and related packages.
Before plotting the data, convert the data to appropriate units. The original units are in 1000th inches, so multiply by 0.001 to get rainfall in inches.
x_se <- x_se * 0.001
x_mw <- x_mw * 0.001
x_hi <- x_hi * 0.001
Here are the three rainfall estimates obtained so far. State boundaries are added for spatial reference.
Rainfall amounts for selected frequency and duration periods - Souhteast USA.
# southeast
plot(x_se, breaks = c(6, 9, 12, 15, 18),
col = c("red", "yellow", "green", "blue"),
main = "100-year 24-hour Rainfall for Southeast USA (inches)")
map('state', region = c('florida', 'arkansas', 'louisiana', 'mississippi',
'alabama', 'georgia'), add = TRUE)
Rainfall amounts for selected frequency and duration periods - Midwest USA.
# midwest
plot(x_mw, breaks = c(2, 5, 10, 15, 20),
col = c("red", "yellow", "green", "blue"),
main = "1000-year 48-hour Rainfall for Midwest USA (inches)")
map('state', region = c('colorado', 'north dakota', 'south dakota', 'nebraska',
'oklahoma', 'minnesota', 'iowa', 'missouri',
'wisconsin', 'michigan'), add = TRUE)
Rainfall amounts for selected frequency and duration periods - Hawaii.
# hawaii
plot(x_hi, breaks = c(1, 3, 6, 9, 12),
col = c("red", "yellow", "green", "blue"),
main = "10-year 6-hour Rainfall for Hawaii (inches)")
extract_freq function's regional selection criterion does not include US territories such as Puerto Rico and Guam. Future updates would incorporate these regions.