This code builds two conditional extreme-flow dependence models (Scenario A/B), selects copulas by AIC, generates joint extreme simulations, derives joint return-period contours, merges them into a conservative envelope, and then ranks where along that envelope joint outcomes are most probable using KDE—yielding both curves and a representative “most likely” design point.
Here is a summary of the steps:
Download daily streamflow from two gages on the mainstem of the Iowa River
Extract a partial duration series at each gage
Compute rank correlation coefficient (Kendall’s tau) between the two locations
Define copula scenarios
Fit copula models and select copula for each scenario based on Akaike Information Criterion (AIC)
Generate samples from each copula and transform to original variable scale
Derive isoline for specified probabilities (return periods)
Combine isolines into envelope curve
Compute relative density of pairs on each isoline using kernel density estimation
Functions
Define custom functions. There are additional functions in later code blocks. I’m sorry for my R sins.
# Computes a time separation criterion, in days, between events in a PDS to ensure independence#Reference: Bulletin 17# DA: drainage area in square miles# returns days between eventsB17_days_between <-function(DA) {ceiling(5+log(DA))}# Madsen annualization of Generalized Pareto distribution to Generalized Extreme Value distribution # gpd_par: Generalized Pareto distribution parameters# lambda: rate of events# returns Generalized Extreme Value (GEV) distribution parameters# Reference: Madsen, H., Rasmussen, P., & Rosbjerg, D. (1997). Comparison of annual maximum series and partial duration series methods for modeling extreme hydrologic events 1. At-site modeling. Water Resources Research, # 33(4), 747-757.madsen_annualization <-function(gpd_par, lambda) { gev_xi <- gpd_par["xi"] + (gpd_par["alpha"] / gpd_par["k"]) * (1- lambda ^-gpd_par["k"]) # scale gev_alpha <- gpd_par["alpha"] * lambda ^-gpd_par["k"] # scale gev_k <- gpd_par["k"] print(gev_xi)print(gev_alpha)print(gev_k)c(xi =as.numeric(gev_xi),alpha =as.numeric(gev_alpha),k =as.numeric(gev_k))}# Langbein adjustment defines the relationship between empirical flood expectancies in the partial duration series and the probability of the correspoding flood as an annual series# p_pds: partial duration non-exceedance probability# lambda: rate of events# returns the annual non-exceedance probability# Reference: Langbein, W. B. (1949). Annual Floods and the Partial-duration Flood Series. American Geophysical Union, 30(6), 879-881.langbein_adjustment <-function(lambda, p_pds) {exp(-lambda * (1- p_pds))}# Gringorten plotting position, non-exceedance probabilitygringorten_pp <-function(x) { n <-length(x) i <-rank(x) (i -0.44) / (n +0.12)}
Download Streamflow Data using dataRetrieval
Download daily streamflow data for:
Iowa River at Marshalltown, IA (05451500)
Iowa River at Marengo, IA (05453100)
NWIS is scheduled for decommission but I’m rolling with these functions for as long as I can.
# Create data frame with obs and time columnspot_marengo <-data.frame(obs = marengo_daily$Flow,time =as.numeric(marengo_daily$Date) # Convert to numeric)da_marengo <-2794# square milestime_sep <-B17_days_between(da_marengo)thresh_marengo <-10000# cfs# Extract partial duration seriespds_marengo <- POT::clust(data = pot_marengo,u = thresh_marengo, # thresholdtim.cond = time_sep, # time separationclust.max =TRUE) # extract cluster maximasummary(pds_marengo)
time obs idx
Min. :-3937 Min. :10100 Min. : 904
1st Qu.: 1908 1st Qu.:11600 1st Qu.: 6749
Median : 7517 Median :13500 Median :12358
Mean : 7849 Mean :15575 Mean :12690
3rd Qu.:14365 3rd Qu.:17100 3rd Qu.:19206
Max. :20300 Max. :48200 Max. :25140
# Create data frame with obs and time columnspot_marshalltown <-data.frame(obs = marshalltown_daily$Flow,time =as.numeric(marshalltown_daily$Date) # Convert to numeric)da_marshalltown <-1532# square milestime_sep <-B17_days_between(da_marshalltown)thresh_marshalltown <-5500# cfs# Extract partial duration seriespds_marshalltown <- POT::clust(data = pot_marshalltown,u = thresh_marshalltown, # thresholdtim.cond = time_sep, # time separationclust.max =TRUE) # extract cluster maximasummary(pds_marshalltown)
time obs idx
Min. :-24324 Min. : 5520 Min. : 241
1st Qu.: -6296 1st Qu.: 6470 1st Qu.:12422
Median : 4825 Median : 8085 Median :23543
Mean : 3577 Mean : 9389 Mean :22459
3rd Qu.: 13985 3rd Qu.:10375 3rd Qu.:32703
Max. : 20548 Max. :39400 Max. :39206
Compute rank correlation between flows at Marshalltown and Marengo using Kendall’s rank correlation coefficient (Kendall’s tau) over a range of lags.
I computed Kendall’s tau under 2 scenarios:
Rank correlation between the PDS of Marengo flows and lagged daily flows at Marshalltown
Rank correlation between the Marshalltown flows and lagged daily flows at Marengo
In both scenarios, the Kendall’s tau value was highest at a lag of 0 (flows on the same date).
Extreme and Coincident Flows
lag <-5# days# Prepare data frame with 3 columns: Date, Variable1, Variable2# For auto-correlation, use same variable twicemarengo_con <-left_join(df_pds_marengo, marshalltown_daily, by ="Date") %>%select(Date, Flow.x, Flow.y) %>%rename(Marengo = Flow.x,Marshalltown = Flow.y) %>%na.omit()head(marengo_con)
# Calculate Kendall's tau at different lagsresult <- MultiHazard::Kendall_Lag(Data = marengo_con,Lags =seq(-lag, lag, 1), # lag is applied to second named variablePLOT =TRUE, # display plotGAP =0.1) # length of y-axis above and below max and min
# Extract resultstau_values <- result$Valuep_values <- result$Test# View tau coefficientsprint(tau_values)
# Prepare data frame with 3 columns: Date, Variable1, Variable2# For auto-correlation, use same variable twicemarshalltown_con <-left_join(df_pds_marshalltown, marengo_daily, by ="Date") %>%select(Date, Flow.x, Flow.y) %>%rename(Marshalltown = Flow.x,Marengo = Flow.y) %>%na.omit()head(marshalltown_con)
# Calculate Kendall's tau at different lagsresult <- MultiHazard::Kendall_Lag(Data = marshalltown_con,Lags =seq(-lag, lag, 1), # lag is applied to second named variablePLOT =TRUE, # display plotGAP =0.1) # length of y-axis above and below max and min
# Extract resultstau_values <- result$Valuep_values <- result$Test# View tau coefficientsprint(tau_values)
The maximum tau in both scenarios occurs with a lag of 0, meaning that the coincident events should be pulled from the same date as the extremal (PDS) events.
Scenario A (Extremal Marshalltown Flows) Marginal Distributions
Scenario A is defined by extremal events (PDS) at the Marshalltown gage and the coincident events at the Marengo gage. The conditioning variable is Marshalltown flows, which are conditioned on Marengo flows.
Kendall’s tau is max at a lag of 0 days. Therefore, use left_join to match coincident events to the date of the PDS events.
# Fit marginal distributionsmarshalltown_gpa <- lmom::pelgpa(lmom::samlmu(marshalltown_con$Marshalltown)) marshalltown_gev <-madsen_annualization(marshalltown_gpa, rate_marshalltown)
xi
6250.71
alpha
3605.862
k
-0.09533262
marengo_glo <- lmom::pelglo(lmom::samlmu(marshalltown_con$Marengo))#lmom::evplot(marshalltown_con$Marshalltown, main = "Marshalltown"); #lmom::evdistq(lmom::quagev, para = marshalltown_gev)marshalltown_con_gev <- marshalltown_con %>%mutate(marshalltown_cdf = lmom::cdfgev(Marshalltown, marshalltown_gev),marengo_cdf = lmom::cdfglo(Marengo, marengo_glo))marshalltown_con_gpa <- marshalltown_con %>%mutate(marshalltown_cdf = lmom::cdfgpa(Marshalltown, marshalltown_gpa),marengo_cdf = lmom::cdfglo(Marengo, marengo_glo))ggplot(marshalltown_con_gev, aes(x = marshalltown_cdf, y = marengo_cdf)) +geom_point() +xlab("CDF of Marshalltown PDS Flow") +ylab("CDF of Marengo Coincident Flow (cfs)") +ggtitle("Iowa River") +theme_bw()
f <-seq(0.001, 0.999, length.out =1000) # nonexceedance probabilitiesq <- lmom::quaglo(f, marengo_glo)plot(qnorm(f), q, type ="l", lwd =2,log ="y",xlab ="Probability scale (qnorm(F))",ylab ="Quantile (x)",main ="Marengo Coincident Events Generalized Logistic Distribution")grid()
x <- marshalltown_con$Marshalltownp_pds <-gringorten_pp(x) # partial duration space non-exceedance probabilitymain_txt <-"Marshalltown GEV Distribution"p <-langbein_adjustment(rate_marshalltown, p_pds)plot(qnorm(p), x, xlab ="Reduced variate, -log(-log(F))", ylab ="Quantile", main ="")title(main = main_txt, line =3) pp <-seq(0.001, 0.999, length.out =200)qfit <- lmom::quagev(pp, marshalltown_gev)lines(qnorm(pp), qfit, col="blue", lwd=2)# Add a top axis with Return Period (years)T_ticks <-c(2, 5, 10, 25, 50, 100, 200) p_ticks <-1-1/ T_ticks # convert T -> nonexceedance probaxis(side =3,at =qnorm(p_ticks),labels = T_ticks)mtext("Return period (years)", side =3, line =2)
Scenario B (Extremal Marengo Flows) Marginal Distributions
Scenario B is defined by extremal events (PDS) at the Marengo gage and the coincident events at the Marshalltown gage. The conditioning variable is Marengo flows, which are conditioned on Marshalltown flows.
# Select the Copula with the minimum AIC (in this scenario, the Joe copula)marengo_cop <- joe_cop
Contour Line Values for Plotting
This code block is for Neptune
n <-60u <- v <-seq(0.02, 0.98, length.out = n)# density on gridz <-outer(u, v, function(a, b)BiCopPDF(a, b,family = marengo_cop$family,par = marengo_cop$par,par2 = marengo_cop$par2))lev <-c(0.2, 0.6, 1, 1.5, 2, 3, 5, 7)contour(u, v, z, levels = lev, col ="blue",xlab ="u", ylab ="v",main =sprintf("BiCop density contours", marengo_cop$family))levcls <-contourLines(u, v, z, levels = lev)df <-do.call(rbind, lapply(seq_along(cls), function(i) {data.frame(level = cls[[i]]$level,x = cls[[i]]$x,y = cls[[i]]$y,line = i, # optional: which polyline (useful if multiple per level)pt =seq_along(cls[[i]]$x) # optional: point index along the line )}))# If you truly want only 3 columns, drop line/pt:# df <- df[, c("level","x","y")]write.csv(df, "iowa_river_contours.csv", row.names =FALSE)
Copula Simulation and Transformation
Generate Samples
Generate a large number of samples). Split the number of samples across the two conditional copula models in proportion to the number of events in each dataset.
N <-10^6n_A <-nrow(marshalltown_con) # Number of coincident events in scenario An_B <-nrow(marengo_con) # Number of coincident events in scenario Bn1 <-round((N * n_A) / (n_A + n_B), 0)n2 <- N - n1# Simulate U(0,1)^2 from each fitted copulas1 <- VineCopula::BiCopSim(n1, obj = marshalltown_cop)s2 <- VineCopula::BiCopSim(n2, obj = marengo_cop)#plot(head(s1, 1000))
Transform to Original Scale using Quantile Functions
Transform the simulated uniforms to original units using quantile functions.
# Scenario Amarshalltown_gev <- lmomco::vec2par(marshalltown_gev, type ="gev")marengo_glo <- lmomco::vec2par(marengo_glo, type ="glo")marshalltown_A <-quagev(s1[,1], marshalltown_gev)marengo_A <-quaglo(s1[,2], marengo_glo)copula_sample_A <-data.frame(Marshalltown = marshalltown_A, Marengo = marengo_A, scenario ="A")# Scenario Bmarengo_gev <-vec2par(marengo_gev, type ="gev")marshalltown_glo <-vec2par(marshalltown_glo, type ="glo")marengo_B <-quagev(s2[,1], marengo_gev)marshalltown_B <-quaglo(s2[,2], marshalltown_glo)copula_sample_B <-data.frame(Marshalltown = marshalltown_B, Marengo = marengo_B, scenario ="B")# Combine samples from both scenarios; used in KDE copula_sample <-rbind(copula_sample_A, copula_sample_B)ggplot() +geom_point(data = dplyr::slice_sample(copula_sample, n =10000), aes(x = Marshalltown, y = Marengo, color = scenario),inherit.aes =FALSE,shape =21 ) +scale_color_manual(name ="Scenario", values =c(A ="blue",B ="red" ) ) +theme_bw() +coord_fixed(1) +xlim(0, 180000) +ylim(0, 180000) +labs(x ="Marshalltown Flow (cfs)", y ="Marengo Flow (cfs)")
Warning: Removed 10 rows containing missing values or values outside the scale range
(`geom_point()`).
Isoline Derivation
Generate Unit Square
# Using "low" resolution here to reduce run timesx <-c(10^(-5), seq(999.9*10^(-5), 1-(1*10^(-6)), 10^(-4)))y <-c(10^(-5), seq(999.9*10^(-5), 1-(1*10^(-6)), 10^(-4)))grid <-list(x = x, y = y)u <-expand.grid(x, y)
Warning in sprintf("", n_x, n_y): 2 arguments not used by format ''
Evaluate Copula CDF
Evaluate the copula CDF at each point on the grid.
# Evaluates CDF of the specified parametric copulauCopA <- VineCopula::BiCopCDF(u[,1], u[,2], marshalltown_cop) uCopB <- VineCopula::BiCopCDF(u[,1], u[,2], marengo_cop)
User-specified Joint Return Periods
# Return periodsrp_vec <-c(10, 50, 100, 200)
Determine Return Period
Determine the return period of each point on the unit square.
This code block defines the sampling frequency for daily data and a function to compute the return period surface matrix for an x, y grid in column blocks.
mu <-365.25# Sampling frequency for daily data# Computes return period surface matrix z on a grid of x and y values in blocks of columnsmake_z_surface_blocked <-function(x, y, Cmat, mean_waiting_time, block =500L){ nx <-length(x); ny <-length(y) z <-matrix(NA_real_, nrow = nx, ncol = ny) # Pre-allocate output matrix z# Iterate over the y-axis in chunksfor (j0 inseq.int(1L, ny, by = block)) { j1 <-min(j0 + block -1L, ny) yb <- y[j0:j1]# denom is nx × nb (nb = block size) denom <-1- x -rep(yb, each = nx) +as.vector(Cmat[, j0:j1, drop=FALSE]) z[, j0:j1] <- mean_waiting_time /matrix(denom, nrow = nx, ncol =length(yb)) } z}
Compute Contour Lines
compute_contour_lines <-function( return_periods, x, y, z, mu,stop_on_missing =TRUE) { purrr::map_dfr(return_periods, function(rp) { xy <- grDevices::contourLines(x, y, z, levels = mu * rp)if (length(xy) ==0) { msg <-paste("No contour found for RP =", rp)if (stop_on_missing) stop(msg) else { warning(msg); return(NULL) } }do.call(rbind, lapply(seq_along(xy), function(i) {data.frame(id = i,level = xy[[i]]$level,marshalltown_cdf = xy[[i]]$x,marengo_cdf = xy[[i]]$y,rp = rp ) })) })}
Scenario A Contour Lines
# Mean waiting time is the average time, in days, between eventsmean_waiting_time_A <- mu / rate_marshalltown # Reshape into a matrix: rows correspond to x, columns correspond to y# IMPORTANT: expand.grid(x,y) varies x fastest, then yCmat <-matrix(uCopA, nrow =length(x), ncol =length(y), byrow =FALSE)z <-make_z_surface_blocked(x, y, Cmat, mean_waiting_time_A, block=200L)contourLines_A <-compute_contour_lines(return_periods = rp_vec,x = x,y = y,z = z,mu = mean_waiting_time_A,stop_on_missing =TRUE# set FALSE to skip missing RPs with a warning instead)
Scenario B Return Period
# Mean waiting time is the average time, in days, between eventsmean_waiting_time_B <- mu / rate_marengo# Reshape into a matrix: rows correspond to x, columns correspond to y# IMPORTANT: expand.grid(x,y) varies x fastest, then yCmat <-matrix(uCopB, nrow =length(x), ncol =length(y), byrow =FALSE)z <-make_z_surface_blocked(x, y, Cmat, mean_waiting_time_B, block=200L)contourLines_B <-compute_contour_lines(return_periods = rp_vec,x = x,y = y,z = z,mu = mean_waiting_time_B,stop_on_missing =TRUE# set FALSE to skip missing RPs with a warning instead)
Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_path()`).
# View 1/100 joint probability contour from scenario A# contourLines[contourLines$scenario == "A" & contourLines$rp == "100", ]
Compute Envelope Curve
Combine Isolines into Envelope Curve
The make_envelope function takes a data frame containing isoline points for one return period, rp, but multiple curves, distinguished by scenario, and produces a new curve that is the upper envelope (for each x-value, it keeps the maximum y-value across all scenarios).
make_envelope <-function(df_rp, n_grid =500L) {# Clean and sort input # Drop rows where marshalltown (x) or marengo (y) are missing/non-finite# Sorts by marshalltown values df_rp <- df_rp %>%filter(is.finite(marshalltown), is.finite(marengo)) %>%arrange(marshalltown)# Create a common x-grid over the range covered by the isolines# Builds xg, a sequence from the minimum to maximum marshalltown in the data, with n_grid points xg <-seq(min(df_rp$marshalltown), max(df_rp$marshalltown), length.out = n_grid)# Interpolate each scenario separately (safer than interpolating across mixed curves)# Splits data by scenario# For each scenario:# Sorts by marshalltown# Removes duplicate x-values so interpolation behaves well# Use approx() to linearly interpolate marengo values at every xg# rule = 1 means no extrapolation (any xout < min(x) or xout > max(x) returns NA) y_by_curve <- df_rp %>%group_by(scenario) %>%group_split() %>%map(function(d) { d <- d %>%arrange(marshalltown)# approx requires unique x to behave well d <- d %>%distinct(marshalltown, .keep_all =TRUE)approx(x = d$marshalltown, y = d$marengo, xout = xg,rule =1, ties ="ordered")$y })# Compute the envelope = pointwise maximum across curves yenv <-do.call(pmax, c(y_by_curve, list(na.rm =TRUE)))# Return a tabletibble(rp =unique(df_rp$rp),aep =unique(df_rp$aep),marshalltown = xg,marengo = yenv) %>% dplyr::filter(marshalltown >=0, marengo >=0)}
Apply make_envelope function to combine isolines in contourLines into an envelope curve.
Sim_Max: truncation multiplier (same intent as your original function)
With very large simulations, heavy-tailed marginals can generate rare but enormous values. Sim_Max applies a guardrail by keeping the KDE focused on the range relevant to your contour/isoline or observed data.
kde_along_isoline <-function(Iso, cop_sample,Data =NULL,Source =c("Sample", "Data"),Sim_Max =10,gridsize =c(200, 200)) { Source <-match.arg(Source)# Keep exactly these rows and order Iso2 <- Iso %>% dplyr::select(marshalltown, marengo)# eval.points must be an n×2 numeric matrix for ks::kde Iso_mat <-as.matrix(Iso2)storage.mode(Iso_mat) <-"double"# ---- Print Iso_mat (or a subset if it's huge) ----#message("Iso_mat preview (first 10 rows):")#print(utils::head(Iso_mat, 10))#message("Iso_mat dimensions: ", paste(dim(Iso_mat), collapse = " x "))# ---------------------------------------------------if (Source =="Sample") { xdat <- cop_sample %>% dplyr::select(Marshalltown, Marengo) %>%as.matrix()storage.mode(xdat) <-"double" ref <-if (!is.null(Data)) Data %>% dplyr::select(Marshalltown, Marengo) else Iso2 remove <-which( xdat[,1] > Sim_Max *max(ref[[1]], na.rm=TRUE) | xdat[,2] > Sim_Max *max(ref[[2]], na.rm=TRUE) )if (length(remove) >0) xdat <- xdat[-remove, , drop =FALSE] } else {if (is.null(Data)) stop("When Source='Data', supply `Data`.") xdat <- Data %>% dplyr::select(Marshalltown, Marengo) %>%na.omit() %>%as.matrix()storage.mode(xdat) <-"double" }# ---- KDE fit (do this ONCE) ---- fit <- ks::kde(x = xdat, eval.points = Iso_mat, gridsize = gridsize) est <- fit$estimateif (length(est) !=nrow(Iso_mat)) {stop("KDE length mismatch: got ", length(est), " but Iso has ", nrow(Iso_mat), " rows.") } rng <-range(est, na.rm =TRUE) scaled <-if (is.finite(rng[1]) &&is.finite(rng[2]) && rng[2] > rng[1]) { (est - rng[1]) / (rng[2] - rng[1]) } elserep(NA_real_, length(est))print(fit$H, digits =6) ev <-eigen(fit$H)#print(ev$values, digits = 6)#print(ev$vectors, digits = 6)list(prediction = est,contour = scaled,H = fit$H,eigen = ev )}
Find Design Event Along Each Isoline
Fine “design events” or “most likely events” on an isoline. The “most likely event” is the isoline point with the maximum KDE.
Iso: data.frame with marshalltown, marengo columns
prediction: KDE values returned by kde_along_isoline()$prediction
N_Ensemble: number of sampled ensemble points along isoline (weighted by prediction)
extract_design_events_2d <-function(Iso, prediction, N_Ensemble =0) { Iso_mat <- Iso %>% dplyr::select(marshalltown, marengo) %>%as.data.frame()# Most likely = max KDE point on isoline idx <-which(prediction ==max(prediction, na.rm =TRUE))[1] most <-data.frame(marshalltown = Iso_mat[idx, 1],marengo = Iso_mat[idx, 2] )# Keep your previous "full dependence" behavior (max x and max y on isoline) full <-data.frame(marshalltown =max(Iso_mat[, 1], na.rm =TRUE),marengo =max(Iso_mat[, 2], na.rm =TRUE) ) ens <-NULL# Initialize ensemble output as NULLif (N_Ensemble >0) {# Find indices of isoline points where the KDE values are finite and strictly positive ok <-which(is.finite(prediction) & prediction >0)# If there are no valid KDE weights, return an empty 0-row matrix/data frame with the same number of columns as Iso_matif (length(ok) ==0) { ens <- Iso_mat[0, , drop =FALSE]# Randomly select N_Ensemble indices from the valid set ok } else {# prob = prediction[ok] makes it weighted sampling: higher KDE density -> more likely to be chosen draw <-sample(ok, size = N_Ensemble, replace =TRUE, prob = prediction[ok])# Extracts the sampled isoline coordinates corresponding to the chosen indices ens <- Iso_mat[draw, , drop =FALSE]# Label the two columnsnames(ens) <-c("marshalltown", "marengo") } }list(MostLikelyEvent = most, FullDependence = full, Ensemble = ens)}
KDE and Design Event
compute_isoline_results <-function( return_periods, envelope_curves, copula_sample,N_Ensemble =200,Source ="Sample",Sim_Max =10,gridsize =c(1000, 1000)) {# ── 1. Filter to requested RPs and split into a list of data frames ────────── iso_list <- envelope_curves %>% dplyr::filter(rp %in% return_periods) %>% dplyr::group_by(rp) %>% dplyr::group_split()# Name each element by its RP value for easy downstream access rp_vals <- purrr::map_dbl(iso_list, ~unique(.x$rp))names(iso_list) <-as.character(rp_vals)# ── 2. For each RP: KDE along isoline + most-likely design event ───────────── purrr::map(iso_list, function(df_iso) { kde_out <-kde_along_isoline(Iso = df_iso,cop_sample = copula_sample,Source = Source,Sim_Max = Sim_Max,gridsize = gridsize )stopifnot(length(kde_out$prediction) ==nrow(df_iso)) design <-extract_design_events_2d( df_iso, kde_out$prediction,N_Ensemble = N_Ensemble )list(isoline = df_iso %>% dplyr::mutate(kde = kde_out$prediction,kde_scaled = kde_out$contour ),kde_out = kde_out,design = design,design_event = design$MostLikelyEvent ) })}# ── Run everything in one call ─────────────────────────────────────────────────iso_results <-compute_isoline_results(return_periods = rp_vec,envelope_curves = envelope_curves,copula_sample = copula_sample,N_Ensemble =200,gridsize =c(400, 400))
Warning in kdde.binned.nd(H = H, deriv.order = r, bin.par = bin.par, verbose =
verbose, : Binning grid may be too coarse for current (small) bandwidth:
consider increasing grid size for dimensions 1, 2
Warning in kdde.binned.nd(H = H, deriv.order = r, bin.par = bin.par, verbose =
verbose, : Binning grid may be too coarse for current (small) bandwidth:
consider increasing grid size for dimensions 1, 2
# ── Flatten all KDE-annotated isolines into one tidy data frame ───────────────contourLines_kde <- purrr::map_dfr( iso_results,~ .x$isoline,.id ="rp"# adds a character column "rp") %>% dplyr::mutate(rp =as.numeric(rp),aep =1/rp)# ── Flatten all design events into one data frame ─────────────────────────────design_events_all <- purrr::map_dfr( iso_results,~ .x$design_event,.id ="rp") %>% dplyr::mutate(rp =as.numeric(rp))
