In afsc-gap-products/akgfmaps: Alaska Groundfish and Ecosystem Survey Area Mapping

Overview

• Advantages and disadvantages of mapping and spatial analysis in R
• Lots and lots of examples
• Older packages: sp, rgdal
• {sf} package
• {raster} package
• {terra} package
• {stars} package

Advantages of R for mapping and spatial analysis

• Reproducibility
• Transparency
• Transferability
• Efficiency through automation

Disadvantages of R for mapping and spatial analysis

• Less intuitive than GIS software (in my opinion)
• Lots of packages and 'ecosystems'
• Pros and cons of packages are ambiguous to new users
• Shifting landscape
• Takes longer to do some things

Load libraries

library(akgfmaps)
library(magrittr)
library(ggplot2)
library(here)
library(sp)
library(rgdal)

library(akgfmaps)
library(magrittr)
library(ggplot2)
library(here)
library(sp)
library(rgdal)
library(shadowtext)

{rgdal} Reading shapefiles, bindings to spatial libraries

This package will be retired by the end of 2023

• Provides bindings to Geospatial Data Abstraction Library (GDAL)
• Provides access to PROJ library

ebs_strata_rgdal <- system.file("data", "bs_grid_w_corners.shp", package = "akgfmaps") |> 
  rgdal::readOGR()

{sp}: Classes and Methods for Spatial Data

• Relies on rgdal for bindings to GDAL and access to PROJ
• Different classes for different data types (e.g. SpatialPolygonDataFrame, SpatialPointsDataFrame)
• Data structure: nested lists
• Not based on ISO standards
• Doesn't play nicely with ggplot (but it has gotten better!)
• Less intuitive (my opinion)

The {sf} package

• Developed by the same folks as sp
• Provides bindings to GDAL and PROJ
• Data structure: Object classes sf and data.frame
• Data structure: uses ISO standard well-known text strings for geometries and coordinate reference system
• Objects work like a data.frame: filter/subset, merge/join, aggregate/summarize
• Faster than sp

{sf}: Reading shapefiles

• For nonraster files (e.g., points, polygons, lines)

ebs_strata <- system.file("data", "ebs_strata.shp", package = "akgfmaps") |> 
  sf::st_read()

{sf}: Data structure

print(ebs_strata)

{sf}: Plotting using base R

plot(ebs_strata)

{sf}: Plotting using ggplot2

• geom_sf(): ggplot2 geometry for sf
• coord_sf(): coordinate system for sf

ggplot() +
  geom_sf(data = ebs_strata,
          aes(fill = factor(Stratum))) +
    geom_sf_text(data = ebs_strata,
          aes(label = Stratum))

{sf}: Creating sf objects

Example using EBS shelf yellowfin sole CPUE from 2017

print(head(akgfmaps:::YFS2017))

yfs_df <- akgfmaps:::YFS2017 |>
  sf::st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = "EPSG:4326")

print(yfs_df)

Plotted using ggplot2:

ggplot() +
  geom_sf(data = yfs_df,
          aes(color = CPUE_KGHA)) +
  scale_color_viridis()

{sf}: Filtering/subsetting

Same as filtering/subsetting from data.frame objects.

Subset values where CPUE > 0:

sebs_strata <- ebs_strata |>
  dplyr::filter(SURVEY == "EBS_SHELF")

ggplot() +
  geom_sf(data = sebs_strata,
          aes(fill = factor(Stratum)))

{sf}: Coordinate reference systems (CRS)

• Supports PROJ4 and PROJ6 strings
• PROJ6 as well-known text stings

sf::st_crs(yfs_df)
sf::st_crs(sebs_strata)

{sf}: Plotting with mis-matched CRS

ggplot2 does an automatic CRS transformation when there is a mismatch among layers.

ggplot() +
  geom_sf(data = sebs_strata, fill = NA, color = "black") +
  geom_sf(data = yfs_df,
          aes(color = CPUE_KGHA)) +
  scale_color_viridis()

{sf}: Spatial operations with mis-matched CRS

{sf} DOES NOT do an automatic CRS transformation for mismatched geodetic coordinate systems. However, operations will work if geodetic coordinate systems match but projected coordinate systems do not.

Attempt to spatially join points where CPUE >0 to the strata they are in.

yfs_cpuegt0_df <- yfs_df |> dplyr::filter(CPUE_KGHA > 0)

print(try(st_join(sebs_strata , 
                  yfs_cpuegt0_df, 
                  join = st_within), 
          silent = TRUE))

{sf}: Transformation/reprojection

yfs_cpuegt0_df_3338 <- sf::st_transform(yfs_cpuegt0_df, crs = 3338) 

print(head(yfs_cpuegt0_df_3338))

Suggestion: Use EPSG standard coordinate reference systems that don't use WGS84 as a pivot

{sf}: Spatial join

Same as joining/merging from data.frame objects.

