Bathymetry comparison

analysis/efh_bathymetry_2028/.gitignore

@@ -16,4 +16,5 @@
 /data/ai_grid_100m/
 /data/Norton_Bathymetry/
 /data/goa_bathy/
-/data/*.nc
+/data/*.nc
+/data/2022_efh_bathy/

analysis/efh_bathymetry_2028/03_evaluate_bathy.R

@@ -0,0 +1,169 @@
+library(akgfmaps)
+
+# Setup
+dir.create(here::here("analysis", "efh_bathymetry_2028", "plots", "bathy"), recursive = TRUE)
+
+region <- c("goa", "ai", "ebs")
+
+theme_2 <- function() {
+  theme_bw() %+replace%
+    theme(legend.title = element_text(size = 7, color = "black"),
+          panel.border = element_rect(color = "black", fill = NA),
+          legend.text = element_text(size = 7, color = "black"),
+          legend.key = element_rect(fill = NA, color = NA),
+          axis.title = element_blank(),
+          axis.text = element_text(size = 6, color = "black"),
+          plot.background = element_rect(fill = "white", color = NA),
+          panel.background = element_rect(fill = "white", color = NA))
+}
+
+# Get bottom depth at tow midpoints
+source(here::here("analysis", "efh_bathymetry_2028", "get_towmid.R"))
+
+
+# Comparison metrics
+mean_relative_error <- function(x, y) {
+  mean(abs((x-y))/y, na.rm = TRUE)
+}
+
+mean_absolute_error <- function(x, y) {
+  mean(abs(log10(x)-log10(y)), na.rm = TRUE)
+}
+
+mean_bias <- function(x, y) {
+  10^mean(log10(x)-log10(y), na.rm = TRUE)
+}
+
+
+# Comparison code
+
+for(ii in 1:length(region)) {
+
+  old_bathy <- terra::rast(here::here("analysis", "efh_bathymetry_2028", "data", "2022_efh_bathy", 
+                                      paste0("Variables_", region[ii], "_1km"), "Bathy.grd"))
+
+  new_bathy <- terra::rast(here::here("analysis", "efh_bathymetry_2028", "output", "efh_bathy", "efh_bathy_1km.tif")) 
+
+  new_bathy_resampled <- terra::resample(new_bathy, old_bathy, method = "bilinear")
+
+  old_bathy_reproj <- terra::project(old_bathy, new_bathy_resampled)
+
+  diff_bathy <- terra::trim(new_bathy_resampled - old_bathy_reproj)
+
+  map_layers <- akgfmaps::get_base_layers(select.region = region[ii], set.crs = terra::crs(diff_bathy))
+
+  plot_cont_bathy_diff <- ggplot() +
+    geom_stars(data = stars::st_as_stars(diff_bathy)) +
+    geom_sf(data = map_layers$akland) +
+    coord_sf(xlim = map_layers$plot.boundary$x, 
+             ylim = map_layers$plot.boundary$y) +
+    scale_fill_distiller(name = expression(Z[new]-Z[2022]), 
+                         palette = "BrBG",
+                         na.value = NA) +
+    scale_x_continuous(breaks = map_layers$lon.breaks) +
+    scale_y_continuous(breaks = map_layers$lat.breaks) +
+    theme_2()
+
+  plot_disc_bathy_diff <- ggplot() +
+    geom_stars(data = stars::st_as_stars(diff_bathy)) +
+    geom_sf(data = map_layers$akland) +
+    coord_sf(xlim = map_layers$plot.boundary$x, 
+             ylim = map_layers$plot.boundary$y) +
+    scale_fill_fermenter(name = expression(Z[new]-Z[2022]), 
+                         palette = "BrBG", 
+                         breaks = c(-Inf, -1000, -100, -10, -1, 1, 10, 100, 1000, Inf),
+                         na.value = NA) +
+    scale_x_continuous(breaks = map_layers$lon.breaks) +
+    scale_y_continuous(breaks = map_layers$lat.breaks) +
+    theme_2()
+
+  y_dim_mult <- diff(map_layers$plot.boundary$y)/diff(map_layers$plot.boundary$x)
+
+  ragg::agg_png(filename = here::here("analysis", "efh_bathymetry_2028", "plots", "bathy", paste0("bathy_diff_2022_current_d_", region[ii], ".png")), 
+                width = 6.5, 
+                height = y_dim_mult*6.5, 
+                units = "in", 
+                res = 180)
+  print(plot_disc_bathy_diff)
+  dev.off()
+
+  ragg::agg_png(filename = here::here("analysis", "efh_bathymetry_2028", "plots", "bathy", paste0("bathy_diff_2022_current_c_", region[ii], ".png")), 
+                width = 6.5, 
+                height = y_dim_mult*6.5, 
+                units = "in", 
+                res = 180)
+  print(plot_cont_bathy_diff)
+  dev.off()
+
+
+  # Plot observed tow depth versus raster tow depth
+
+  depth_bins <- c(0, 50, seq(100, 1000, 100), 1200, 1400, Inf)
+
+  towmid <- sf::st_read(here::here("analysis", "efh_bathymetry_2028", "data", "bts_tows", paste0(region[ii], "_towmid.shp"))) |>
+    sf::st_transform(crs = map_layers$crs) |>
+    dplyr::mutate(DEPTH_BIN = cut(DEPTH, 
+                                  breaks = depth_bins,
+                                  labels = paste(depth_bins[1:length(depth_bins)-1], 
+                                                 depth_bins[2:length(depth_bins)], 
+                                                 sep = "-")))
+
+  old_bathy_points <- terra::as.points(old_bathy_reproj) |>
+    sf::st_as_sf()
+
+  new_bathy_points  <- terra::as.points(new_bathy) |>
+    sf::st_as_sf()
+
+  towmid$OLD_BATHY_DEPTH <- old_bathy_points[[1]][sf::st_nearest_feature(towmid, old_bathy_points)]
+  towmid$NEW_BATHY_DEPTH <- new_bathy_points[[1]][sf::st_nearest_feature(towmid, new_bathy_points)]  
+
+  xy_range <- range(c(towmid$OLD_BATHY_DEPTH, towmid$NEW_BATHY_DEPTH, towmid$DEPTH), na.rm = TRUE)
+
+  plot_depth_scatter <- cowplot::plot_grid(
+    ggplot() +
+      geom_abline(intercept = 0, slope = 1, color = "black", linetype = 2) + 
+      geom_point(data = towmid, 
+                 mapping = aes(x = OLD_BATHY_DEPTH, y = DEPTH), color = "grey50", alpha = 0.3) +
+      scale_x_continuous(name = "2022 EFH Depth (m)", 
+                         limits = xy_range) +
+      scale_y_continuous(name = "Observed tow midpoint depth (m)\n2000-2023", 
+                         limits = xy_range) +
+      theme_bw(),
+    ggplot() +
+      geom_abline(intercept = 0, slope = 1, color = "black", linetype = 2) +
+      geom_point(data = towmid, 
+                 mapping = aes(x = NEW_BATHY_DEPTH, y = DEPTH), color = "grey50", alpha = 0.3) +
+      scale_x_continuous(name = "New Compilation Depth (m)", 
+                         limits = xy_range) +
+      scale_y_continuous(name = "Observed tow midpoint depth (m)\n2000-2023", 
+                         limits = xy_range) +
+      theme_bw()
+  )
+
+
+  ragg::agg_png(filename = here::here("analysis", "efh_bathymetry_2028", "plots", "bathy", paste0("bathy_raster_vs_observed_", region[ii], ".png")), 
+                width = 6.5, 
+                height = 3.25, 
+                units = "in", 
+                res = 300)
+  print(plot_depth_scatter)
+  dev.off()
+
+  bias_df <- towmid |>
+    sf::st_drop_geometry() |>
+    dplyr::group_by(DEPTH_BIN) |>
+    dplyr::summarise(NEW_MRE = mean_relative_error(DEPTH, NEW_BATHY_DEPTH),
+                     NEW_MAE = mean_absolute_error(DEPTH, NEW_BATHY_DEPTH),
+                     NEW_BIAS = mean_bias(DEPTH, NEW_BATHY_DEPTH),
+                     OLD_MRE = mean_relative_error(DEPTH, OLD_BATHY_DEPTH),
+                     OLD_MAE = mean_absolute_error(DEPTH, OLD_BATHY_DEPTH),
+                     OLD_BIAS = mean_bias(DEPTH, OLD_BATHY_DEPTH)) |>
+    dplyr::mutate(REGION = toupper(region[ii]))
+
+  write.csv(bias_df, 
+            file = here::here("analysis", "efh_bathymetry_2028", "plots", "bathy", paste0("bathy_raster_accuracy_", region[ii], ".csv")),
+            row.names = FALSE)
+
+}
+
+

analysis/efh_bathymetry_2028/assets/bpi_inner_outer.R → analysis/efh_bathymetry_2028/04_evaluate_bpi.R

@@ -2,6 +2,8 @@ library(akgfmaps)
 library(MultiscaleDTM)
 source(here::here("R", "bbox_to_polygon.R")) # Function to be included in an upcoming akgfmaps release
 
+dir.create(here::here("analysis", "efh_bathymetry_2028", "plots", "BPI"), recursive = TRUE)
+
 locations <- read.csv(file = here::here("analysis", "efh_bathymetry_2028", "data", "terrain_example_zones.csv"))
 
 theme_1 <- function() {
@@ -146,6 +148,7 @@ for(ii in 1:nrow(locations)) {
   ragg::agg_png(filename = here::here("analysis", 
                                       "efh_bathymetry_2028", 
                                       "plots", 
+                                      "BPI",
                                       paste0("BPI_Bathy_", 
                                              gsub(x = set_loc$label, 
                                                   pattern = " ", 
@@ -261,6 +264,7 @@ for(ii in 1:nrow(locations)) {
   ragg::agg_png(filename = here::here("analysis", 
                                       "efh_bathymetry_2028", 
                                       "plots", 
+                                      "BPI",
                                       paste0("BPI_Bathy_", 
                                              gsub(x = set_loc$label, 
                                                   pattern = " ", 
@@ -367,6 +371,7 @@ for(ii in 1:nrow(locations)) {
   ragg::agg_png(filename = here::here("analysis", 
                                       "efh_bathymetry_2028", 
                                       "plots", 
+                                      "BPI",
                                       paste0("BPI_Bathy_", 
                                              gsub(x = set_loc$label, 
                                                   pattern = " ", 

analysis/efh_bathymetry_2028/get_towmid.R

@@ -0,0 +1,106 @@
+library(akgfmaps)
+
+
+# Get tow midpoints
+st_line_midpoints <- function(sf_lines = NULL) {
+
+  g <- sf::st_geometry(sf_lines)
+
+  g_mids <- lapply(g, function(x) {
+
+    coords <- as.matrix(x)
+
+    # this is just a copypaste of View(maptools:::getMidpoints):
+    get_mids <- function (coords) {
+      dist <- sqrt((diff(coords[, 1])^2 + (diff(coords[, 2]))^2))
+      dist_mid <- sum(dist)/2
+      dist_cum <- c(0, cumsum(dist))
+      end_index <- which(dist_cum > dist_mid)[1]
+      start_index <- end_index - 1
+      start <- coords[start_index, ]
+      end <- coords[end_index, ]
+      dist_remaining <- dist_mid - dist_cum[start_index]
+      mid <- start + (end - start) * (dist_remaining/dist[start_index])
+      return(mid)
+    }
+
+    mids <- st_point(get_mids(coords))
+  })
+
+  geometry <- sf::st_sfc(g_mids, 
+                    crs = sf::st_crs(sf_lines))
+
+  out <- sf::st_drop_geometry(sf_lines)
+
+  out <- cbind(out[, which(names(out) != "geometry")], sf::st_sf(geometry))
+
+  return(out)
+}
+
+channel <- navmaps::get_connected(schema = "AFSC")
+
+tow_region <- c("BS", "AI", "GOA")
+label_region <- c("ebs", "ai", "goa")
+
+for(ii in 1:length(tow_region)) {
+
+  haul_data <- RODBC::sqlQuery(
+    channel = channel,
+    query = paste0(
+      "select 
+        vessel, 
+        cruise, 
+        haul, 
+        performance, 
+        start_latitude, 
+        end_latitude, 
+        start_longitude, 
+        end_longitude, 
+        bottom_depth 
+      from 
+        racebase.haul 
+      where 
+        region = '", tow_region[ii], "' 
+        and haul_type = 3 
+        and performance >= 0 
+        and abundance_haul = 'Y' 
+        and bottom_depth > 0
+        and cruise > 200000
+        and gear in (30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 172)")) |>
+      dplyr::filter(!is.na(START_LONGITUDE), !is.na(START_LATITUDE), !is.na(END_LONGITUDE), !is.na(END_LATITUDE))
+
+  towpaths <- dplyr::bind_rows(
+    haul_data |>
+      dplyr::select(-END_LONGITUDE, -END_LATITUDE) |>
+      tidyr::pivot_longer(cols = c(START_LONGITUDE, START_LATITUDE)) |>
+      dplyr::mutate(name = if_else(name == "START_LATITUDE", "LATITUDE", "LONGITUDE")) |>
+      tidyr::pivot_wider(names_from = name, values_from = value),
+    haul_data |>
+      dplyr::select(-START_LONGITUDE, -START_LATITUDE) |>
+      tidyr::pivot_longer(cols = c(END_LONGITUDE, END_LATITUDE)) |>
+      dplyr::mutate(name = if_else(name == "END_LATITUDE", "LATITUDE", "LONGITUDE")) |>
+      tidyr::pivot_wider(names_from = name, values_from = value)
+  ) |>
+    dplyr::filter(!is.na(LATITUDE), !is.na(LONGITUDE)) |>
+    sf::st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = "WGS84") |>
+    dplyr::arrange(VESSEL, CRUISE, HAUL) |>
+    dplyr::group_by(VESSEL, CRUISE, HAUL, BOTTOM_DEPTH) |>
+    dplyr::summarise(do_union = FALSE) |>
+    sf::st_cast(to = "LINESTRING") |>
+    st_line_midpoints() |>
+    dplyr::rename(DEPTH = BOTTOM_DEPTH)
+
+  sf::st_write(obj = towpaths, 
+               dsn = here::here("analysis", 
+                                "efh_bathymetry_2028", 
+                                "data", 
+                                "bts_tows", 
+                                paste0(
+                                  label_region[ii], 
+                                  "_", 
+                                  "towmid.shp")
+               ),
+               append = FALSE
+  )
+
+}

analysis/efh_bathymetry_2028/plots/bathy/bathy_raster_accuracy_ai.csv

@@ -0,0 +1,7 @@
+"DEPTH_BIN","NEW_MRE","NEW_MAE","NEW_BIAS","OLD_MRE","OLD_MAE","OLD_BIAS","REGION"
+"0-50",0.130700631416062,0.0617165337525394,0.927137962302289,0.130679431105905,0.0584301056535382,0.964053418800325,"AI"
+"50-100",0.0504103540567059,0.0213218812837841,1.01837288396307,0.0616163640529554,0.028713241726905,1.00161314549256,"AI"
+"100-200",0.0442685242740121,0.0196529089311664,0.993107304780834,0.0593482797235181,0.028628142827316,0.972992473361545,"AI"
+"200-300",0.291980054602265,0.0285012371404104,1.01650801554026,0.0646337884626028,0.028325323401616,0.99332771656431,"AI"
+"300-400",0.0617915068495647,0.0254249264939416,1.01801980204704,0.0654458703072229,0.0275725016443836,1.00969103037322,"AI"
+"400-500",0.0480878628953349,0.0206544597667234,0.999466069154683,0.0343506161272969,0.0146068409551939,1.01244429056195,"AI"

analysis/efh_bathymetry_2028/plots/bathy/bathy_raster_accuracy_ebs.csv

@@ -0,0 +1,13 @@
+"DEPTH_BIN","NEW_MRE","NEW_MAE","NEW_BIAS","OLD_MRE","OLD_MAE","OLD_BIAS","REGION"
+"0-50",0.0713704382024566,0.0294144718504243,1.05472273490826,0.14356605983855,0.0592144887690071,1.01313434992439,"EBS"
+"50-100",0.0171833483969421,0.0073646550500967,1.00869903791271,0.116303425980895,0.048561734996928,1.01442036366489,"EBS"
+"100-200",0.00780619992351972,0.00339461948474713,0.999754777968775,0.0202338878928625,0.00894136991277215,0.995010057596975,"EBS"
+"200-300",0.0980606699157047,0.067788839525398,0.876937271282848,0.0467842297224619,0.0192811919606618,1.0125405892609,"EBS"
+"300-400",0.08303872161253,0.0445294665257916,0.940298558314596,0.0573571208135955,0.0265270023027767,1.00775811577057,"EBS"
+"400-500",0.0578631037522344,0.027038658396349,0.969287764184541,0.0324880597783217,0.0139822720978878,1.00116353815899,"EBS"
+"500-600",0.0378499961391256,0.0165878283918904,0.988478818921721,0.0252682166425052,0.0109708896511535,0.999356290028403,"EBS"
+"600-700",0.0533014011108076,0.024835399115991,0.97281038624731,0.0286609881351489,0.0123921452388682,0.998888292689543,"EBS"
+"700-800",0.04760386968974,0.0214413467851657,0.99742182495497,0.0361834070761754,0.0151104151087796,1.00674584938589,"EBS"
+"800-900",0.0657494707532494,0.0325930598525242,0.952116381350328,0.0487268984449387,0.0199729852147607,1.0203751536223,"EBS"
+"900-1000",0.0912297562788419,0.0427098163654889,0.97749034724347,0.0345840424224312,0.0138194752556817,1.01239533418843,"EBS"
+"1000-1200",0.0659364930906344,0.0301667321947186,0.972159962774111,0.0412052269760598,0.0174356861470149,1.02058069806395,"EBS"

analysis/efh_bathymetry_2028/plots/bathy/bathy_raster_accuracy_goa.csv

@@ -0,0 +1,12 @@
+"DEPTH_BIN","NEW_MRE","NEW_MAE","NEW_BIAS","OLD_MRE","OLD_MAE","OLD_BIAS","REGION"
+"0-50",0.100315783149378,0.0435178596998945,1.02674295640701,0.136054019283673,0.0630694873621466,1.0015766193335,"GOA"
+"50-100",0.0444696451421345,0.0189582870185426,1.01903394775846,0.0843287758190755,0.0511348444663457,0.951500313422852,"GOA"
+"100-200",0.0327420003294576,0.0142273238041102,1.00515620777472,0.088265700898032,0.0440486320788723,0.942015295667179,"GOA"
+"200-300",0.0317382877924005,0.0136991148250317,1.00271492267274,0.0739745290841734,0.0399614004063242,0.938101775916513,"GOA"
+"300-400",0.0577172300123324,0.0236177135890422,0.994404514490344,0.100452078107725,0.0456358425557637,0.950811896911196,"GOA"
+"400-500",0.0363009517985758,0.0157970924451541,0.997783483397166,0.0410919413676141,0.0181215122625672,0.989718008846834,"GOA"
+"500-600",0.0440328556538761,0.0196317021976391,0.984683623349239,0.064783431522889,0.0290337468767565,0.985675848682744,"GOA"
+"600-700",0.0405348605405812,0.0179934018272694,0.977653846693839,0.075956968827303,0.0357336609629971,0.939251890908243,"GOA"
+"700-800",0.025185222617622,0.011045418373701,0.985790792166459,0.0362929799610986,0.0148952210398448,0.998803672471987,"GOA"
+"800-900",0.0259148888040429,0.0113998630622571,0.989161940939797,0.0354173160486659,0.0155614601882287,0.987802666515987,"GOA"
+"900-1000",0.0156708061464795,0.00672579347083535,1.01557474982531,0.0191293923954764,0.00820160608330911,1.01816599271901,"GOA"