The rainerosr package provides tools for calculating
rainfall intensity and erosivity indices from breakpoint rainfall data.
This vignette demonstrates the main functions and typical workflows.
I30 is the maximum 30-minute rainfall intensity, expressed in mm/hr. It is calculated by finding the 30-minute period within a storm that has the highest average rainfall intensity. This value is a key component in erosion prediction models like the Universal Soil Loss Equation (USLE) and its revisions.
EI30 is the rainfall erosivity index, calculated as the product of: - E: Total storm kinetic energy (MJ ha⁻¹) - I30: Maximum 30-minute intensity (mm hr⁻¹)
The result (MJ mm ha⁻¹ hr⁻¹) quantifies the erosive potential of a rainfall event.
library(rainerosr)Let’s create a sample storm event:
storm_data <- data.frame(
datetime = as.POSIXct(c(
"2024-01-15 14:00:00",
"2024-01-15 14:10:00",
"2024-01-15 14:20:00",
"2024-01-15 14:30:00",
"2024-01-15 14:40:00",
"2024-01-15 14:50:00",
"2024-01-15 15:00:00"
)),
rainfall_mm = c(2.5, 5.8, 8.2, 6.1, 3.4, 2.0, 1.5)
)
print(storm_data)
#> datetime rainfall_mm
#> 1 2024-01-15 14:00:00 2.5
#> 2 2024-01-15 14:10:00 5.8
#> 3 2024-01-15 14:20:00 8.2
#> 4 2024-01-15 14:30:00 6.1
#> 5 2024-01-15 14:40:00 3.4
#> 6 2024-01-15 14:50:00 2.0
#> 7 2024-01-15 15:00:00 1.5i30_result <- calculate_i30(storm_data, "datetime", "rainfall_mm")
print(i30_result)
#> $i30
#> [1] 49.2
#>
#> $total_rainfall
#> [1] 29.5
#>
#> $duration
#> [1] 70
#>
#> $start_time
#> [1] "2024-01-15 14:10:00 IST"
#>
#> $end_time
#> [1] "2024-01-15 14:20:00 IST"
#>
#> $interval_type
#> [1] "equal"The result shows: - i30: Maximum 30-minute intensity (mm/hr) - total_rainfall: Total storm depth (mm) - duration: Storm duration (minutes) - start_time and end_time: Period of maximum intensity
ei30_result <- calculate_ei30(storm_data, "datetime", "rainfall_mm")
print(ei30_result)
#> $ei30
#> [1] 347.37
#>
#> $total_energy
#> [1] 7.06
#>
#> $i30
#> [1] 49.2
#>
#> $total_rainfall
#> [1] 29.5
#>
#> $duration
#> [1] 70
#>
#> $ke_equation
#> [1] "brown_foster"The EI30 result includes: - ei30: Erosivity index - total_energy: Total kinetic energy - i30: Maximum 30-minute intensity - Plus other storm characteristics
The package supports three kinetic energy equations:
# Brown & Foster (default)
ei30_bf <- calculate_ei30(storm_data, "datetime", "rainfall_mm",
ke_equation = "brown_foster")
# Wischmeier & Smith
ei30_ws <- calculate_ei30(storm_data, "datetime", "rainfall_mm",
ke_equation = "wischmeier")
# McGregor & Mutchler
ei30_mm <- calculate_ei30(storm_data, "datetime", "rainfall_mm",
ke_equation = "mcgregor_mutchler")
# Compare results
comparison <- data.frame(
Equation = c("Brown & Foster", "Wischmeier", "McGregor & Mutchler"),
EI30 = c(ei30_bf$ei30, ei30_ws$ei30, ei30_mm$ei30),
Energy = c(ei30_bf$total_energy, ei30_ws$total_energy, ei30_mm$total_energy)
)
print(comparison)
#> Equation EI30 Energy
#> 1 Brown & Foster 347.37 7.06
#> 2 Wischmeier 358.85 7.29
#> 3 McGregor & Mutchler -4511.82 -91.70When you have data containing multiple storms, use
process_storm_events():
# Create data with two storms
multi_storm <- data.frame(
datetime = as.POSIXct(c(
# Storm 1
"2024-03-10 09:00:00", "2024-03-10 09:15:00",
"2024-03-10 09:30:00", "2024-03-10 09:45:00",
# Storm 2 (10 hours later)
"2024-03-10 19:00:00", "2024-03-10 19:20:00",
"2024-03-10 19:40:00", "2024-03-10 20:00:00"
)),
rainfall_mm = c(
4.5, 7.2, 9.1, 5.3, # Storm 1
6.8, 10.2, 8.5, 4.9 # Storm 2
)
)
# Process all events
events_summary <- process_storm_events(
multi_storm,
"datetime",
"rainfall_mm",
min_gap_hours = 6,
min_rainfall_mm = 10,
calculate_ei30 = TRUE
)
print(events_summary)
#> event_id start_time end_time duration_min
#> 1 1 2024-03-10 09:00:00 2024-03-10 09:45:00 60
#> 2 2 2024-03-10 19:00:00 2024-03-10 20:00:00 80
#> total_rainfall_mm i30 n_breakpoints ei30 total_energy
#> 1 26.1 36.4 4 223.48 6.14
#> 2 30.4 30.6 4 209.80 6.86Always validate your data before analysis:
validation <- validate_rainfall_data(storm_data, "datetime", "rainfall_mm")
if (validation$valid) {
cat("✓ Data validation passed\n")
} else {
cat("✗ Data validation failed:\n")
print(validation$issues)
}
#> ✓ Data validation passed
if (length(validation$warnings) > 0) {
cat("\nWarnings:\n")
print(validation$warnings)
}plot_rainfall_pattern(storm_data, "datetime", "rainfall_mm", plot_type = "both",
title = "Storm Event - January 15, 2024")
plot_intensity_profile(storm_data, "datetime", "rainfall_mm",
highlight_i30 = TRUE,
title = "Rainfall Intensity Profile")
If you have irregular time intervals or want to specify them explicitly:
data_with_intervals <- data.frame(
time = as.POSIXct(c(
"2024-01-01 10:00:00",
"2024-01-01 10:05:00",
"2024-01-01 10:15:00",
"2024-01-01 10:30:00"
)),
depth = c(3.2, 5.1, 4.8, 2.9),
interval = c(5, 5, 10, 15) # minutes
)
result <- calculate_i30(data_with_intervals, "time", "depth", "interval")
print(result)
#> $i30
#> [1] 61.2
#>
#> $total_rainfall
#> [1] 16
#>
#> $duration
#> [1] 35
#>
#> $start_time
#> [1] "2024-01-01 10:00:00 IST"
#>
#> $end_time
#> [1] "2024-01-01 10:05:00 IST"
#>
#> $interval_type
#> [1] "explicit"If you have your own event IDs:
data_with_event_id <- data.frame(
datetime = as.POSIXct(c(
"2024-01-01 10:00", "2024-01-01 10:15",
"2024-01-01 11:00", "2024-01-01 11:15"
)),
rainfall_mm = c(6, 8, 7, 9),
event_id = c(1, 1, 2, 2)
)
events <- process_storm_events(
data_with_event_id,
"datetime",
"rainfall_mm",
event_col = "event_id"
)
#> Warning in calculate_i30(event_data, time_col, depth_col, interval_col, : Data validation warnings:
#> Very few data points (< 3) may produce unreliable results
#> Warning in calculate_i30(event_data, time_col, depth_col, interval_col, : Data validation warnings:
#> Very few data points (< 3) may produce unreliable results
print(events)
#> event_id start_time end_time duration_min
#> 1 1 2024-01-01 10:00:00 2024-01-01 10:15:00 30
#> 2 2 2024-01-01 11:00:00 2024-01-01 11:15:00 30
#> total_rainfall_mm i30 n_breakpoints ei30 total_energy
#> 1 14 32 2 107.05 3.35
#> 2 16 36 2 142.88 3.97Here’s a complete workflow for analyzing rainfall data:
# 1. Load your data
my_data <- read.csv("rainfall_data.csv")
my_data$datetime <- as.POSIXct(my_data$datetime)
# 2. Validate the data
validation <- validate_rainfall_data(my_data, "datetime", "rainfall_mm")
if (!validation$valid) {
stop("Data validation failed")
}
# 3. Process multiple events
events <- process_storm_events(
my_data,
"datetime",
"rainfall_mm",
min_gap_hours = 6,
min_rainfall_mm = 12.7,
calculate_ei30 = TRUE,
ke_equation = "brown_foster"
)
# 4. Analyze results
summary(events$ei30)
mean_erosivity <- mean(events$ei30)
# 5. Visualize individual storms
for (evt_id in unique(events$event_id)) {
evt_data <- my_data[my_data$event_id == evt_id, ]
p <- plot_intensity_profile(evt_data, "datetime", "rainfall_mm")
print(p)
}
# 6. Export results
write.csv(events, "erosivity_results.csv", row.names = FALSE)Brown, L.C., & Foster, G.R. (1987). Storm erosivity using idealized intensity distributions. Transactions of the ASAE, 30(2), 379-386.
Wischmeier, W.H., & Smith, D.D. (1978). Predicting rainfall erosion losses: A guide to conservation planning. USDA Agriculture Handbook No. 537.
Renard, K.G., et al. (1997). Predicting soil erosion by water: A guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). USDA Agriculture Handbook No. 703.