diff --git a/.gitignore b/.gitignore index d1978e7..671b963 100644 --- a/.gitignore +++ b/.gitignore @@ -38,9 +38,11 @@ inst/doc doc Meta -*.png *.lock /vignettes/*.html /vignettes/*.R -/vignettes/*.md \ No newline at end of file +/vignettes/*.md +/vignettes/*.png +/vignettes/figure +/plots/ \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index aad4fcd..be1f7e5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: akgfmaps Type: Package Title: Alaska Groundfish and Ecosystem Survey Area Mapping -Version: 3.2.0 +Version: 3.3.0 Authors@R: c(person("Sean", "Rohan", email = "sean.rohan@noaa.gov", role = c("aut", "cre")), person("Jason", "Conner", email = "jason.conner@noaa.gov", role = "ctb"), person("Liz", "Dawson", email = "liz.dawson@noaa.gov", role = "ctb"), diff --git a/NEWS b/NEWS index 0ba70ba..7d4850d 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,58 @@ +akgfmaps 3.3.0 (October 19, 2023) +---------------------------------------------------------------- + +NEW FEATURE + +- Added an include_grid_cell argument to make_2d_grid() that + makes the function return columns with the center coordinates + (lon_plot, lat_plot) of each grid cell prior to intersecting + with the user-provided sf object (passed to obj argument). + Centroids for intersected grid polygons are still returned by + the function for use in spatial analyses and index production + (e.g. VAST, sdmTMB) because they are the center of grid cells + after intersecting with stratum shapefiles. The centroid prior + to intersection is more useful for plotting gridded data with- + out prior transformation. + + EXAMPLE + + Grid cell centroids before intersection (regularly-spaced): + + +-----------+-----------+ + | | | + | | | + | x | x | + | | | + | | | + +-----------+-----------+ + + Grid cell centroids after intersection (irregularly-spaced; + useful for analyses that require accurate distances): + + +-----------+-----------+ + | /| | + | x / | | + | / | x | + | / x | | + | / | | + +-----------+-----------+ + + Grid cell centroids for plotting after intersection (still + regularly-spaced): + + +-----------+-----------+ + | /| | + | / | | + | X / | x | + | / | | + | / | | + +-----------+-----------+ + + USAGE + + See https://github.com/afsc-gap-products/akgfmaps/blob/main/assets/make_plot_grid.md + + akgfmaps 3.2.0 (October 6, 2023) ---------------------------------------------------------------- diff --git a/R/make_2d_grid.R b/R/make_2d_grid.R index 349b102..d563ca9 100644 --- a/R/make_2d_grid.R +++ b/R/make_2d_grid.R @@ -7,9 +7,10 @@ #' @param output_type Type of output as a 1L character vector ("point", "raster", "polygon). #' @param bbox Optional. A 4L numeric vector to define the edges of the interpolation grid c(xmin, ymin, xmax, ymax). #' @param model Character vector indicating the geometry model to use for st_intersection boundaries. Default = "semi-open"; for details see ?s2::s2_options +#' @param include_tile_center Should the center point of each tile be included in the output? This helps with plotting data because the center point for grid locations after intesecting with the stratum shapefile is not always the center of a grid cell since intersected cells are not always squares. #' @export -make_2d_grid <- function(obj, resolution = c(3704, 3704), output_type = "point", bbox = NULL, model = "semi-open") { +make_2d_grid <- function(obj, resolution = c(3704, 3704), output_type = "point", bbox = NULL, model = "semi-open", include_tile_center = FALSE) { # Extract CRS obj_srid <- sf::st_crs(obj, parameters = TRUE)$srid @@ -62,7 +63,22 @@ make_2d_grid <- function(obj, resolution = c(3704, 3704), output_type = "point", # Convert to polygon and find intersection with obj interp_polygons <- terra::as.polygons(interp_grid) |> - sf::st_as_sf(crs = obj_srid) |> + sf::st_as_sf(crs = obj_srid) + + # Add tile center coordinates for plotting + if(include_tile_center) { + + coords_for_plots <- sf::st_centroid(interp_polygons) |> + sf::st_coordinates() |> + as.data.frame() |> + dplyr::rename(lon_plot = X, lat_plot = Y) + + interp_polygons <- interp_polygons |> + dplyr::bind_cols(coords_for_plots) + + } + + interp_polygons <- interp_polygons |> sf::st_intersection(obj, model = model) # Convert to square kilometers diff --git a/assets/ex_2d_grid_with_plot_centroids.png b/assets/ex_2d_grid_with_plot_centroids.png new file mode 100644 index 0000000..760d805 Binary files /dev/null and b/assets/ex_2d_grid_with_plot_centroids.png differ diff --git a/assets/make_plot_grid.md b/assets/make_plot_grid.md new file mode 100644 index 0000000..13f8bcf --- /dev/null +++ b/assets/make_plot_grid.md @@ -0,0 +1,34 @@ +# Return coordinates for plotting gridded data using make_2d_grid() + +### _Feature added in 3.3.0 (see [NEWS](/NEWS))_ + +This example show how to use the make_2d_grid() function to produce grid cell centroid +coordinates that facilitate plotting data from 2d grids that are used for spatial +analyses (e.g. using VAST or sdmTMB). Otherwise this would require post-processing +manipulation of gridded data. + +The plotting coordinate columns (lon_plot and lat_plot) should only be used for +visualization. + +``` r +library(akgfmaps) + +ebs_layers <- akgfmaps::get_base_layers(select.region = "sebs", + set.crs = "EPSG:3338") + +ebs_grid <- akgfmaps:::make_2d_grid(obj = ebs_layers$survey.strata, + resolution = res, + bbox = vast_bbox, + output_type = "point", + include_tile_center = TRUE) + +p1 <- ggplot() + + geom_tile(data = ebs_grid, + mapping = aes(x = lon_plot, + y = lat_plot, + fill = Stratum)) + +p1 +``` + +![](/assets/ex_2d_grid_with_plot_centroids.png) diff --git a/man/make_2d_grid.Rd b/man/make_2d_grid.Rd index f996a49..2036cef 100644 --- a/man/make_2d_grid.Rd +++ b/man/make_2d_grid.Rd @@ -9,7 +9,8 @@ make_2d_grid( resolution = c(3704, 3704), output_type = "point", bbox = NULL, - model = "semi-open" + model = "semi-open", + include_tile_center = FALSE ) } \arguments{ @@ -22,6 +23,8 @@ make_2d_grid( \item{bbox}{Optional. A 4L numeric vector to define the edges of the interpolation grid c(xmin, ymin, xmax, ymax).} \item{model}{Character vector indicating the geometry model to use for st_intersection boundaries. Default = "semi-open"; for details see ?s2::s2_options} + +\item{include_tile_center}{Should the center point of each tile be included in the output? This helps with plotting data because the center point for grid locations after intesecting with the stratum shapefile is not always the center of a grid cell since intersected cells are not always squares.} } \description{ Creates a raster (terra 'SpatRaster'), point grid (sf 'POINT'), or polygon mesh/'fishnet' (sf 'POLYGON') from an sf polygon. diff --git a/test_3_0_0.R b/test_3_0_0.R deleted file mode 100644 index d7988ba..0000000 --- a/test_3_0_0.R +++ /dev/null @@ -1,221 +0,0 @@ -# Testing updated functions for akgfmaps version 3.0 -# Updates for 3.0: -# - Removed raster and sp as dependencies and eplaced raster functions with terra functions ahead of the 2023 retirement of rgdal and sp. Functions that previously generated RasterLayer and RasterBrick objects (raster package) now produce SpatRaster objects (terra package). -# - Added a make_idw_stack() function that creates multilayer SpatRaster objects where interpolations are conducted by grouping variables (such as YEAR). Intended to reduce unnecessary run-times and post-hoc manipulation when generating maps for multiple years. -# - Updated vignettes and demos to reflect changes. - - -library(akgfmaps) - -test1 <- make_idw_map( - x = akgfmaps::YFS2017, - COMMON_NAME = "yellowfin sole", - LATITUDE = NA, - LONGITUDE = NA, - CPUE_KGHA = NA, - region = "bs.south", - extrap.box = NULL, - extrapolation.grid.type = "stars", - set.breaks = "jenks", - grid.cell = c(1e4, 1e4), - in.crs = "+proj=longlat", - out.crs = "EPSG:3338", - key.title = "auto", - log.transform = FALSE, - idw.nmax = 4, - use.survey.bathymetry = TRUE, - return.continuous.grid = TRUE -) - -test1$plot - -test2 <- make_idw_map( - x = akgfmaps::YFS2017, - COMMON_NAME = "yellowfin sole", - LATITUDE = NA, - LONGITUDE = NA, - CPUE_KGHA = NA, - region = "bs.south", - extrap.box = NULL, - extrapolation.grid.type = "sf", - set.breaks = "jenks", - grid.cell = c(1e4, 1e4), - in.crs = "+proj=longlat", - out.crs = "EPSG:3338", - key.title = "auto", - log.transform = FALSE, - idw.nmax = 4, - use.survey.bathymetry = TRUE, - return.continuous.grid = TRUE -) - -test2$plot - -test3 <- make_idw_map( - x = akgfmaps::YFS2017, - COMMON_NAME = "yellowfin sole", - LATITUDE = NA, - LONGITUDE = NA, - CPUE_KGHA = NA, - region = "bs.south", - extrap.box = NULL, - extrapolation.grid.type = "sf.simple", - set.breaks = "jenks", - grid.cell = c(1e4, 1e4), - in.crs = "+proj=longlat", - out.crs = "EPSG:3338", - key.title = "auto", - log.transform = FALSE, - idw.nmax = 4, - use.survey.bathymetry = TRUE, - return.continuous.grid = TRUE -) - -test3$plot - -plot(test1$continuous.grid) -plot(test2$continuous.grid) -plot(test3$continuous.grid) - - -# Testing 2D grid generation -test4 <- make_2d_grid(obj = test3$map_layers$survey.area, - resolution = c(1e4, 1e4), - output_type = "raster", - bbox = NULL, - model = "semi-open") - -test5 <- make_2d_grid(obj = test3$map_layers$survey.area, - resolution = c(1e4, 1e4), - output_type = "polygon", - bbox = NULL, - model = "semi-open") - -test6 <- make_2d_grid(obj = test3$map_layers$survey.area, - resolution = c(1e4, 1e4), - output_type = "point", - bbox = NULL, - model = "semi-open") - - -plot(test4) -plot(test5) -plot(test6) - -test_input <- dplyr::bind_rows( - akgfmaps::YFS2017 |> - dplyr::mutate(YEAR = 2017), - akgfmaps::YFS2017 |> - dplyr::mutate(YEAR = 1999, - CPUE_KGHA = CPUE_KGHA * 0.1) -) - -test7 <- make_idw_stack(x = test_input, - 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 = "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(test7$extrapolation.stack['1999']) -plot(test7$extrapolation.stack['2017']) - -test8 <- make_idw_stack(x = test_input, - 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 = "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 = test8$extrapolation.stack, - mapping = aes(fill = var1.pred)) + - facet_grid(~YEAR) - -names(test8$extrapolation.stack) - -test9 <- make_idw_stack(x = test_input, - 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 = "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 = test9$extrapolation.stack, - mapping = aes(fill = var1.pred)) + - facet_grid(~YEAR) - -map_layers <- akgfmaps::get_base_layers(select.region = "nbs", set.crs = "EPSG:3338") - - -test10 <- reproject_gebco(x = system.file("./extdata/test_files/gebco_2023_n66.0_s60.0_w-177.0_e-159.5.nc", package = "akgfmaps"), - z_varname = "elevation", - z_direction = -1, - raster_template = NULL, - extent = map_layers$survey.area, - resolution = 1e4, - set_crs = NULL, - return_slope_aspect = FALSE, - interpolation_method = "bilinear", - slope_nn = 8, - slope_units = "degrees") - -# Note that default raster::mask() and terra::mask() functions have different behaviors -ggplot() + - geom_stars(data = terra::mask(x = test10, - mask = map_layers$survey.area, - touches = FALSE) |> - stars::st_as_stars()) + - geom_sf(data = map_layers$survey.area, - color = "black", - fill = NA) + - ggtitle(label = "terra::mask(touches = FALSE)") + - scale_fill_distiller(palette = "Spectral", na.value = NA) - -ggplot() + - geom_stars(data = terra::mask(x = test10, - mask = map_layers$survey.area, - touches = TRUE) |> - stars::st_as_stars()) + - geom_sf(data = map_layers$survey.area, - color = "black", - fill = NA) + - ggtitle(label = "terra::mask(touches = TRUE)") + - scale_fill_distiller(palette = "Spectral", na.value = NA) - - -# Generate region guide pdfs -generate_layer_guide() -generate_region_guide() -