Extreme value distribution parameters were sourced from Canute 3 with reference to O’Grady et al 2019 to estimate Annual Exceedance Probabilities (AEPs) and AEP multiplication factors. Tide gauge estimates have been updated to the latest version of GESLAv3. Multiplication factors are explained in Canute 3

The code pipeline generates pdf images in the plots dir of this git repo, and csv files in the data dir. 10th and 90th percentile estimates were made using 1.28 times the sqrt of the standard error of the scale parameter fit.

File names are in the format <metric>_<domain>_<multimodel_allpercentiles>_<SLR>.csv:

CSV file columns represent:

GESLAv3 csv file columns also include the TG name and the number of years with at least 80% of hourly observations (nCompleteYr).

The NCRA will explore SLR increments of 0.06, 0.1, 0.2, 0.38, 0.6 and 1.0 m

A value of SLR 0.06 extrapolates the current SLR trajectory from the IPCC AR6 baseline to 2020. Below is a table of global warming level (GWL) match ups to SSP to inform the SLR increments, based on IPCC Tables 9.9 and 9.10 in Fox-Kemper et al 2021.

Year GWL Table 9.10 (Indicative SSP for GWL) Table 9.9 (Closest SLR and Year) Table 9.9 Values Used for SLR Inc.
2030 1.5 SSP1-2.6 SSP1-2.6 0.09 (0.08–0.12)
2050 2 SSP1-2.6/SSP2-4.5 SSP2-4.5 0.20 (0.17–0.26)
2090 2 SSP1-2.6/SSP2-4.5 SSP1-2.6 0.39 (0.30–0.54)
2090 3 SSP2-4.5/SSP3-7.0 SSP5-8.5 0.63 (0.52–0.83)
2100 5 SSP5-8.5 SSP5-8.5 0.70 (0.63–1.01) 
#install.packages("terra")
require(terra)
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.4.1
## terra 1.7.78
#choose input SLR increments from AR6 
SLRs <- c(0.06, 0.1, 0.2, 0.38, 0.6, 1.0)
#also need increments form NCRA baseline
SLRs <- c(SLRs,SLRs[-1]-SLRs[1],0,0.9)
#SLRs <- c(0.09,0.20)

outdir = "plots/"
colz = rev(c("#000078","#0000BB","#0000FF","#0C3AFC","#1975FA","#25B0F7","#32EBF5","#73EEE6","#B4F1D7","#F5F5C8",
             "#F5F564","#F5E032","#F6CC00","#F7A300","#F97A00","#FB5100","#FF0000","#DA0006","#B5000D","#910014",
             "#7B003D","#660066","#88008C","#AB00B2","#CD00D8","#F000FF","#CAA6DE","#DBA7DD","#EDA9D2"))#,"#FFABC7") )


ncolz = length(colz)
colzRamp = colorRampPalette(rev(colz))
roms = vect("ROMS_gumbel_parameters_SWLandMTWL_Australia")

tg = vect("gesla3_fgumbel_shp")
tg = tg[tg$nCompleteY >= 20]

#Align the ROMs 1 year ARI (63% AEP) location parameter to the nearest tide gauge value following O'Grady et al 2019 and Haigh et al. 2014.
crs(roms) = crs(tg)
nb = nearby(roms,tg)
roms$muOff = tg$mle_1[nb[,2]]
#Define a region for each point
#regions = project(vect("Coastal_Waters_areas_(AMB2020)"),roms)
#regions$state = c("VIC","TAS","QLD","NT","WA","SA","NSW")

regionsb = project(vect("NCRA_regions_coastal_waters_GDA94"),roms)
regionsb$short_name = regionsb$NCRA
regionsc = project(vect("ncra_regions"),roms)

regionsl = lapply(regionsb$short_name, function(x) {
  y = aggregate(union(regionsb[regionsb$short_name == x],regionsc[regionsc$short_name == x]));
  y$state = x;y})
regions = vect(regionsl)
regions$state
## [1] "NSW"      "NT"       "QLDNorth" "QLDSouth" "SA"       "TAS"      "VIC"     
## [8] "WANorth"  "WASouth"
plot(regions,"state")

ints = relate(roms,regions,"intersects",pairs=TRUE)
roms$state = NA
roms$state[ints[,1]] = regions$state[ints[,2]]
plot(roms,"state")