yfs_cpuegt0_df_3338 <- sf::st_join(yfs_cpuegt0_df_3338,
                    sebs_strata,
                    join = st_within)

ggplot() +
  geom_sf(data = sebs_strata) +
  geom_sf(data = yfs_cpuegt0_df_3338,
          aes(color = factor(Stratum)))

{sf}: Spatial join

Lots of spatial join functions:

• st_contains_properly
• st_contains
• st_covered_by
• st_covers
• st_crosses
• st_disjoint
• st_equals_exact
• st_equals
• st_is_within_distance
• st_nearest_feature
• st_overlaps
• st_touches
• st_within

{sf}: Non-spatial join

Works just like a data frame. Note: This won't work if both objects are sf.

npac_grid <- system.file("data", "bs_grid_w_corners.shp", package = "akgfmaps") |> 
  sf::st_read()

yfs_grid <- npac_grid |>
  dplyr::inner_join(akgfmaps:::YFS2017, by = "STATIONID")

ggplot() +
  geom_sf(data = yfs_grid,
          aes(fill = CPUE_KGHA),
          color = "black") +
  scale_fill_viridis()

{sf}: Intersection

Clip features to the extent of other features.

Survey grid:

npac_grid <- system.file("data", "bs_grid_w_corners.shp", package = "akgfmaps") |> 
  sf::st_read()

ggplot() +
  geom_sf(data = npac_grid)

{sf}: Intersection (cont.)

Clip the extent of npac_grid to the survey area:

npac_grid <- system.file("data", "bs_grid_w_corners.shp", package = "akgfmaps") |> 
  sf::st_read()

ebs_grid <- sf::st_intersection(npac_grid, ebs_strata)

ggplot() +
  geom_sf(data = sf::st_intersection(npac_grid, ebs_strata))

{sf}: Distance

Calculate great circle or Euclidean distance.

# Great Circle Distance
yfs_4326 <- akgfmaps:::YFS2017 |> 
  sf::st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = "EPSG:4326")
dist_4326 <- sf::st_distance(yfs_4326)

st_is_longlat(akgfmaps:::YFS2017 |> 
  sf::st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = "EPSG:4326") |> 
  sf::st_transform(crs = 3338))

yfs_3338 <- akgfmaps:::YFS2017 |> 
  sf::st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = "EPSG:4326") |> 
  sf::st_transform(crs = 3338)

dist_3338 <- sf::st_distance(yfs_3338)

print(dim(dist_4326))

{sf}: Distance (cont.)

Different distances. Why?

print(dist_4326[1:3,1:3])
print(dist_3338[1:3,1:3])

Because 4326 is a geodetic coordinate system and 3338 is a projected coordinate system. st_distance calculates great circle distance for 4326 and Euclidean distance for 3338.

st_is_longlat(yfs_4326)
st_is_longlat(yfs_3338)

{sf}: Area

Calculate area based on Euclidean or great circle geometry.

sf::st_area(ebs_strata |> sf::st_transform(crs = 3338))[1:3]
sf::st_area(ebs_strata |> sf::st_transform(crs = 4269))[1:3]
sf::st_area(ebs_strata |> sf::st_transform(crs = 4326))[1:3]

{sf}: Length and perimeter

• Use st_length
• Operates on LINESTRING and MULTILINESTRING well-known text.
• Use the st_cast function to convert a POLYGON to a MULTILINESTRING.

st_length(st_cast(ebs_strata, to = "MULTILINESTRING") |> sf::st_transform(crs = 4269))[1:3] 
st_length(st_cast(ebs_strata, to = "MULTILINESTRING") |> sf::st_transform(crs = 3338))[1:3]

{sf}: Centroid

Find the center of geometries.

stratum_center <- sf::st_centroid(ebs_strata)
lab_coords <- sf::st_coordinates(stratum_center) |>
  as.data.frame() |>
  dplyr::mutate(Stratum = ebs_strata$Stratum)


ggplot() +
  geom_sf(data = ebs_strata,
          aes(fill = factor(Stratum))) +
  shadowtext::geom_shadowtext(data = lab_coords,
                              aes(x = X, y = Y, label = Stratum),
                              color = "black", bg.color = "white")

{akgfmaps}: Access to shapefiles

• Provides simple access to many common survey shapefiles
• Applies transformations to layers to simplify transformation
• Shapefiles shared among regions when possible; filters to retrieve subsets (when possible).
• Current support: AI, GOA, EBS shelf, NBS

sebs_layers <- akgfmaps::get_base_layers(select.region = "sebs")
names(sebs_layers)

{akgfmaps}: Features

• akgfmaps::get_base_layers(): Retrieve layers and default plotting options for a region
• akgfmaps::generate_layer_guide(): Create a pdf that shows all of the layers in the package
• akgfmaps::generate_region_guide(): Create a pdf that shows all of the layers for a region

afsc-gap-products/akgfmaps documentation built on April 14, 2025, 7:13 p.m.