afsc-gap-products / akgfmaps

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

facet_warp/_grid-friendly idw grid outputs #28

Closed EmilyMarkowitz-NOAA closed 1 year ago

EmilyMarkowitz-NOAA commented 1 year ago

The make_idw_map() function is great for producing one map's worth of map elements and plot. However, I am often tasked with plotting multiple maps in a multi-panel plot (by year example below). To do this, I currently run the make_idw_map() function for a subset of data, and extract the extrapolation.grid and continuous.grid from each run's output so I can plot of the all grids using facet_wrap( ~ year).

It may be helpful to develop a function that produces the IDW for each data subset outside of make_idw_map() (essentially replacing this section of the function code). Additionally, it could be useful to allow the user to define a subset or group variable in this new function that (though it would typically be year) could allow the user to systematically iterate through different and multiple criteria (e.g., species, length class bin).

This may be overexplaining, a proposed workflow may look something like this (not run):

map_layers <- akgfmaps::get_base_layers(select.region = region, set.crs = out.crs)

# 1 group -----------------------------------------------------
extrap.grid <- akgfmaps::make_idw_df(..., # new function
                                       group = c("year")) 
p1 <- ggplot2::ggplot() +
          ggplot2::geom_sf(data = extrap.grid,
                                       mapping = aes(fill = var1.pred)) %>%
          ggplot2::facet_wrap( ~ group) +
          ggplot2::geom_sf(data = map_layers$survey.area, fill = NA)

# 2 groups -----------------------------------------------------
extrap.grid <- akgfmaps::make_idw_df(..., # new function
                                       group = c("year", "length_class")) 
p2 <- ggplot2::ggplot() +
          ggplot2::geom_sf(data = extrap.grid,
                                       mapping = aes(fill = var1.pred)) %>%
          ggplot2::facet_wrap(length_class ~ group) +
          ggplot2::geom_sf(data = map_layers$survey.area, fill = NA)

fig-dist-fish-walleye-pollock-below

sean-rohan-NOAA commented 1 year ago

@EmilyMarkowitz-NOAA I added a function to the slope branch called make_idw_stack(). The changes were not as simple as modifying the code in the section you highlighted in make_idw_map() because I needed to maintain the automatic breaks functionality (causes more issues than using manual breaks).

The new function only returns a subset of the objects returned by make_idw_map becuase map_layers, extrapolation.stack, continuous.stack, region, and crs, but can return both sf and stars objects.

I've included an example of how to run the function below. Can you check that it does what you wanted for this feature? I still need to test the function, clean the code up a bit, and potentially add some additional documentation before merging w/ main.

library(akgfmaps)

stars_stack <- make_idw_stack(
  x = akgfmaps::YFS2017 |>
    dplyr::mutate(grp = plyr::round_any(HAUL, 100)),
  COMMON_NAME = NA,
  LATITUDE = NA,
  LONGITUDE = NA,
  CPUE_KGHA = NA,
  region = "bs.south",
  extrap.box = NULL,
  extrapolation.grid.type = "stars",
  set.breaks = "jenks",
  grouping.vars = c("grp", "YEAR"),
  grid.cell = c(5000,5000),
  in.crs = "+proj=longlat",
  out.crs = "EPSG:3338",
  log.transform = FALSE,
  idw.nmax = 4,
  use.survey.bathymetry = TRUE)

plot(stars_stack$continuous.stack[1])
plot(stars_stack$extrapolation.stack[1])
plot(stars_stack$continuous.stack[2])
plot(stars_stack$extrapolation.stack[2])
plot(stars_stack$continuous.stack[3])
plot(stars_stack$extrapolation.stack[3])

sf_stack <- make_idw_stack(
  x = akgfmaps::YFS2017 |>
    dplyr::mutate(grp = plyr::round_any(HAUL, 100)),
  COMMON_NAME = NA,
  LATITUDE = NA,
  LONGITUDE = NA,
  CPUE_KGHA = NA,
  region = "bs.south",
  extrap.box = NULL,
  extrapolation.grid.type = "sf",
  set.breaks = "jenks",
  grouping.vars = c("grp", "YEAR"),
  grid.cell = c(5000,5000),
  in.crs = "+proj=longlat",
  out.crs = "EPSG:3338",
  log.transform = FALSE,
  idw.nmax = 4,
  use.survey.bathymetry = TRUE)

ggplot() +
  geom_sf(data = sf_stack$extrapolation.stack, 
          mapping = aes(fill = var1.pred)) +
  facet_wrap(~grp)

plot(sf_stack$continuous.stack[1])

sf_simple_stack <- make_idw_stack(
  x = akgfmaps::YFS2017 |>
    dplyr::mutate(grp = plyr::round_any(HAUL, 100)),
  COMMON_NAME = NA,
  LATITUDE = NA,
  LONGITUDE = NA,
  CPUE_KGHA = NA,
  region = "bs.south",
  extrap.box = NULL,
  extrapolation.grid.type = "sf.simple",
  set.breaks = "jenks",
  grouping.vars = c("grp", "YEAR"),
  grid.cell = c(5000,5000),
  in.crs = "+proj=longlat",
  out.crs = "EPSG:3338",
  log.transform = FALSE,
  idw.nmax = 4,
  use.survey.bathymetry = TRUE)

ggplot() +
  geom_sf(data = sf_simple_stack$extrapolation.stack, 
          mapping = aes(fill = var1.pred)) +
  facet_wrap(~grp)

plot(sf_simple_stack$continuous.stack[1])
sean-rohan-NOAA commented 1 year ago

Closed by #31