lgas = vect("LGA_2024_AUST_GDA2020")
lgas$id = 1:length(lgas)
roms$id = 1:length(roms)
romsp = project(roms,lgas)
lgas_coast = buffer(as.lines(buffer(terra::aggregate(lgas),0)),5)
lgas_coast
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 0  (geometries, attributes)
##  extent      : 96.81691, 167.9981, -43.74054, -9.142117  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat GDA2020 (EPSG:7844)
plot(lgas_coast,col=2)

ints = intersect(lgas,lgas_coast)
lga_coast_cp = centroids(ints,inside = TRUE)
gg = geom(lga_coast_cp)
out = data.frame(lga = lga_coast_cp$LGA_NAME24,longitude = gg[,3],latitude = gg[,4])
write.csv(out,file = paste0("data/Coastal_lga_5m_centroid.csv"))
require(sf)
## Loading required package: sf
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
nearest_lgas <- function(romsp, lgas) {
  # Convert SpatVector objects to sf objects
  romsp_sf <- st_as_sf(romsp)
  lgas_sf <- st_as_sf(lgas)
  
  # Find the nearest LGAs for each romsp feature
  nearest_indices <- st_nearest_feature(romsp_sf, lgas_sf)
  
  # Return the corresponding nearest LGA geometries or data
  nearest_lgas <- lgas_sf[nearest_indices, ]
  
  return(nearest_lgas)
}

result <- nearest_lgas(romsp, lgas)
roms$lga = result$LGA_NAME24
plot(roms,"lga")
lines(lgas)

#length(unique(roms$lga))

#writeVector(roms,"roms/roms.shp",overwrite=TRUE)

Confidence intervals

#90 percent confidence intervals around the maximum likelihood estimate (MLE) GEV or MC EVD were computed from the quantiles of 100 Monte Carlo simulations of the fitted model.
#larger nMonteCarloSims results in smaller confidence bounds.
gumbel_mleAndConfIntvals = function(loc,scale,ARI = 1:500,nMonteCarloSims = 100,myseed = 19811019){
  require(evd)
  set.seed(myseed)
  if(length(ARI) == 1) ARIs = c(ARI-.9,ARI-.5,ARI,ARI+.5,ARI+.9)
  if(length(ARI) > 1) ARIs = ARI
  ARI_max = max(ARIs)
  nSampes = nMonteCarloSims
  samp <- rgev(nSampes * 99, loc = loc, scale = scale,shape = 0)
  samp <- matrix(samp, nSampes, 99)
  samp <- apply(samp, 2, sort)
  samp <- apply(samp, 1, sort)
  env <- t(samp[c(5,50, 95), ])
  xx <- -1/log((1:nSampes)/(nSampes + 1))
  newdata <-  -1/log((ARIs)/(ARI_max + 1))

  fits = data.frame(logxx = log(xx),env1=env[, 1],env2=env[, 2],env3=env[, 3])
  #probablity fitted with loess function
  p05 = predict(loess(env1 ~ logxx,fits,control = loess.control(surface = "direct")),newdata = data.frame(logxx = log(newdata)))
  p50 = predict(loess(env2 ~ logxx,fits,control = loess.control(surface = "direct")),newdata = data.frame(logxx = log(newdata)))
  p95  = predict(loess(env3 ~ logxx,fits,control = loess.control(surface = "direct")),newdata = data.frame(logxx = log(newdata)))
  #exceedence probablity (1-probability), linear interpolate to input ARIs

  AEPs = 1-exp(-1/ARI)
  mle  = qgumbel(1-AEPs,loc =  loc, scale =  scale)

  ep95 = approx(newdata,p05,ARI)$y
  ep50 = approx(newdata,p50,ARI)$y
  ep05 = approx(newdata,p95,ARI)$y

  return(data.frame(ARIs = ARI,ep05 = ep05-ep50+mle,ep50=mle,ep95=ep95-ep50+mle))
}

SWL_rls = t(sapply(1:length(roms), function(x) t(gumbel_mleAndConfIntvals(loc = roms$muOff[x],scale = roms$SWLlam[x],ARI = 100,nMonteCarloSims = 100))))
## Loading required package: evd
MTWL_rls = t(sapply(1:length(roms), function(x) t(gumbel_mleAndConfIntvals(loc = roms$muOff[x]-roms$SWLmu[x]+roms$MTWLmu[x],scale = roms$MTWLlam[x],ARI = 100,nMonteCarloSims = 100))))

