diff --git a/DESCRIPTION b/DESCRIPTION
index e421e4a..dd775c6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: akgfmaps
 Type: Package
 Title: Alaska Groundfish and Ecosystem Survey Area Mapping
-Version: 4.0.0
-Date: 2025-01-03
+Version: 4.0.1
+Date: 2025-01-07
 Authors@R: c(person("Sean", "Rohan", email = "sean.rohan@noaa.gov", role = c("aut", "cre")),
     person("Emily", "Markowitz", email = "emily.markowitz@noaa.gov", role = "ctb"),
     person("Zack", "Oyafuso", email = "zack.oyafuso@noaa.gov", role = "ctb"),
diff --git a/NEWS b/NEWS
index c9309aa..3263f2d 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,12 @@
+akgfmaps 4.0.1 (January 7, 2025)
+----------------------------------------------------------------
+
+BUG FIX
+
+- Change STATION fields to use GRID_X-GRID_Y-STRATUM naming
+  convention.
+
+
 akgfmaps 4.0.0 (January 3, 2025)
 ----------------------------------------------------------------
 
diff --git a/analysis/akgfmaps4/make_bt_gpkg.R b/analysis/akgfmaps4/make_bt_gpkg.R
index 3e68781..8105c71 100644
--- a/analysis/akgfmaps4/make_bt_gpkg.R
+++ b/analysis/akgfmaps4/make_bt_gpkg.R
@@ -81,14 +81,13 @@ goa_stratum_2025 |>
                append = TRUE,
                delete_dsn = FALSE)
 
-
 goa_station_grid_2025 <-
