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

rmarkdown::render("oisst_to_grid26.Rmd")

knitr::opts_chunk$set(echo = TRUE)

Optimal Interoplation Sea Surface Temperature (OISST)

This is for Philina's request (issue https://github.com/pbs-assess/pacea/issues/60 ) to project the OISST data onto the same grid (grid26) as used for the BCCM model results. Can adapt some of this into functions if needed later. Doing this quickly here for now.

The OISST data are saved on the original 1/4° grid that NOAA provides them as. So we have not done any projection to get them into pacea; so try that here.

Basing this on approach used in data-raw/roms/roms-data-interpolation.R for the BCCM model results, and recently adopted in data-raw/hotssea/hotssea-data-interpolation.R for the HOTSSea model results.

Specifically, want to transform the oisst_month data, which is an sf object of monthly means in the original spatial 1/4° longitude x 1/4° latitude grid and WGS 84 projection and masked to show only the BC Exclusive Economic Zone area, onto the grid26 grid. The coordinates in the geometry column of the oisst_month object are the centroid of each grid cell.

This means that exactly the same spatiotemporal analyses can be done on both datasets.

load_all() # library(pacea)  # or load_all() when developing
library(dplyr)
library(sf)
library(ggplot2)

In addition to the mean temperature value (weekly or monthly mean) for each grid cell, there are also other statistical measures, such as the standard deviation and number of observations. There is also a start- and end-date column that indicates the temporal period for which the mean is calculated from. Focussing here just on the monthly values.

oisst_month
# plot(oisst_month)
bccm_sst <- bccm_surface_temperature()
bccm_sst

Need to convert the former (means only) into the format of the latter, so with:

1. year-month headings

2. different co-ordinate system

3. and POLYGON instead of POINT geometry:

head(oisst_month$geometry)
head(bccm_sst$geometry)

First generate a subset of data to use, filtering is for speed (now ignoring), selecting here is just what we need:

# sub <- dplyr::filter(oisst_month,
#                     year == 1981) %>%
sub <- dplyr::select(oisst_month,
                     year,
                     month,
                     sst,
                     geometry)
# plot.pacea_oi(sub)    this works

Step 1 - year-month headings

sub_wide <- sub %>%
  tidyr::pivot_wider(names_from = c("year", "month"),
                     values_from = "sst")
sub_wide

Step 2 - different co-ordinate system

sub_trans <- sf::st_transform(sub_wide,
                              crs = "EPSG:3005")

Step 3 - map onto grid26

Attempt 1.

This would be ideal, but requires data to be a terra spatraster. TODO ask Kelsey when she's back. Copying from Kelsey's code in depth.R and editing as necessary.

So this is not run:

#template raster to transform bathy data to match pacea grid
ext2 <- st_bbox(grid26)

temp <- terra::rast(resolution = c(500, 500), #have to resample to lower res than input data
                    xmin = ext2[1] - 500, #match extent + buffer of pacea grid
                    ymin = ext2[2] - 500,
                    xmax = ext2[3] + 500,
                    ymax = ext2[4] + 500,
                    crs = paste0(terra::crs(grid26, describe=TRUE)[,2],
                                 ":",
                                 terra::crs(grid26, describe=TRUE)[,3]))
                                 #crs of pacea grid for zonal statistics

#reproject to new crs, using bilinear for continuous data
sub_trans_proj <- terra::project(sub_trans,
                                 temp,
                                 method="bilinear")
# Error in (function (classes, fdef, mtable)  :
#  unable to find an inherited method for function 'project' for signature '"sf"'

sub_trans_proj

Attempt 2:

Copying from roms-data-interpolation.R and editing:

# interpolate data
# 2 km res
sub_2 <- point2rast(data = sub_trans,
                    spatobj = inshore_poly,
                    loc = c("x", "y"),
                    cellsize = 2000,
                    as = "SpatRast")
# 6 km res
sub_6 <- point2rast(data = sub_trans,
                    spatobj = offshore_poly,
                    loc = c("x", "y"),
                    cellsize = 6000,
                    as = "SpatRast")

# crop out grid cells with polygon masks
sub_sf2 <- sub_2 %>%
  terra::mask(bccm_eez_poly) %>%
  terra::mask(inshore_poly) %>%
  stars::st_as_stars() %>%  ## check here for converting to points (not raster)
  st_as_sf()
sub_sf6 <- sub_6 %>%
  terra::mask(bccm_eez_poly) %>%
  terra::mask(offshore_poly) %>%
  stars::st_as_stars() %>%
  st_as_sf()