gg = terra::geom(roms)
out = data.frame(state = roms$state,lga = roms$lga,longitude = gg[,3],latitude = gg[,4],
                 SWL_pc5 = SWL_rls[,4], SWLpc50 = SWL_rls[,3], SWLpc95 = SWL_rls[,2],
                 MTWL_pc5 = MTWL_rls[,4], MTWLpc50 = MTWL_rls[,3], MTWLpc95 = MTWL_rls[,2])

saout = out[out$state == "SA",]
write.csv(saout,"data/SA_ROMS_CAWCR_OGRADY_ET_AL_2019_1pcAEP.csv",row.names = FALSE)

tgsa = crop(tg,ext(roms[roms$state == "SA"]))

tg_SWL_rls = t(sapply(1:length(tgsa), function(x) t(gumbel_mleAndConfIntvals(loc = tgsa$mle_1[x],scale = tgsa$mle_2[x],ARI = 100,nMonteCarloSims = 100))))

gg = terra::geom(tgsa)
tgout = data.frame(name = tgsa$Name,nCompleteYears = tgsa$nCompleteY,longitude = gg[,3],latitude = gg[,4],
                 SWL_pc5 = tg_SWL_rls[,4], SWLpc50 = tg_SWL_rls[,3], SWLpc95 = tg_SWL_rls[,2])

write.csv(tgout,"data/SA_GESLA3_TIDEGAUGES_OGRADY_ET_AL_2019_1pcAEP.csv",row.names = FALSE)
#predict return levels
AEP = 0.01
return_level <- function(mu,lam,AEP){
  z = mu+lam*(-log(-log(1-AEP)))
  return(z)
}


#storm tide return level
roms$SWL_RL = return_level(roms$muOff,roms$SWLlam,AEP)
#storm tide plus wave setup
roms$MTWL_RL = return_level(roms$muOff-roms$SWLmu+roms$MTWLmu,roms$MTWLlam,AEP)

tg$SWL_RL = return_level(tg$mle_1,tg$mle_2,AEP)

#compute MTWL
ats = seq(0,7,0.5)
plotit = function(){
  plot(roms,"SWL_RL",breaks=ats,range = range(ats), main = paste0("1% AEP extreme SWL (reletive to MSL of that year)"),type = "continuous",col = colzRamp(length(ats)-1))
  plot(tg,"SWL_RL",add = TRUE,legend=FALSE,col = colzRamp(length(ats)-1),cex = 1.5,breaks=ats,range = range(ats))
  plot(tg,cex = 1.5,pch=1,add=TRUE)
}

plotit()

pdf(paste0(outdir,"One_Percent_AEP_SLW.pdf"),width=12)
plotit()
dev.off()
## png 
##   2

Figure : Height (in m relative to mean sea level in a given year) of the 1% annual exceedance probability (AEP) extreme still water level (SWL) storm tide (Source: O’Grady et al 2019) overlayed here with GESLA3 tide gauge record estimates (Source Canute 3). This shows the height of a water level from the combination of storm surge and tide with a 1% chance of occurring, these heights can be offset with a mean sea level reference to a given year and vertical datum (e.g. AHD) and can be generated for any SLR or AEP for NCRA. Heights can also be generated with the inclusion of damaging wave effect (wave setup O’Grady et al 2019) for sandy beach coastlines. Here, levelled heights can be integrated with land height data to map coastal flooding. Similar plots based on SCHISM data will be developed by ACS by June 2025.

ats = seq(0,7,0.5)
plotit = function() plot(roms,"MTWL_RL",breaks=ats,range = range(ats), main = paste0("1% AEP extreme MTWL (reletive to MSL of that year)"),type = "continuous",col = colzRamp(length(ats)-1))

plotit()

pdf(paste0(outdir,"One_Percent_AEP_MTWL.pdf"),width=12)
plotit()
dev.off()
## png 
##   2

Figure : Height (in m relative to mean sea level in a given year) of the 1% annual exceedance probability (AEP) extreme mean total water level (MTWL) storm tide plus wave setup (Source: O’Grady et al 2019, Canute 3). This shows the height of a water level from the combination of storm surge and tide and wave setup for sandy beaches with a 1% chance of occurring, these heights can be offset with a mean sea level reference to a given year and vertical datum (e.g. AHD) and can be generated for any SLR or AEP for NCRA. Heights can also be generated without the inclusion of damaging wave effect (still water level O’Grady et al 2009). Here, levelled heights can be integrated with land height data to map coastal flooding. Similar plots based on SCHISM data will be developed by ACS by June 2025.