-  sf::st_read(here::here("analysis", "goa_strata_2025", "goaai_grid_2025.shp")) |>
-  sf::st_intersection(goa_stratum_2025) |>
-  dplyr::mutate(AREA_M2 = as.numeric(sf::st_area(geometry)),
-                STATION = paste0(AREA_ID, "-", GRIDID),
+  sf::st_read(here::here("analysis", "goa_strata_2025", "goa_stations_2025.gpkg")) |>
+  dplyr::mutate(AREA_M2 = AREA_KM2*1e6,
+                SURVEY_DEFINITION_ID = 47,
+                DESIGN_YEAR = 2025,
                 AREA_TYPE = "STATION") |>
-  dplyr::select(AREA_TYPE, SURVEY_DEFINITION_ID, DESIGN_YEAR, GRID_ID = id, STATION, AREA_M2)
+  dplyr::select(AREA_TYPE, SURVEY_DEFINITION_ID, DESIGN_YEAR, GRID_ID = GRIDID, STATION, AREA_M2)
 
 goa_station_grid_2025 |>
   sf::st_write(dsn = here::here("inst", "extdata", "afsc_bottom_trawl_surveys.gpkg"),
diff --git a/assets/demo/akgfmaps_SCRUG_20250122/ex1_installing_akgfmaps.R b/assets/demo/akgfmaps_SCRUG_20250122/ex1_installing_akgfmaps.R
new file mode 100644
index 0000000..6e4e269
--- /dev/null
+++ b/assets/demo/akgfmaps_SCRUG_20250122/ex1_installing_akgfmaps.R
@@ -0,0 +1,29 @@
+# 1. Install akgfmaps using the devtools package ---------------------------------------------------
+
+install.packages("devtools", "dplyr")
+devtools::install_github(repo = "afsc-gap-products/akgfmaps")
+
+# 2. Load akgfmaps/verify that you have version 4 installed ------------------------------------
+library(akgfmaps)
+
+# Verify that you have akgfmaps version >= 4.0.0
+(session <- sessionInfo())
+
+session$loadedOnly$akgfmaps$Version > 4
+
+# Troubleshooting ----------------------------------------------------------------------------------
+
+# Problem: API Rate Limit error
+# Solution: Setup a GitHub Personal Access Token in R Studio
+# Instructions: https://carpentries.github.io/sandpaper-docs/github-pat.html
+
+# Workaround without API access
+#
+# 1. Install these packages from CRAN:
+#   ggplot2, sf, stars, terra, gstat (>= 2.0-1), stats, utils, methods, colorspace, RColorBrewer,
+#   here, ggspatial, rmapshaper, shadowtext, ggthemes, classInt (>= 0.4-1), units
+#
+# 2. Download the akgfmaps tarball from:
+#
+#
+# 3. Install the akgmaps tarball using install.packages()
diff --git a/assets/demo/akgfmaps_SCRUG_20250122/ex2_get_layers.R b/assets/demo/akgfmaps_SCRUG_20250122/ex2_get_layers.R
new file mode 100644
index 0000000..a82b313
--- /dev/null
+++ b/assets/demo/akgfmaps_SCRUG_20250122/ex2_get_layers.R
@@ -0,0 +1,346 @@
+# Example: Retrieving data for survey areas
+# SCRUG, akgfmaps, January 22, 2025
+# Created by Sean Rohan <sean.rohan@noaa.gov>
+
+library(akgfmaps)
+
+# 1. get_base_layers() for bottom trawl ------------------------------------------------------------
+# Access AFSC bottom trawl survey regions, handle coordinate reference system (CRS)
+# transformations, create plot breaks and axis labels.
+# Outputs work with any graphical interface that supports 'sf' object classes from the sf package,
+# but is optimized to work with ggplot2.
+
+?get_base_layers
+
+# Get layers for one survey region
+ebs_layers <- akgfmaps::get_base_layers(
+  select.region = "sebs", # EBS shelf bottom traw survey regions
+  set.crs = "EPSG:3338", # Coordinate reference system, in this case Alaska Albers Equal Area
+  high.resolution.coast = FALSE, # Set to TRUE when coast is used for analyses, e.g. raster masking
+  design.year = NULL # Return the most recent/current survey design
+)
+
+# Inspect Layers returned
+names(ebs_layers)
+
+ggplot() +
+  geom_sf(data = ebs_layers$akland) +
+  geom_sf(data = ebs_layers$survey.grid, # Survey station grid
+          fill = NA) +
+  geom_sf(data = ebs_layers$bathymetry) +
+  geom_sf(data = ebs_layers$graticule,
+          alpha = 0.3) +
+  scale_x_continuous(limits = ebs_layers$plot.boundary$x,
+                     breaks = ebs_layers$lon.breaks) +
+  scale_y_continuous(limits = ebs_layers$plot.boundary$y,
+                     breaks = ebs_layers$lat.breaks) +
+  theme_bw() +
+  theme(panel.grid = element_blank())
+
+
+# NEW: Return layers for multiple surveys ----
+all_bt_layers <- akgfmaps::get_base_layers(
+  select.region = c("ebs", "goa", "ai", "ebs.slope", "ecs"), # EBS, EBS slope, NBS, GOA, AI, and ECS
+  set.crs = "EPSG:3338", # Coordinate reference system, in this case Alaska Albers Equal Area
+  high.resolution.coast = FALSE # Set to TRUE when coast is used for analyses, e.g. raster masking
+)
+
+ggplot() +
+  geom_sf(data = all_bt_layers$akland) +
+  geom_sf(data = all_bt_layers$survey.area,
+          mapping = aes(fill = SURVEY_NAME)) +
+  geom_sf(data = all_bt_layers$bathymetry) +
+  geom_sf(data = all_bt_layers$graticule,
+          alpha = 0.3) +
+  scale_x_continuous(limits = all_bt_layers$plot.boundary$x,
+                     breaks = all_bt_layers$lon.breaks) +
+  scale_y_continuous(limits = all_bt_layers$plot.boundary$y,
+                     breaks = all_bt_layers$lat.breaks) +
+  scale_fill_viridis_d(name = "Survey Name") +
+  theme_bw() +
+  theme(panel.grid = element_blank())
+
+
+# NEW: All bottom trawl survey regions have the same fields. Fields match GAP_PRODUCTS. ----
+head(all_bt_layers$survey.area)
+head(all_bt_layers$survey.grid)
+head(all_bt_layers$survey.strata) # STRATUM field matches AREA_ID field in GAP_PRODUCTS and STRATUM
+# field in RACEBASE
+
+# NEW: Retrieve layers for a specific design year for a bottom trawl survey ----
+#   - Returns the most recent design year when NULL.
+#   - Numerous caveats may apply when dealing with old survey designs (e.g., area calculation issues)
+#   - Many historical designs are not available and may never be.
+
+goa_1984_layers <- akgfmaps::get_base_layers(
+  select.region = "goa.west", # Setup plot for WGOA
+  design.year = 1984,
+  set.crs = "EPSG:3338"
+)
+
+goa_2025_layers <- akgfmaps::get_base_layers(
+  select.region = "goa.west",
+  design.year = 2025,
+  set.crs = "EPSG:3338"
+)
+
+# Combine survey area strata
+goa_1984_2025_strata <- rbind(
+  goa_1984_layers$survey.strata,
+  goa_2025_layers$survey.strata
+)
+
+# Plot 1984 and 2025 survey strata for the WGOA
+ggplot() +
+  geom_sf(data = goa_1984$akland) +
+  geom_sf(data = goa_1984_2025_strata,
+          mapping = aes(color = factor(DESIGN_YEAR)),
+          fill = NA,
+          alpha = 0.5) +
+  scale_color_manual(name = "Design year",
+                     values = c("red", "black")) +
+  scale_x_continuous(limits = goa_1984$plot.boundary$x,
+                     breaks = goa_1984$lon.breaks) +
+  scale_y_continuous(limits = goa_1984$plot.boundary$y,
+                     breaks = goa_1984$lat.breaks) +
+  theme_bw()
+
+goa_1984_2025_grid <- rbind(
+  goa_1984_layers$survey.grid,
+  goa_2025_layers$survey.grid
+)
+
+# Plot 1984 and 2025 survey grid for the WGOA
+ggplot() +
+  geom_sf(data = goa_1984$akland) +
+  geom_sf(data = goa_1984_2025_grid,
+          mapping = aes(color = factor(DESIGN_YEAR)),
+          fill = NA,
+          alpha = 0.5) +
+  scale_color_manual(name = "Design year",
+                     values = c("red", "black")) +
+  scale_x_continuous(limits = goa_1984$plot.boundary$x,
+                     breaks = goa_1984$lon.breaks) +
+  scale_y_continuous(limits = goa_1984$plot.boundary$y,
+                     breaks = goa_1984$lat.breaks) +
+  theme_bw()
+
+# EBS shelf survey with corner stations for the 2023 survey
+# The EBS survey dropped corner stations for the 2024 survey. Using design.year = 2023 (last year
+# with corner stations) returns the latest available design year for the survey grid and strata.
+#
+# Note that there are warning when selecting a design.year that does not exist; the closest
+# preceding year is returned instead.
+
+ebs_2023 <- akgfmaps::get_base_layers(
+  select.region = "sebs",
+  design.year = 2023,
+  set.crs = "EPSG:3338"
+)
+
+# DESIGN_YEAR for survey.grid is 2010, the year the grid was redefined and extended to the NBS
+head(ebs_2023$survey.grid)
+
+# DESIGN_YEAR for survey.strata is 2022, the year EBS shelf and NBS survey were snapped together and
+# areas were recalculated using the Albers equal area projection
+head(ebs_2023$survey.strata)
+
+
+# Error when trying to retrieve layers for a design.year that's earlier than the earliest available
+# in the package. It's likely that some of these years cannot be recovered or may not have existed
+# in digital format.
+akgfmaps::get_base_layers(
+  select.region = "sebs",
+  design.year = 1982,
+  set.crs = "EPSG:3338"
+)
+
+ggplot() +
+  geom_sf(data = ebs_2023$survey.grid)
+
+# Aleutian Islands bottom trawl survey removed Bowers Ridge in 1991 (start of the standardized
+# survey time series)
+
+ai_1980 <- get_base_layers(
+  select.region = "ai",
+  design.year = 1980,
+  set.crs = "EPSG:3338"
+)
+
+ai_1991 <- get_base_layers(
+  select.region = "ai",
+  design.year = 1991,
+  set.crs = "EPSG:3338"
+)
+
+ai_1980_1991_strata <- rbind(ai_1980$survey.strata, ai_1991$survey.strata)
+
+ggplot() +
+  geom_sf(data = ai_1980_1991_strata,
+          mapping = aes(fill = factor(STRATUM)),
+          color = NA) +
+  facet_wrap(~paste0("DESIGN_YEAR = ", DESIGN_YEAR)) +
+  scale_fill_viridis_d(option = "H") +
+  scale_x_continuous(limits = ai_1980$plot.boundary$x,
+                     breaks = ai_1980$lon.breaks) +
+  scale_y_continuous(limits = ai_1980$plot.boundary$y,
+                     breaks = ai_1980$lat.breaks) +
+  theme_bw() +
+  theme(legend.position = "none")
+
+
+# 2. NEW: get_base_layers() for longline survey ----------------------------------------------------
+# Access AFSC longline survey layers
+# - Currently no access to historical designs using design.year
+
+ll_goa_layers <- akgfmaps::get_base_layers(
+  select.region = "ll.goa",
+  set.crs = "EPSG:3338",
+  design.year = NULL
+)
+
+ggplot() +
+  geom_sf(data = ll_goa_layers$akland) + # Same land polygon
+  geom_sf(data = ll_goa_layers$survey.grid) + # Survey stations rather than a grid
+  geom_sf_text(data = ll_goa_layers$survey.grid,
+               mapping = aes(label = STATION_NUMBER),
+               hjust = 0.5,
+               vjust = -0.5) +
+  scale_x_continuous(limits = ll_goa_layers$plot.boundary$x,
+                     breaks = ll_goa_layers$lon.breaks) +
+  scale_y_continuous(limits = ll_goa_layers$plot.boundary$y,
+                     breaks = ll_goa_layers$lat.breaks) +
+  theme_bw() +
+  theme(axis.title = element_blank())
+
+# Example of cleaner station labels using ggrepel
+
+install.packages("ggrepel")
+
+library(ggrepel)
+
+ll_goa_layers$survey.grid[, c('x', 'y')] <- sf::st_coordinates(ll_goa_layers$survey.grid)
+
+ggplot() +
+  geom_sf(data = ll_goa_layers$akland) + # Same land polygon
+  geom_sf(data = ll_goa_layers$survey.grid) + # Survey stations rather than a grid
+  geom_text_repel(data = ll_goa_layers$survey.grid,
+                  mapping = aes(x = x,
+                                y = y,
+                                label = STATION_NUMBER),
+                  hjust = 0.5,
+                  vjust = -0.5) +
+  scale_x_continuous(limits = ll_goa_layers$plot.boundary$x,
+                     breaks = ll_goa_layers$lon.breaks) +
+  scale_y_continuous(limits = ll_goa_layers$plot.boundary$y,
+                     breaks = ll_goa_layers$lat.breaks) +
+  theme_bw() +
+  theme(axis.title = element_blank())
+
+
+# 3. NEW: EBS crab stratum boundaries --------------------------------------------------------------
+# Access stratum boundaries for EBS and NBS crab stocks
+# - Select by stock or return all crab
+# - Does not return additional objects, e.g.,
+
+?get_crab_strata
+
+# Layers for a single stock, Bristol Bay red king crab
+pirkc_strata <- akgfmaps::get_crab_strata(
+  select.stock = "pirkc",
+  select.region = "ebs",
+  set.crs = "EPSG:3338")
+
+# Layers for all EBS stocks
+all_ebs_crab <- akgfmaps::get_crab_strata(
+  select.stock = NULL,
+  select.region = "ebs",
+  set.crs = "EPSG:3338")
+
+head(all_ebs_crab)
+
+ggplot() +
+  geom_sf(data = ebs_layers$akland) +
+  geom_sf(data = pirkc_strata,
+          mapping = aes(fill = STRATUM)) +
+  scale_fill_viridis_d() +
+  facet_wrap(~STOCK) +
+  theme_bw()
+
+# 4. NMFS, ADFG, ESR, and BSIERP layers ------------------------------------------------------------
+# - These layers are also available in Matt Callahan's akmarineareas2 package
+#   (https://github.com/MattCallahan-NOAA/akmarineareas2)
+
+install.packages("shadowtext")
+library(shadowtext)
+
+# Alaska Department of Fish and Game (ADFG) Areas
+adfg_areas <- get_adfg_areas(
+  set.crs = "EPSG:3338",
+  subset.fields = TRUE # Only return a subset of columns
+  )
+
+head(adfg_areas)
+
+ggplot() +
+  geom_sf(data = adfg_areas)
+
+# National Marine Fisheries Service (NMFS) Statistical Areas
+nmfs_areas <- get_nmfs_areas(set.crs = "EPSG:3338")
+
+nmfs_area_labels <- sf::st_centroid(nmfs_areas)
+nmfs_area_labels[, c('x', 'y')] <- sf::st_coordinates(nmfs_area_labels)
+
+ggplot() +
+  geom_sf(data = nmfs_areas, color = "grey40") +
+  geom_shadowtext(data = nmfs_area_labels,
+          mapping = aes(x = x, y = y, label = REP_AREA),
+          bg.color = "white",
+          color = "black") +
+  theme_bw()
+
+# Ecosystem Status Report (ESR) regions
+esr_areas <- get_esr_regions(
+  select.region = "esr_area",
+  set.crs = "EPSG:3338"
+)
+
+ggplot() +
+  geom_sf(data = esr_areas,
+          mapping = aes(fill = AREA_NAME)) +
+  scale_fill_viridis_d(option = "A") +
+  theme_bw()
+
+esr_subareas <- get_esr_regions(
+  select.region = "esr_subarea",
+  set.crs = "EPSG:3338"
+)
+
+ggplot() +
+  geom_sf(data = esr_subareas,
+          mapping = aes(fill = AREA_NAME)) +
+  scale_fill_viridis_d(option = "B") +
+  theme_bw()
+
+esr_subareas_inside <- get_esr_regions(
+  select.region = "esr_subarea_inside",
+  set.crs = "EPSG:3338"
+)
+
+ggplot() +
+  geom_sf(data = esr_subareas_inside,
+          mapping = aes(fill = AREA_NAME)) +
+  scale_fill_viridis_d(option = "E") +
+  theme_bw()
+
+# Bering Sea Integrated Ecosystem Research Program (BSIERP) regions
+bsierp_areas <- get_bsierp_regions(set.crs = "EPSG:3338")
+
+head(bsierp_areas)
+
+ggplot() +
+  geom_sf(data = bsierp_areas) +
+  geom_sf_text(data = sf::st_centroid(bsierp_areas),
+          mapping = aes(label = BSIERP_ID)) +
+  theme_bw() +
+  theme(axis.title = element_blank())
diff --git a/assets/demo/akgfmaps_SCRUG_20250122/ex3_make_idw_map.R b/assets/demo/akgfmaps_SCRUG_20250122/ex3_make_idw_map.R
new file mode 100644
index 0000000..7d58f71
--- /dev/null
+++ b/assets/demo/akgfmaps_SCRUG_20250122/ex3_make_idw_map.R
@@ -0,0 +1,3 @@
+# Example: Access data directly from a geopackage
+# SCRUG, akgfmaps, January 22, 2025
+# Created by Sean Rohan (GitHub: sean-rohan-noaa)
diff --git a/assets/demo/akgfmaps_SCRUG_20250122/geopackages.R b/assets/demo/akgfmaps_SCRUG_20250122/ex4_geopackages.R
similarity index 88%
rename from assets/demo/akgfmaps_SCRUG_20250122/geopackages.R
rename to assets/demo/akgfmaps_SCRUG_20250122/ex4_geopackages.R
index b3ae0aa..5e32252 100644
--- a/assets/demo/akgfmaps_SCRUG_20250122/geopackages.R
+++ b/assets/demo/akgfmaps_SCRUG_20250122/ex4_geopackages.R
@@ -1,6 +1,12 @@
 # Example: Access data directly from a geopackage
 # SCRUG, akgfmaps, January 22, 2025
-# Created by Sean Rohan <sean.rohan@noaa.gov>
+# Created by Sean Rohan (GitHub: sean-rohan-noaa)
+
+# akgfmaps installs layer files locally
+# Layers used to include a variety of formats; starting with akgfmaps v4, layers will be included in
+# geopackage (.gpkg) files.
+list.files(system.file("extdata", package = "akgfmaps"),
+           full.names = TRUE)
 
 # Loading layers directly from built-in Geopackage
 # Geopackages are an open, standards-based, platform-independent format for storing geospatial data.
diff --git a/assets/demo/akgfmaps_SCRUG_20250122/ex5_get_layers_old_version.R b/assets/demo/akgfmaps_SCRUG_20250122/ex5_get_layers_old_version.R
new file mode 100644
index 0000000..6e24903
--- /dev/null
+++ b/assets/demo/akgfmaps_SCRUG_20250122/ex5_get_layers_old_version.R
@@ -0,0 +1,75 @@
+# Example: Using the old get_base_layers() function
+# SCRUG, akgfmaps, January 22, 2025
+# Created by Sean Rohan (GitHub: sean-rohan-noaa)
+
+# Through the end of 2025, an option will be provided to use the old interface and layers using
+# get_base_layers_v3(). The function is limited to the functionality and layers that were avaialble
+# in akgfmaps v 3.6.2
+#
+# The old version has the following limitations:
+# - No access to design years (e.g., the most recent GOA BT survey stratum shapefile will be 1984)
+# - Cannot select multiple regions at the same time
+# - Field names differ among shapefiles
+# - No correspondence with GAP_PRODUCTS
+# - Access to longline
+
+# Recommendation: Don't use this version unless you can't replace with the new version due to time
+# and reliance on fields in specific shapefiles.
+
+library(akgfmaps)
+
+# Arguments differ between function versions
+?get_base_layers
+?get_base_layers_v3
+
+# Retrieve and plot GOA and EBS/NBS strata
+goa_1984_v3 <- akgfmaps::get_base_layers_v3(
+  select.region = "goa",
+  set.crs = "EPSG:3338",
+  high.resolution.coast = FALSE
+)
+
+ebs_nbs_2024_v3 <- akgfmaps::get_base_layers_v3(
+  select.region = "ebs",
+  set.crs = "EPSG:3338",
+  high.resolution.coast = FALSE
+)
+
+ggplot() +
+  geom_sf(data = goa_1984_v3$akland) +
+  geom_sf(data = goa_1984_v3$survey.strata,
+          mapping = aes(
+            fill = factor(STRATUM)
+          )
+  ) +
+  scale_x_continuous(limits = goa_1984_v3$plot.boundary$x,
+                     breaks = goa_1984_v3$lon.breaks) +
+  scale_y_continuous(limits = goa_1984_v3$plot.boundary$y,
+                     breaks = goa_1984_v3$lat.breaks) +
+  theme_bw() +
+  scale_fill_viridis_d(option = "H") +
+  theme(legend.position = "none")
+
+ggplot() +
+  geom_sf(data = ebs_nbs_2024_v3$akland) +
+  geom_sf(data = ebs_nbs_2024_v3$survey.strata,
+          mapping = aes(
+            fill = factor(Stratum)
+          )
+  ) +
+  scale_x_continuous(limits = ebs_nbs_2024_v3$plot.boundary$x,
+                     breaks = ebs_nbs_2024_v3$lon.breaks) +
+  scale_y_continuous(limits = ebs_nbs_2024_v3$plot.boundary$y,
+                     breaks = ebs_nbs_2024_v3$lat.breaks) +
+  scale_fill_viridis_d(option = "H") +
+  theme_bw() +
+  theme(legend.position = "none")
+
+# Note that the field names do not match within or between regions.
+head(goa_1984_v3$survey.area)
+head(goa_1984_v3$survey.strata)
+head(goa_1984_v3$survey.grid)
+
+head(ebs_nbs_2024_v3$survey.area)
+head(ebs_nbs_2024_v3$survey.strata)
+head(ebs_nbs_2024_v3$survey.grid)
diff --git a/assets/demo/akgfmaps_SCRUG_20250122/installing_akgfmaps.R b/assets/demo/akgfmaps_SCRUG_20250122/installing_akgfmaps.R
deleted file mode 100644
index 56327a8..0000000
--- a/assets/demo/akgfmaps_SCRUG_20250122/installing_akgfmaps.R
+++ /dev/null
@@ -1,4 +0,0 @@
-# Installing akgfmaps
-
-install.packages("devtools", "dplyr")
-devtools::install_github(repo = "afsc-gap-products/akgfmaps")
diff --git a/demo/eval_plot_breaks.R b/demo/eval_plot_breaks.R
deleted file mode 100644
index a650d78..0000000
--- a/demo/eval_plot_breaks.R
+++ /dev/null
@@ -1,84 +0,0 @@
-# Contributed by Sean Rohan on November 24, 2023
-
-# eval_plot_breaks helps users select breaks for make_idw_map() and 
-# make_idw_stack(). However, users are encouraged to select breaks 
-# based on the properties of their data. This demo shows how to use 
-# eval_plot_breaks() for exploring options for breaks.
-
-library(akgfmaps)
-library(tidyr)
-
-# Example 1: Breaks with samples drawn from a normal distribution
-dat <- data.frame(x = abs(rnorm(450, mean = 100, sd = 50)))
-
-set.seed(1251)
-
-algo_breaks <- akgfmaps::eval_plot_breaks(CPUE = dat$x, 
-                                          n.breaks = 4)
-
-# Breaks are spread out pretty evenly among bins. If these were spatial data.
-# If these were spatial data, it's pretty likely jenks would work pretty well and we would see 
-# clear distributional patterns. 
-table(cut(dat$x, 
-          breaks = c(0, algo_breaks[algo_breaks$style == "jenks" , 2:6]),
-          right = TRUE))
-
-brks <- algo_breaks |>
-  tidyr::pivot_longer(cols = 2:6)
-
-ggplot() +
-  geom_density(data = dat,
-               mapping = aes(x = x)) +
-  geom_vline(data = brks,
-             mapping = aes(xintercept = value, color = name)) +
-  scale_x_log10() +
-  facet_wrap(~style)
-
-ggplot() +
-  stat_ecdf(data = dat,
-            mapping = aes(x = x)) +
-  geom_vline(data = brks,
-             mapping = aes(xintercept = value, color = name)) +
-  scale_x_log10() +
-  facet_wrap(~style)
-
-# Example 2: Breaks with a multi-modal distribution where modes cover multiple orders of magnitude.
-
-set.seed(1251)
-
-dat <- data.frame(x = abs(c(rnorm(200, mean = 2, sd = 0.4),
-                            rnorm(200, mean = 5, sd = 1),
-                            rnorm(30, mean = 10, sd = 2),
-                            rnorm(17, mean = 100, sd = 20),
-                            rnorm(3, mean = 1e3, sd = 2e4))))
-
-algo_breaks <- akgfmaps::eval_plot_breaks(CPUE = dat$x, 
-                                          n.breaks = 4)
-
-# Jenks breaks would result in >99% of data falling in one bin (1-163)
-# If these were spatial data, the whole map would look the same except for a couple of hot spots.
-table(cut(dat$x, 
-          breaks = c(0, algo_breaks[algo_breaks$style == "jenks" , 2:6]),
-          right = FALSE))
-
-brks <- algo_breaks |>
-  tidyr::pivot_longer(cols = 2:6)
-
-# It's unclear which would be the best automatic option based on the distribution of the data, 
-# although all of them perform pretty poorly. Manual break selection would be needed to provide
-# a reasonsable representation of species distributions.
-ggplot() +
-  geom_density(data = dat,
-               mapping = aes(x = x)) +
-  geom_vline(data = brks,
-             mapping = aes(xintercept = value, color = name)) +
-  scale_x_log10() +
-  facet_wrap(~style)
-
-ggplot() +
-  stat_ecdf(data = dat,
-            mapping = aes(x = x)) +
-  geom_vline(data = brks,
-             mapping = aes(xintercept = value, color = name)) +
-  scale_x_log10() +
-  facet_wrap(~style)
\ No newline at end of file
diff --git a/demo/get_base_layers.R b/demo/get_base_layers.R
deleted file mode 100644
index a8b9f01..0000000
--- a/demo/get_base_layers.R
+++ /dev/null
@@ -1,36 +0,0 @@
-# Created by Sean Rohan on November 24, 2023
-
-library(akgfmaps)
-
-# Get layers for the EBS survey area with Alaska Albers Equal area projection (EPSG:3338). Refer to function documentation for valid survey regions.
-map_layers <- akgfmaps::get_base_layers(select.region = "bs.south", 
-                                  set.crs = "EPSG:3338")
-
-(map_layers_plot <- ggplot() +
-   geom_sf(data = map_layers$akland) +
-   geom_sf(data = map_layers$bathymetry) +
-   geom_sf(data = map_layers$survey.area, fill = NA) +
-   geom_sf(data = map_layers$graticule, color = "grey70", alpha = 0.5) +
-   coord_sf(xlim = map_layers$plot.boundary$x, 
-            ylim = map_layers$plot.boundary$y) +
-   scale_x_continuous(name = "Longitude", 
-                      breaks = map_layers$lon.breaks) + 
-   scale_y_continuous(name = "Latitude", 
-                      breaks = map_layers$lat.breaks) + 
-   theme_bw())
-
-
-## Add a label at 58 N, -165 W by specifying coordinates and converting to
-# Alaska AEA projection
-my_label <- data.frame(x = -165, 
-                       y = 58, 
-                       label = "Thar be\ndragons") |> 
-  akgfmaps::transform_data_frame_crs(in.crs = "+proj=longlat", 
-                                     out.crs = map_layers$crs)
-
-map_layers_plot + 
-  geom_text(data = my_label, 
-            mapping = aes(x = x, 
-                          y = y, 
-                          label = label), 
-            color = "red")
\ No newline at end of file
diff --git a/inst/extdata/afsc_bottom_trawl_surveys.gpkg b/inst/extdata/afsc_bottom_trawl_surveys.gpkg
index 3272bbf..e2a6753 100644
Binary files a/inst/extdata/afsc_bottom_trawl_surveys.gpkg and b/inst/extdata/afsc_bottom_trawl_surveys.gpkg differ