# mask 2k grid with 6k grid, then combine grids
sub_sf26a <- sub_sf2[!st_intersects(st_union(sub_sf6),
                                    sub_sf2,
                                    sparse=FALSE,
                                    prepared=TRUE),] %>%
  rbind(sub_sf2[st_intersects(st_union(sub_sf6),
                              sub_sf2,
                              sparse=FALSE,
                              prepared=TRUE),]) %>%
  rbind(sub_sf6)


# Ideally want to make it bigger? For now Philina just wants same as BCCM.
# Need to use same outline as for bccm values, i.e. roms_cave:

## snc_lat <- as.vector(ncvar_get(snc_dat, "lat_rho"))
## svar <- as.vector(ncvar_get(snc_dat, "temp", count = c(-1, -1, 1)))

## sdat <- data.frame(x = snc_lon, y = snc_lat, value = svar) %>%
##   st_as_sf(coords = c("x", "y"), crs = "EPSG:4326") %>%
##   st_transform(crs = "EPSG:3005")

## sroms_cave <- sdat %>%
##   na.omit() %>%
##   concaveman::concaveman()
## sroms_buff <- sdat %>%
##   na.omit() %>%
##   st_geometry() %>%
##   st_buffer(dist = 2000) %>%
##   st_union() %>%
##   st_as_sf()



# Think can ignore roms_buff and roms_cave, which were data-specific
# (e.g. bottom temperature), and we're already using the surface values
sub_sf26b <- sub_sf26a[roms_cave,]

    # 2. use roms_buff to get haida gwaii outline and shore
    sub_sf26b <- sub_sf26b[roms_buff,]

# TODO may need these if domains don't match
    # 3. use default surface roms_cave
#    sub_sf26b <- sub_sf26b[sroms_cave,]

    # 4. use default surface roms_buff
#    sub_sf26 <- sub_sf26b[sroms_buff,]

expect_equal(nrow(sub_sf26b), 41288)   # Passes!
    # data should have 41,288 grid cells

    # assign column names as year_month. TODO seems to already be there, might
    # be formatted slightly differently?
#    names(sub_sf26)[1:(ncol(sub_sf26) - 1)] <- cnames

    # covert to long format data - Don't do long format as it is too big
    # t3_sf26 <- t2_sf26 %>%
    #   tidyr::pivot_longer(cols = !last_col(), cols_vary = "slowest", names_to = "date", values_to = "value")  %>%
    #   mutate(year = substr(date, 1, 4),
    #          month = substr(date, 6, 7)) %>%
    #   dplyr::select(-date) %>%
    #   relocate(value, .after = last_col()) %>%
    #   relocate(geometry, .after = last_col())

    # round to 6 decimal places to reduce file size
    # t3_sf26[, "value"] <- t3_sf26$value %>% # for long format
    #   round(digits = 6)
    oisst_month_grid26 <- sub_sf26b %>%
      st_drop_geometry() %>%
      round(digits = 6) %>%
      st_as_sf(geometry = st_geometry(sub_sf26b))

    # assign pacea class
    class(oisst_month_grid26) <- c("pacea_st", "sf", "tbl_df", "tbl", "data.frame")

    # assign units attribute
    attr(oisst_month_grid26, "units") <- attr(oisst_month, "units")

   save(oisst_month_grid26, file = "oisst_month_grid26.rds", compress = "xz")
# To load (anywhere) before I've saved into package:
# load(paste0(here::here(), "/vignettes/oisst_month_grid26.rds"))
# TODO use_data

  # filename <- paste0("../pacea-data/data/",objname, "_", version, ".rds")
    #filename <- paste0("../pacea-data/data/",objname, ".rds")
    #assign(objname, t3_sf26)

    #do.call("save", list(as.name(objname), file = filename, compress = "xz"))

Test the values. Compare with oisst_month. Plots look different but I don't know how to fix the colour scales:

plot.pacea_st(oisst_month_grid26,     # being explicit about plot here, to help understanding
     year = 1981,
     month = 10)

Original:

plot.pacea_oi(oisst_month,
              year = 1981,
              month = 10)

Figures come out differently.

Figuring out ranges for colour bars, as done slightly differently in each function, but actually the global ranges are a bit different.

c(floor(min(st_drop_geometry(oisst_month_grid26))), ceiling(max(st_drop_geometry(oisst_month_grid26))))
# 3 19
c(floor(min(st_drop_geometry(oisst_month$sst))), ceiling(max(st_drop_geometry(oisst_month$sst))))
# 0 21

Makes sense that interpolated one covers a narrow range, but it also covers a smaller area.

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