Spatial Areas

Spatial areas available in pacea, and how to create your own

We have some spatial areas saved in pacea, technically as sf (spatial features) polygon or multipolygon data objects. We hope to include many fisheries management areas soon, as being discussed in our GitHub Issue.

We first show how these can be used, and then show how to create your own custom areas using our helper functions.

Current areas

The current areas are:

See the respective help files for details. We plot each of these areas later, but use one of them to explain plotting options. (Developers - see the hidden code here to check for new areas).

Using areas

We can easily restrict spatial data within our plotting functions, or using restrict_to_area() for doing analyses.

We will first use the Optimum Interpolation Sea Surface Temperature (OISST) data described in the OISST vignette,

To just restrict a plot to an area, as an example we use the elusive rarely-used fisheries management area called sasquatch and then plot each March to September from 2014 to 2016:

plot(oisst_month,
     months = 3:9,
     years = 2014:2016,
     area = sasquatch)

We can add options to restrict the plot to the reduced area, and to add an outline (the default is the last month of data):

plot(oisst_month,
     area = sasquatch,
     restrict_plot = TRUE,
     area_outline = sasquatch)

This apparently whimsical plot is instructive in showing that the restriction to the area for oisst_month is based on whether the centre of each cell is within the sasquatch or not. The geometry of the oisst_month data is POINT (in sf speak), not polygons, to save filesize:

sf::st_geometry(oisst_month)
#> Geometry set for 500425 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -138.625 ymin: 46.625 xmax: -123.125 ymax: 54.625
#> Geodetic CRS:  WGS 84
#> First 5 geometries:
#> POINT (-129.125 46.625)
#> POINT (-129.125 46.625)
#> POINT (-129.125 46.625)
#> POINT (-129.125 46.625)
#> POINT (-129.625 46.875)

We can similarly restrict the area of plotting of the BCCM output (the first line either downloads the output or loads it into R; see the BCCM vignette):

primary_production <- bccm_primaryproduction_full()   

plot(primary_production, 
     area = sasquatch,
     restrict_plot = FALSE,
     area_outline = sasquatch)

Looking more closely, we can see that, unlike the oisst_month example above, any cells that are at all inside the area are included:

plot(primary_production, 
     area = sasquatch,
     restrict_plot = TRUE,
     area_outline = sasquatch)

This is because the geometry of the BCCM results are POLYGONS:

sf::st_geometry(primary_production)
#> Geometry set for 177247 features 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 108882.9 ymin: -165260 xmax: 1282883 ymax: 1236740
#> Projected CRS: NAD83 / BC Albers
#> First 5 geometries:
#> POLYGON ((670882.9 1236740, 672882.9 1236740, 6...
#> POLYGON ((672882.9 1236740, 674882.9 1236740, 6...
#> POLYGON ((674882.9 1236740, 676882.9 1236740, 6...
#> POLYGON ((668882.9 1234740, 670882.9 1234740, 6...
#> POLYGON ((670882.9 1234740, 672882.9 1234740, 6...

Internally, plot.pacea_st(), which the plot() calls above are using, is using sf::st_intersects() to calculate intersection locations between the main data set (e.g. primary_production) and the area. The sf::st_intersects() function returns TRUE when features intersect, or even just touch.

For the oisst_month example that uses points as the centre of grid cells, the cell is included in the plot (i.e. is within the area) if the point representing it is within the area. Hence, if the point is outside the area (even if some of the grid cell is within it), then the cell is excluded and we see the white areas inside the sasquatch area.

But for the BCCM primary_production values, the geometry is polygons, and so if any part of a grid cell (represented by a polygon) touches the area then that cell is included. Hence there is no white within the area in this case. Note the differing resolution of the underlying grid cells between the two examples (but this is not relevant here).

You can use restrict_to_areas() to reduce the data to the area of interest, so you can examine the data of interest:

oisst_sasquatch <- restrict_to_area(oisst_month, 
                                    sasquatch)