# Set initial Sea Level Rise value
SLR <- 0.6

# Function to plot multiplication factors
plotMF <- function(v, att, ats, buf = 1, SLR = 0.6, ...) {
  v$SLR <- SLR
  
  if ("SWLlam" %in% names(v)) {
    v$SWL_MFSLR_50 <- exp(SLR / v$SWLlam)
    v$SWL_MFSLR_10 <- exp(SLR / (v$SWLlam + 1.28 * v$SWLlam_se))
    v$SWL_MFSLR_90 <- exp(SLR / (v$SWLlam - 1.28 * v$SWLlam_se))
    
    v$MTWL_MFSLR_50 <- exp(SLR / v$MTWLlam)
    v$MTWL_MFSLR_10 <- exp(SLR / (v$MTWLlam + 1.28 * v$MTWLlam_se))
    v$MTWL_MFSLR_90 <- exp(SLR / (v$MTWLlam - 1.28 * v$MTWLlam_se))
  }
  
  if ("mle_2" %in% names(v)) {
    se = sqrt(tg$cov_22)
    v$SWL_MFSLR_50 <- exp(SLR / tg$mle_2)
    v$SWL_MFSLR_10 <- exp(SLR / tg$mle_2 + 1.28 * se)
    v$SWL_MFSLR_90 <- exp(SLR / tg$mle_2 - 1.28 * se)
    
  }
  
  logats <- log(ats)
  
  extu <- as.numeric(as.vector(ext(v)))
  dat <- cut(as.numeric(unlist(values(v[att]))), breaks = ats)
  v$plot <- dat
  levs <- levels(dat)
  
  dumy <- terra::vect(array(0, dim = c(length(ats), 2)), "points")
  dumy$plot <- cut(ats, breaks = c(ats))
  dumy <- rbind(dumy, v)
  dumy <- dumy[order(dumy$plot, decreasing = TRUE)]
  
  terra::plot(dumy, "plot", levels = levs, sort = FALSE, 
              col = rev(viridis::viridis(length(ats - 1))),
              xlim = extu[1:2] + c(-buf, buf), ylim = extu[3:4] + c(-buf, buf), ...)
  
  return(v)
}

# Function to create plots for specific SLR values
plotit <- function(SLR) {
  par(mfrow = c(1, 2))
  
  dat <- plotMF(v = roms, att = "SWL_MFSLR_50", ats = ats, SLR = SLR, 
                main = paste0("Extreme SWL multiplication factor for ", SLR, "m of SLR"), 
                mar = c(2, 2, 2, 4) + 0.1)
  
  ge <- geom(dat)
  
  out <- data.frame(
    state = dat$state,
    lga = dat$lga,
    longitude = ge[, 3], 
    latitude = ge[, 4], 
    SLR = dat$SLR, 
    SWL_AEP1PC = round(dat$SWL_RL + dat$SLR, 2), 
    MTWL_AEP1PC = round(dat$MTWL_RL + dat$SLR, 2), 
    SWL_MFSLR_50 = round(dat$SWL_MFSLR_50, 2), 
    SWL_MFSLR_10 = round(dat$SWL_MFSLR_10, 2), 
    SWL_MFSLR_90 = round(dat$SWL_MFSLR_90, 2), 
    MTWL_MFSLR_50 = round(dat$MTWL_MFSLR_50, 2), 
    MTWL_MFSLR_10 = round(dat$MTWL_MFSLR_10, 2), 
    MTWL_MFSLR_90 = round(dat$MTWL_MFSLR_90, 2)
  )
  write.csv(out, file = paste0(outdir, "../data/MFSLR_ACSSH_MMALL_SLR", SLR * 100, "cm.csv"), row.names = FALSE)
  
  dattg <- plotMF(v=tg, att="SWL_MFSLR_50", ats = ats, SLR = SLR, cex = 1.5, add = TRUE, legend = FALSE)
  plot(tg, cex = 1.5, SLR = SLR, pch = 1, add = TRUE)
  plotMF(v = roms, att = "MTWL_MFSLR_50", ats = ats, SLR = SLR, 
         main = paste0("Extreme MTWL multiplication factor for ", SLR, "m of SLR"), 
         mar = c(2, 2, 2, 4) + 0.1)

  ge <- geom(dattg)
  dattg <- values(dattg)
  outtg <- data.frame(
    name <- dattg$Name,
    nCompleteY <- dattg$nCompleteY,
    longitude = ge[, 3], 
    latitude = ge[, 4], 
    SLR = dattg$SLR, 
    SWL_AEP1PC = round(dattg$SWL_RL + dattg$SLR, 2), 
    SWL_MFSLR_50 = round(dattg$SWL_MFSLR_50, 2), 
    SWL_MFSLR_10 = round(dattg$SWL_MFSLR_10, 2), 
    SWL_MFSLR_90 = round(dattg$SWL_MFSLR_90, 2) 
  )
  write.csv(outtg, file = paste0(outdir, "../data/MFSLR_GESLA3_MMALL_SLR", SLR * 100, "cm.csv"), row.names = FALSE)
    
}

