In pbs-assess/PACea: An R package of Pacific ecosystem information to help facilitate an ecosystem approach to fisheries management

knitr::opts_chunk$set(echo = TRUE)

Create just a partition??????mapping for test_grid_20

**Okay, this is only for helping with extra partitions now (doing the initial grid separately). Not doing yet, just copied from Create_BC_Partition_New.Rmd, but finally realising we don't need to do this yet (want to test putting data sets onto the grid).

Script for Creating the spatial Partition of BC

This markdown document creates the regular grid of BC's coastline. This regular spatial grid defines how the data are stored internally within PACea. To obtain covariates on other spatial resolutions, the function ffff is used to aggregate/project the values from the grid onto the other regions (typically with polygons of a larger area).

All other regions available for use in PACea are also defined here, with a mapping matrix created for projecting the covariate values from the grid to the regions.

First, we read in all the different regions.

library(dplyr)
load_all()

make_grid_here <- FALSE    # Whether to remake the main BC Grid

# Shapefile of the entire BC Coastline for plotting
Coastline <- sf::st_as_sf(PACea::Coastline)
# This loads in Coastline.rda but converts it, and then changes it further down
# - dangerous! Moving Coastline stuff to coastline.Rmd, and saving a final
# object as coastline.rda probably.

# Very coarse, may not need:
# Coastline_Simp <- sf::st_simplify(sf::st_buffer(Coastline, dist=20),
#                                  dTolerance=25)

# Coarse polygons defining 4 Major regions: QCS, WCVI, HS, WCHG.
BC_Major_Area_Boundaries <- sf::st_as_sf(PACea::BC_Major_Area_Boundaries)
BC_Major_Area_Boundaries <- sf::st_transform(BC_Major_Area_Boundaries,
                                             sf::st_crs(Coastline))

# Name the individual polygons
BC_Major_Area_Boundaries$Poly_Name <- c('HS',
                                        'QCS',
                                        'WCHG',
                                        'WCVI')

# Coarse polygons defining the Pacific Herring Sections (from SpawnIndex package)
Pacific_Herring_Sections <- sf::st_transform(PACea::Pacific_Herring_Sections,
                                             sf::st_crs(Coastline))

# Convert Section to Statistical Area (first two digits)
Pacific_Herring_Sections$Statistical_Area <-
  factor(substr(Pacific_Herring_Sections$Section, 1, 2))

# Merge by Statistical Area
Pacific_Herring_Sections <- Pacific_Herring_Sections %>%
  dplyr::group_by(Statistical_Area) %>%
  dplyr::summarize()

# Name the individual polygons
Pacific_Herring_Sections$Poly_Name <-
  paste0('Herring_SA_',
         Pacific_Herring_Sections$Statistical_Area)

which.max(diff(sf::st_coordinates(Coastline)[,2]))
sf::st_coordinates(Coastline)[1333:1337,]
sf::st_bbox(Coastline)

# expand the polygon by 50km Westwards
poly_expansion <-
  sf::st_polygon(list(matrix(c(390,1036.5509,
                               339,1036.5509,
                               1030,253.6858,
                               1080,253.6858,
                               390,1036.5509),
                             ncol=2,
                             byrow = T)))
poly_expansion <- sf::st_sf(sf::st_sfc(poly_expansion),
                            crs = sf::st_crs(Coastline))

Coastline <- sf::st_union(poly_expansion,
                          Coastline)
plot(Coastline)

if(make_grid_here)
{
 ----
moving to make-grid.R, keeping here for now to test.
  # create large pixel grid over the simplified coastline
  # Rotate the grid to maximize the spatial overlap
  rotang = 318.5

  center <- sf::st_centroid(sf::st_union(Coastline))
  # This is the function call which actually creates the master grid
  # To increase the resolution (almost certainly desired), increase the
  # n argument (e.g. n=c(100,100)).
  Coastline_Grid <- sf::st_sf(
                          sf::st_make_grid(tran(sf::st_geometry(Coastline),
                                                -rotang, center),
                                           n=c(10,10)))

  Coastline_Grid <- tran(sf::st_geometry(Coastline_Grid),
                         rotang,
                         center)

  # Keep only the polygons that fall within or touch the original Coastline to save storage space
  Coastline_Grid <- sf::st_as_sf(sf::as_Spatial(Coastline_Grid))
  sf::st_crs(Coastline_Grid) <- sf::st_crs(Coastline)

  # Name the individual polygons by row-column method
  # Note that sf works from bottom-left to top-right
  # The first grid cell is the bottom-left grid cell (1,1)
  # The final grid cell is the top-right grid cell (n,n)
  # NEED TO CHANGE THE VALUES 10 TO MATCH THE CHOICE OF n=c(,) above!!
  Coastline_Grid$Poly_Name <- paste0('Grid_Index_(',
                                     rep(1:10, each=10), # repeat each value of 1:cols, nrows
                                     ',',
                                     rep(1:10, times=10), # repeat sequence 1:nrows, ncols times
                                     ')')

  Coastline_Grid  <- Coastline_Grid[Coastline, ]

----
  plot(Coastline_Grid)
  plot(Coastline, add=T) # , col='red')

  usethis::use_data(Coastline_Grid, overwrite = T)
}
if(!make_grid_here)
{
  Coastline_Grid <- PACea::Coastline_Grid
}

