Add HOTSSea physical oceanographic model for the Salish Sea

[andrew-edwards]

Key bits of discussion with Greig Oldford, @goldford. I think @elisekeppel is one person who would be interested.

Hi Andrew,
We have the HOTSSea model paper almost through peer review and now some DFO salmon researchers want to make use of the outputs. Do you think I could process them like Angelica's BCCM model and get them loaded into PACEA? I could do the processing on my end if there's space there.

HOTSSea v1 is a NEMO-based model for the Salish Sea which has been evaluated for Canadian waters. It fills a niche between the high resolution SalishSeaCast and the lower resoluton BCCM. horizontal resolution is 1.5 km^2 and it's been run from 1980 - 2018. We plan an update from 2019 to present and an extension backwards eventually to 1950.

GMDD - HOTSSea v1: a NEMO-based physical Hindcast of the Salish Sea (1980–2018) supporting ecosystem model development (copernicus.org)

--
Hi Greig, absolutely, this sounds excellent.... We can use the same infrastructure in pacea that we have for Angelica's model, though we set up a grid for that that did not go into Salish Sea. Yes, it would likely take some wrangling on your end, depending on the format of your output. My postdoc Travis did this part of pacea (he doesn't work with me now and is currently away), but we have all the code. I need to do some similar things anyway (people have asked for Angelica's model output extending further south, which we have but just cut if off at the US border).

Be awesome to have the model extended to present and backwards to 1950. Having it in pacea would definitely make it more accessible to the biologists (since they generally use R).

--
Hi Andrew,
Our outputs are probably the same or similar but differ in the grid used. We don't use anything weird like an irregular triangulated grid - just the close-to-square grid cells (~1.5 km horiz).

I likely could interpret and revise the code he used to process the raw model outputs.

I do currently have several fields exported on our model grid, e.g., 3 day means of temps (depth integrated 0 - 4 m - and other depth strata).

--
We could maybe just stick with your ~1.5 km grid, to save any interpolation. Our current grid doesn't actually go into the Salish Sea, and it might not matter if the BCCM grid and yours do not exactly mesh up.

--
Next steps:

• Andy finish reading through manuscript
• Andy work out best approach and what variables (and what temporal resolution) we would like. Think it will be fine to stick with Greig's spatial and depth resolution, and adapt Travis's code for importing it.

[andrew-edwards]

Workflow and next steps (summarising chat with Greig):

We are sure the interpolation procedure (in roms-data-interpolation.R) is doing this:
(1) read NC file data as points
(2) points to raster (Travis's point2rast(), which uses inshore_poly, which we need a hotssea version, see below)
(3) clipping
(4) raster to points (think it might be polygons, check output2)
(5) store as a pacea_st data class - but the only geometry stored is the polygons

Travis makes a grid in make-mask-layer.R and create-grid.R with the key method being st_make_grid() using a polygon created from the model's extent (the buffered hull, or similar). These are based on the ROMS output.

So, need to do something similar for hotssea, and do a 1500 m grid.

• adapt make-mask-layer.R which currently gives inshore_poly, offshore_poly, bccm_eez_poly. That gives a basic polygon of the area (can just plot each one to check). So need a similar poly based on hotssea domain.
• adapt create-grid.R to create a grid also based on hotssea domain, and do a 1500 m resolution.

[andrew-edwards]

Almost got the plotting working, think the coverting into R is good.
See HERE in R/plot-pacea-st-hotssea.R and interpolation one in latest commit.

• fix the axis ranges plus the filled in land
• Greig updating .nc files
• Andy fetch them
• check loading and interpolating works on updated surface temperature, do an .Rmd file with some summaries
• loop round to do all files, figuring out consistent variable names to match bccm
• check plotting works
• upload to Zenodo
• adapt Travis's functions for downloading, from Zenodo instead of from pacea-data
• check plots (maybe need an .Rmd file) and values, maybe reproduce somewhat Figure 9
• add other variables
• write help - check 0 to 4 m detail in next comment
• write vignette, based on bccm one. No need to repeat it all though, just confirm the main parts work.
• add restrict_plotting_range attribute to calc_clim() and calc_anom() based on the input (if it exists), and maybe any other functions. See plot.pacea_st(). Made into new Issue Make calc_clim() and calc_anom() pass on attributes so plotting includes BC coast etc #64.

[andrew-edwards]

More notes from Teams before they helpfully disappear:

Lats and lons in .nc files are centres of the grid.

Figuring out Travis's method, we think:
I was a bit confused but now I'm sure the interpolation procedure is:
(1) read NC file data as points (centre of each grid cell, cells are rotated)
(2) points to raster
(3) clipping
(4) raster to points
(5) store as a special data class - but the only geometry stored is sf as polygons)

Have to trace through code to see where grids get made.

Greig: Until I get the interpolation stuff working more precisely, 0 to 4 m will actually be more like 0 to 4.11 metres and 0 to 30 will be 0 to 31.2 m (for help file, but may get corrected anyway).

[andrew-edwards]

• Add to help/vignette for BCCM: the surface values are at the center of the top grid and varies with the bathymetry, from 0.2 to 1.6 m (from Angelica).

[andrew-edwards]

Regarding the new grid, we had said (top of this issue) that " it might not matter if the BCCM grid and yours do not exactly mesh up.".

Turns out, talking to @ecophilina, that it would really help if the grid does line up with the BCCM grid, for doing spatiotemporal analyses.
Given our 1.5 km grid is somewhat arbitrarily defined, think about whether we can define the new grid to:

• be consistent with the grid_26 grid, by:
• thinking about and defining the start location
• making the cells 1 km square, because then four of them would fit into the 2 x 2 km cells for the BCCM output.

[goldford]

I think it's fine to coarsen the hotssea outputs to 2 km and try to match the exact grid spacing as bccm. If they have to match then I think coarsening is better than down-scaling.

I think there may be a couple ways to achieve making processed hotssea outputs be on the exact same grid_26 grid, and @andrew-edwards might have a better approach than this. But one way might be to create an alt point2rast function (eg point2rast_tbb) that can accept an alternative SpatialObj with which to derive an alternative bounding box (var is 'tbb' in point2rast). My guess is that this SpatialObj would need to be derived from the raw bccm data, ie a 'tdat_sf' object in roms-data-interpolation.R ... but maybe it could be the grid_26 itself.

This would set the bounding box for terra::rast() to be artificially large while passing this function only the data from hotssea. The resulting raster would be on same grid as grid_26. Then, the hotssea data gets cropped (as the code already does) before exporting as a pacea data class.

Again, maybe there is an easier way.

[goldford]

I've revised my thinking here and want to double-check re-interpolating to a different grid resolution is necessary. @ecophilina has two issues: (1) points for oisst vs polys for bccm and (2) different grid resolutions.

I think resolving (1) is a good idea. I think changing the resolution for all datasets to a common one (2) is not a great approach. Instead, a vignette to show users how they may re-interpolate to a grid they choose as the need arises would suffice.

[goldford]

That said, there are good reasons for a common (and nested) grid setup! So, I'm happy to help get that done.