# Break points for the plots
ats <- c(0, 1, 2, 5, 10, 20, 50, 100, 200, Inf)

# Different SLR values to plot
plotit(SLR = 0.6)
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter

# Loop over SLR values and generate plots and CSV files
for (SLR in SLRs) {
  pdf(paste0(outdir, "MFSLR_ACSSH_MM50_SLR_", SLR * 100, "cm.pdf"), width = 14, height = 5)
  plotit(SLR=SLR)
  dev.off()
}
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter
## Warning in graphics::plot.xy(list(x = g[, 1], y = g[, 2]), type = "p", pch =
## pch, : "SLR" is not a graphical parameter

Figure : Multiplication/amplification factors (Hunter 2012) of extreme storm tide water level and mean total water level (storm tide plus wave setup) AEPs (e.g. 1% AEP in Figure A1) for an amount of SLR (see figure title), overlayed here with GESLA3 tide gauge record estimates (Source Canute 3). Here, a AEP with an associated return level is mutiplied by these factors (e.g., 1% x 10) to yeild a new AEP (i.e., 10%) for the associated return level when considering sea level rise. These return levels can be aligned with coastal structure design levels (or DEMs) to understand the change in frequency of exceedence and can be generated for any SLR or AEP for the NCRA.

# Load necessary libraries
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stringr)
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:terra':
## 
##     extract
# Define the sea level rise increments
#SLRs <- c(0.06, 0.1, 0.2, 0.38, 0.6, 1.0)
#SLRs <- SLRs[-1] - SLRs[1]

# Define the input file paths
infiles <- paste0("data/MFSLR_ACSSH_MMALL_SLR", SLRs * 100, "cm.csv")

# Initialize empty lists to store results for SWL and MTWL
results_swl <- list()
results_mtwl <- list()

# Loop over each file
for (file in infiles) {
  # Load the dataset
  data <- read.csv(file)
  
  # Extract the SLR value from the file name
  slr_value <- as.numeric(str_extract(file, "\\d+")) / 100  # Convert to float
  
  # Find the row with the median SWL_MFSLR_50 for each state
  median_rows_swl <- data %>%
    group_by(state) %>%
    slice(which.min(abs(SWL_MFSLR_50 - median(SWL_MFSLR_50)))) %>%
    mutate(SLR = slr_value)  # Use the extracted and converted SLR value
  
  # Append the SWL result to the results list
  results_swl[[file]] <- median_rows_swl
  
  # Find the row with the median MTWL_MFSLR_50 for each state
  median_rows_mtwl <- data %>%
    group_by(state) %>%
    slice(which.min(abs(MTWL_MFSLR_50 - median(MTWL_MFSLR_50)))) %>%
    mutate(SLR = slr_value)  # Use the extracted and converted SLR value
  
  # Append the MTWL result to the results list
  results_mtwl[[file]] <- median_rows_mtwl
}

# Combine all results into single data frames
final_results_swl <- bind_rows(results_swl)
final_results_mtwl <- bind_rows(results_mtwl)

# Optionally, you can save the final results to CSV files
write.csv(final_results_swl, "data/Median_SWL_MFSLR_50_by_State.csv", row.names = FALSE)
write.csv(final_results_mtwl, "data/Median_MTWL_MFSLR_50_by_State.csv", row.names = FALSE)

