afsc-gap-products / akgfmaps

Make AFSC bottom trawl survey maps and retrieve map layers
Other
17 stars 6 forks source link

Add NMFS and INPFC areas #7

Closed sean-rohan-NOAA closed 1 year ago

JaneSullivan-NOAA commented 1 year ago

Hi Sean,

First, I love this package. I find mapping projections and shapefiles really challenging, and I appreciate how straightforward you've made it here.

I used the source code recently to make a map of NMFS areas. I thought I'd leave an issue to potentially add this functionality, but looks like you already have one open. I don't have INPFC areas unfortunately but could help track it down if you want.

image


set.crs <- 'EPSG:3338'

akland <- sf::st_read(system.file("extdata", "alaska_canada_dcw.shp", package = "akgfmaps"), quiet = TRUE)
# akland <- sf::st_read("data/P3_AK_All.shp", quiet = TRUE)
nmfsareas <- sf::st_read("data/NMFS_Zones_Clean.shp", quiet = TRUE)
nmfscentroids <- sf::st_read("data/NMFS_Area_Centroid.shp", quiet = TRUE)

lon.breaks <- c(170, 175, seq(-180, -130, 5))
lat.breaks <- seq(42, 64, 2)

plot.boundary <- sf::st_bbox(nmfsareas)
plot.boundary <- data.frame(x = c(plot.boundary['xmin'], plot.boundary['xmax']),
                            y = c(plot.boundary['ymin'], plot.boundary['ymax']))

# Set CRS
if(tolower(class(set.crs)) != "crs") {
  set.crs <- sf::st_crs(set.crs)
}

# Make graticule
graticule <- sf::st_graticule(lat = lat.breaks, 
                              lon = lon.breaks, 
                              margin = 1e-5)

# Set CRS for layers 
akland <- akland %>% sf::st_transform(crs = set.crs)
nmfsareas <- nmfsareas %>% sf::st_transform(crs = set.crs)
nmfscentroids <- nmfscentroids %>% sf::st_transform(crs = set.crs)

baselay <- list(akland = akland,
                nmfsareas = nmfsareas, 
                nmfscentroids = nmfscentroids,
                crs = set.crs,
                plot.boundary = plot.boundary,
                lon.breaks = lon.breaks,
                lat.breaks = lat.breaks)

# clean up the labels...
baselay$nmfscentroids <- baselay$nmfscentroids %>%
  mutate(YCoord = ifelse(REP_AREA == 519, 
                         YCoord + 5e4,
                         ifelse(REP_AREA %in% c(514, 541, 542, 543),
                                YCoord - 1.8e5,
                                ifelse(REP_AREA %in% c(620, 630),
                                       YCoord - 8e4,
                                       ifelse(REP_AREA c(550),
                                              YCoord - 5e4,
                                              ifelse(REP_AREA %in%  c(659),
                                                     YCoord + 1.2e5,
                                                     ifelse(REP_AREA %in%  c(649),
                                                            YCoord + 1e5,
                                                     YCoord)))))),
         XCoord = ifelse(REP_AREA == 659,
                         XCoord + 5e4, XCoord))
ak_plot <- ggplot() +
  geom_sf(data = baselay$nmfsareas, fill = NA, col = 'grey70') +
  geom_sf(data = baselay$graticule, color = "black", alpha = 0.5) +
  geom_sf(data = baselay$akland) +
  coord_sf(xlim = baselay$plot.boundary$x, 
           ylim = baselay$plot.boundary$y) +
  scale_x_continuous(name = "Longitude", 
                     breaks = baselay$lon.breaks) + 
  scale_y_continuous(name = "Latitude", 
                     breaks = baselay$lat.breaks) + 
  geom_text(data = baselay$nmfscentroids, aes(x = XCoord, y = YCoord, label = REP_AREA),
            check_overlap = TRUE, size=3, fontface = "bold") +
  theme_bw() 
ak_plot
sean-rohan-NOAA commented 1 year ago

@JaneSullivan-NOAA Glad you're finding the package useful and thank you for chiming in about NMFS areas. I'll aim to include NMFS areas in the next update, which should happen within the next month.

I could likely add an INPFC shapefile at the same time and would greatly appreciate it if you have a suggestion for a shapefile. I've seen one floating around but I'm not entirely sure where it came from.

sean-rohan-NOAA commented 1 year ago

I think this should now be addressed by #33 (version 3.1.0), but please let me know if you notice any issues.

A few code examples:

INPFC

# Retrieve GOA INPFC strata and set the crs to Alaska Albers Equal Area (functions do the same thing)
get_base_layers(select.region = "inpfc.goa", set.crs = "EPSG:3338)
get_inpfc_strata(select.region = "goa", set.crs = "EPSG:3338)

# AI INPFC strata with the default CRS (functions do the same thing)
get_base_layers(select.region = "inpfc.ai")
get_inpfc_strata(select.region = "ai")

NMFS

# Functions do the same thing
get_base_layers(select.region = "nmfs", set.crs = "EPSG:3338)
get_nmfs_areas(set.crs = "EPSG:3338)