oisst_sasquatch
#> Simple feature collection with 68707 features and 7 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -131.625 ymin: 47.625 xmax: -128.125 ymax: 52.125
#> Geodetic CRS:  WGS 84
#> # A tibble: 68,707 × 8
#>     year month   sst sst_sd sst_n start_date end_date            geometry
#>    <dbl> <dbl> <dbl>  <dbl> <int> <date>     <date>           <POINT [°]>
#>  1  1981     9 16.1   0.618    30 1981-09-01 1981-09-30 (-130.875 47.625)
#>  2  1981    10 13.9   0.625    31 1981-10-01 1981-10-31 (-130.875 47.625)
#>  3  1981    11 12.0   0.846    30 1981-11-01 1981-11-30 (-130.875 47.625)
#>  4  1981    12  9.15  0.974    31 1981-12-01 1981-12-31 (-130.875 47.625)
#>  5  1981     9 16.1   0.654    30 1981-09-01 1981-09-30 (-130.625 47.625)
#>  6  1981    10 14.0   0.656    31 1981-10-01 1981-10-31 (-130.625 47.625)
#>  7  1981    11 12.0   0.856    30 1981-11-01 1981-11-30 (-130.625 47.625)
#>  8  1981    12  9.23  0.984    31 1981-12-01 1981-12-31 (-130.625 47.625)
#>  9  1981     9 16.1   0.729    30 1981-09-01 1981-09-30 (-130.375 47.625)
#> 10  1981    10 14.1   0.675    31 1981-10-01 1981-10-31 (-130.375 47.625)
#> # ℹ 68,697 more rows

plot(oisst_sasquatch,
     restrict_plot = TRUE)

The plot is the same as that above (for which we set the area at the plotting stage), except the colour scale is slightly different.

For making consequent indices, e.g. the average monthly temperature over the area, see calculate_anomalies(oisst_month, area = sasquatch) in the Ecosystem Summary Heatmaps vignette.

Plots of each area

The easiest way to see an area is to plot its outline over some OISST data (because the scale of the map will be automatically suitable and the coastline is included), such as for Pacific Fisheries Management Area 126:

plot(oisst_month,
     area_outline = pfma_area_126_temp)

All the areas from the table at the top are:

for(i in 1:nrow(areas_table)){
  print(areas_table$`Area name`[i])
  plot(oisst_month,
                    area_outline = get(areas_table$`Area name`[i]))
}
#> [1] "bc_eez"

#> [1] "bccm_eez_poly"

#> [1] "bccm_hotssea_poly"

#> [1] "hotssea_poly"

#> [1] "inshore_poly"

#> [1] "offshore_poly"

#> [1] "pfma_area_126_temp"

#> [1] "sasquatch"

Create your own customised area

You can create your own area (as an sf object) from a matrix or tibble of latitude and longitude values, using our create_area_sf_polygon() function. If the first and last rows are different then the function repeats the first row as the last to make an enclosed polygon. Using a matrix or tibble gives the same result:

my_area_matrix <- matrix(c(-130, -132, -130, -128,
                           48, 50, 52, 50),
                         ncol = 2)
my_area_matrix
#>      [,1] [,2]
#> [1,] -130   48
#> [2,] -132   50
#> [3,] -130   52
#> [4,] -128   50
my_area_1 <- create_area_sf_polygon(my_area_matrix)
my_area_1
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -132 ymin: 48 xmax: -128 ymax: 52
#> Geodetic CRS:  WGS 84
#>                                x
#> 1 POLYGON ((-130 48, -132 50,...

my_area_tibble <- tibble::tibble(
  lon = c(-130, -132, -130, -128),
  lat = c(48, 50, 52, 50))
my_area_tibble
#> # A tibble: 4 × 2
#>     lon   lat
#>   <dbl> <dbl>
#> 1  -130    48
#> 2  -132    50
#> 3  -130    52
#> 4  -128    50
my_area_2 <- create_area_sf_polygon(my_area_tibble)
expect_equal(my_area_1,
             my_area_2)