# Create a new table for SWL with states as rows and SLR values as columns
formatted_table_swl <- final_results_swl %>%
  mutate(SWL_MFSLR_50_Range = ifelse(SWL_MFSLR_50 > 200, 
                                     ">200 [>200 - >200]", 
                                     paste0(SWL_MFSLR_50, " [", 
                                            ifelse(SWL_MFSLR_10 > 200, ">200", SWL_MFSLR_10), " - ", 
                                            ifelse(SWL_MFSLR_90 > 200, ">200", SWL_MFSLR_90), "]"))) %>%
  select(state, SLR, SWL_MFSLR_50_Range) %>%
  pivot_wider(names_from = SLR, values_from = SWL_MFSLR_50_Range)

# Display the formatted SWL table
print(formatted_table_swl)
## # A tibble: 10 × 14
## # Groups:   state [10]
##    state    `0.06`    `0.1` `0.2` `0.38` `0.6` `1`   `0.04` `0.14` `0.32` `0.54`
##    <chr>    <chr>     <chr> <chr> <chr>  <chr> <chr> <chr>  <chr>  <chr>  <chr> 
##  1 NSW      5.83 [4.… 18.8… >200… >200 … >200… >200… 3.24 … 61.14… >200 … >200 …
##  2 NT       3.11 [2.… 6.64… 44.1… >200 … >200… >200… 2.13 … 14.17… >200 … >200 …
##  3 QLDNorth 2.79 [2.… 5.55… 30.7… >200 … >200… >200… 1.98 … 11 [7… >200 … >200 …
##  4 QLDSouth 4.52 [3.… 12.3… 152.… >200 … >200… >200… 2.73 … 33.76… >200 … >200 …
##  5 SA       1.58 [1.… 2.14… 4.57… 17.91… 95.1… >200… 1.35 … 2.9 [… 11.36… 60.36…
##  6 TAS      2.55 [2.… 4.77… 22.7… >200 … >200… >200… 1.87 … 8.92 … 148.5… >200 …
##  7 VIC      2.25 [2.… 3.87… 14.7… 170.3… >200… >200… 1.71 … 6.58 … 75.69… >200 …
##  8 WANorth  2.57 [2.… 4.81… 23.1… >200 … >200… >200… 1.87 … 9 [6.… 151.7… >200 …
##  9 WASouth  2.01 [1.… 3.21… 10.2… 83.66… >200… >200… 1.59 … 5.11 … 41.59… >200 …
## 10 <NA>     3.12 [2.… 6.66… 44.3… >200 … >200… >200… 2.16 … 14.22… >200 … >200 …
## # ℹ 3 more variables: `0.94` <chr>, `0` <chr>, `0.9` <chr>
# Save the formatted SWL table to a CSV file
write.csv(formatted_table_swl, "data/Formatted_SWL_MFSLR_50_Table.csv", row.names = FALSE)

# Create a new table for MTWL with states as rows and SLR values as columns
formatted_table_mtwl <- final_results_mtwl %>%
  mutate(MTWL_MFSLR_50_Range = ifelse(MTWL_MFSLR_50 > 200, 
                                      ">200 [>200 - >200]", 
                                      paste0(MTWL_MFSLR_50, " [", 
                                             ifelse(MTWL_MFSLR_10 > 200, ">200", MTWL_MFSLR_10), " - ", 
                                             ifelse(MTWL_MFSLR_90 > 200, ">200", MTWL_MFSLR_90), "]"))) %>%
  select(state, SLR, MTWL_MFSLR_50_Range) %>%
  pivot_wider(names_from = SLR, values_from = MTWL_MFSLR_50_Range)