Next, we bind the competing region definitions into a single spatial object called BC_Partition. We name the BC grid, which will eventually store the data, as 'BC Grid' within the variable Regions_Name. We then define a matrix to perform the aggregations from the BC Grid to the other regions' polygons

The basic idea of the mapping matrix is as follows. The number of rows of the matrix equals the number of cells of the BC grid (in which the data is stored). The number of columns of the matrix equals the number of polygons across all the regions stored in PACea (e.g. Herring Sections). Then, a given column of the matrix has entries equal to the fraction of its total area which is contained within each grid cell. These fractions sum to 1.

AE: Maybe better to do that as region_name_partition for each one, doesn't have to be a giant matrix

Thus, to evaluate the average value of a covariate (e.g. SST) in a specific region, simply sum the average values within each grid cell, multiplied by the fraction of the total area of intersection. Because these sum to 1, the quantity is a spatially-weighted average. Of course, this only works if the data are not missing in any grid cell, else the fractions need to be re-weighted.

Furthermore, if a sum or total of a quantity over a region is desired instead of a mean concentration or density, then we simply scale by the area of the region - stored as Ocean_Intersection_Areas, assuming the quantity can only exist in the water.

# If additional regions are added, simply add them to the bind_rows
# function call. Note that BC Grid MUST GO FIRST
All_regions <- dplyr::bind_rows(
                        Coastline_Grid %>% dplyr::mutate(Regions_Name =
                                                         'BC Grid'),
                        BC_Major_Area_Boundaries %>% dplyr::mutate(Regions_Name =
                                                                   'BC Major Areas'),
                        Pacific_Herring_Sections %>% dplyr::mutate(Regions_Name = 'Pacific Herring Spawn Statistical Areas'))

#BC_Partition <-
#  sf::st_intersection(All_regions)

# Convert the GEOMETRYCOLLECTION to multipolygon
#BC_Partition <-  sf::st_collection_extract(BC_Partition, type = c("POLYGON"))

# Create a Poly_ID variable
BC_Partition <- All_regions %>%
  mutate(Poly_ID = row_number()) %>%
  select(Poly_Name,
         Regions_Name,
         Poly_ID)

# Compute the areas of each region
region_Areas <- sf::st_area(BC_Partition)

# Compute the area of intersection between the (approximate) ocean polygon and each of the regions
Ocean_Intersection_Areas <- rep(0,
                                dim(BC_Partition)[1])

for(i in 1:dim(BC_Partition)[1]){
  Ocean_Intersection_Areas[i] <-
    ifelse(is.null(sf::st_area(sf::st_intersection(BC_Partition[i,],
                                                   Coastline,
                                                   drop=F))),
           0,
           sf::st_area(sf::st_intersection(BC_Partition[i,],
                                           Coastline,
                                           drop=F)))
}

# Check that all the regions are POLYGON or MULTIPOLYGON type
sf::st_geometry_type(BC_Partition)

plot(BC_Partition)

# Define the mapping between each region with the BC Grid
# Note - the matrix will also be a useful concept in future, when
# users wish to define new regions which are union of existing regions
Mapping_Matrix <- #diag(nrow=dim(BC_Partition)[1])
  matrix(0,
         nrow = dim(Coastline_Grid)[1],#dim(BC_Partition)[1],
         ncol = dim(BC_Partition)[1])

for(i in 1:dim(BC_Partition)[1]){
  # loop through the partitioned regions and compute the area of
  # intersection with Coastline Grid.
  Mapping_Matrix[sf::st_intersects(Coastline_Grid,
                                   BC_Partition[i,],
                                   sparse = F),
                 i] <-
    as.numeric(sf::st_area(sf::st_intersection(Coastline_Grid,
                                               BC_Partition[i,])))

  # Divide by rowSums to normalize to 1 (to create spatially-weighted mean)
  if(sum(Mapping_Matrix[,i])>0){
    Mapping_Matrix[,i] <- Mapping_Matrix[,i] / sum(Mapping_Matrix[,i])}

  if(sum(Mapping_Matrix[,i])==0){
    Mapping_Matrix[,i] <- 0}
}

# Check the mapping has been made correctly
Mapping_Matrix
colSums(Mapping_Matrix)

# create a named list with each entry equal to a index vector
# each index vector points the set of regions to the correct rows
# Mapping Matrix  (e.g. WCHG, WCVI, HS, QCS)
index_vectors <-
  list(Coastline = which(BC_Partition$Regions_Name == "BC Grid"),
       BC_Major_Area_Boundaries = which(BC_Partition$Regions_Name == "BC Major Areas"),
       Pacific_Herring_Sections = which(BC_Partition$Regions_Name == "Pacific Herring Spawn Statistical Areas"))


# save the results for future use
BC_Partition_Objects <-
  list(BC_Partition = BC_Partition,
       region_Areas = region_Areas,
       Ocean_Intersection_Areas = Ocean_Intersection_Areas,
       Mapping_Matrix = Mapping_Matrix,
       index_vectors = index_vectors)

usethis::use_data(BC_Partition_Objects,
                  overwrite = T)

pbs-assess/PACea documentation built on April 17, 2025, 11:36 p.m.