# Display the formatted MTWL table
print(formatted_table_mtwl)
## # A tibble: 10 × 14
## # Groups:   state [10]
##    state    `0.06`    `0.1` `0.2` `0.38` `0.6` `1`   `0.04` `0.14` `0.32` `0.54`
##    <chr>    <chr>     <chr> <chr> <chr>  <chr> <chr> <chr>  <chr>  <chr>  <chr> 
##  1 NSW      1.3 [1.2… 1.56… 2.42… 5.39 … 14.2… 83.3… 1.19 … 1.86 … 4.13 … 10.89…
##  2 NT       1.45 [1.… 1.86… 3.46… 10.58… 41.4… >200… 1.28 … 2.38 … 7.29 … 28.56…
##  3 QLDNorth 1.5 [1.4… 1.95… 3.82… 12.77… 55.7… >200… 1.31 … 2.55 … 8.51 … 37.32…
##  4 QLDSouth 1.49 [1.… 1.95… 3.81… 12.72… 55.4… >200… 1.31 … 2.55 … 8.51 … 37.11…
##  5 SA       1.25 [1.… 1.45… 2.1 … 4.09 … 9.23… 40.6… 1.16 … 1.68 … 3.27 … 7.39 …
##  6 TAS      1.33 [1.… 1.6 … 2.57… 6.03 … 17.0… 112.… 1.21 … 1.94 … 4.54 … 12.84…
##  7 VIC      1.34 [1.… 1.63… 2.65… 6.36 … 18.5… 130.… 1.21 … 1.98 … 4.75 … 13.85…
##  8 WANorth  1.41 [1.… 1.77… 3.12… 8.74 … 30.3… >200… 1.26 … 2.22 … 6.17 … 21.79…
##  9 WASouth  1.26 [1.… 1.47… 2.16… 4.33 … 10.1… 47.3… 1.17 … 1.72 … 3.44 … 8.03 …
## 10 <NA>     1.31 [1.… 1.58… 2.46… 5.54 … 14.9… 90.6… 1.2 [… 1.88 … 4.23 … 11.4 …
## # ℹ 3 more variables: `0.94` <chr>, `0` <chr>, `0.9` <chr>
# Save the formatted MTWL table to a CSV file
write.csv(formatted_table_mtwl, "data/Formatted_MTWL_MFSLR_50_Table.csv", row.names = FALSE)
# Load necessary libraries
library(ggplot2)
library(dplyr)

new_data <- read.csv("data/Median_MTWL_MFSLR_50_by_State.csv")

# Create a PDF file to save the plots
pdf("plots/SLR_vs_MTWL_MFSLR_by_State.pdf")

# Get a list of unique states
states <- unique(new_data$state)
states = states[!is.na(states)]

par(mfrow = c(3,3),mar = c(4,4,1,1))
# Loop over each state and create a plot
for (state in states) {
  # Filter data for the current state and remove NA rows
  state_data <- new_data[new_data$state == state, ]
  state_data <- na.omit(state_data)
  
  # Set up the plot with the y-axis fixed between 0 and 200
  plot(state_data$SLR, state_data$MTWL_MFSLR_50, type = "l", col = "blue",
       ylim = c(0, 200), xlab = "SLR (Sea Level Rise)",
       ylab = "MTWL Magnification Factor", 
       main = paste(" ", state))
  
  # Add the ribbon as a shaded polygon
  polygon(c(state_data$SLR, rev(state_data$SLR)),
          c(state_data$MTWL_MFSLR_10, rev(state_data$MTWL_MFSLR_90)),
          col = rgb(0, 0, 1, alpha = 0.2), border = NA)
  
  # Add the 50th percentile line on top
  lines(state_data$SLR, state_data$MTWL_MFSLR_50, col = "blue", lwd = 2)
}


# Close the PDF device
dev.off()
## png 
##   2
# Load the dataset
new_data <- read.csv("data/Median_MTWL_MFSLR_50_by_State.csv")

# Create a PDF file to save the plots
pdf("plots/SLR_vs_SWL_MFSLR_by_State.pdf")

# Get a list of unique states
states <- unique(new_data$state)
states <- states[!is.na(states)]

# Set up the plotting area for multiple plots
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))

# Loop over each state and create a plot
for (state in states) {
  # Filter data for the current state and remove NA rows
  state_data <- new_data[new_data$state == state, ]
  state_data <- na.omit(state_data)
  
  # Set up the plot with the y-axis fixed between 0 and 200
  plot(state_data$SLR, state_data$SWL_MFSLR_50, type = "l", col = "blue",
       ylim = c(0, 200), xlab = "SLR (Sea Level Rise)",
       ylab = "SWL Magnification Factor", 
       main = paste(" ", state))
  
  # Add the ribbon as a shaded polygon
  polygon(c(state_data$SLR, rev(state_data$SLR)),
          c(state_data$SWL_MFSLR_10, rev(state_data$SWL_MFSLR_90)),
          col = rgb(0, 0, 1, alpha = 0.2), border = NA)
  
  # Add the 50th percentile line on top
  lines(state_data$SLR, state_data$SWL_MFSLR_50, col = "blue", lwd = 2)
}

# Close the PDF device
dev.off()
## png 
##   2