generalize slice extraction

[mdsumner]

There is now roms_xz, roms_xy etc. to do the obvious slicing from 4D, but romscoords is still purely XY.

This code is a start

#' convert the depth ramp Cs_r, h (bottom depth), and cell number
#' to a correctly oriented layer of depth values
roms_z <- function(Cs_r, h, cell) {
  out <- flip(raster(matrix(rep(extract(h, cell), each = length(Cs_r)) *  rep(Cs_r, length(cell)), 
         length(Cs_r))), "y")
  setExtent(out, extent(0, ncol(out), 0, nrow(out)))
}
## important to readAll here, else extract is very slow in the loop
h <- readAll(raster(file_db$fullname[1], varname = "h"))
## Cs_r is the S-coord stretching
Cs_r <- rawdata(file_db$fullname[1], "Cs_r")

Can we generalize roms_z to work in any orientation? An example workflow, I think this is right but a bit more checking needed.

It would be good if roms_xz etc had coords_xz arbitrarily, then that matches most needed extractions I think.

## read the 180th (reading up) latitude
## which happens to cut the coast a few times
zz <- roms_xz(f, "temp", slice = c(180, 1))

## read the "height" for this slice (in xz)
## note the top-down conversion for raster cell
zh <- roms_z(Cs_r, h, cellFromRow(romsdata(f, "temp"), 392 - 180 + 1))

plot(romsdata(f, "temp"), asp = "")
plot(zz, asp = "")
 plot(zh, asp = "")
 par(mfrow = c(3, 1))
 plot(romsdata(f, "temp"), asp = "")
 abline(h = 180)
 plot(zz, asp = "")
 contour(zh, add = TRUE)
 plot(zh, asp = "")
 contour(zz, add = TRUE)

All panels in index space, with data or contours providing geographic context.

First panel: surface temperature in (X-Y), horizontal line is Y-slice for next panels at depth
Second panel: temperature at depth (X-Z) - with staggerered depth grid as contours
Third panel: staggered depth (X-Z) with temperature as contours

[mdsumner]

now have functions romsdata2d (replaces romsdata()) and romsdata3d, so then we only need use romshcoords, and raster::stackSelect to get a nearest-layer.

f <- raadtools::cpolarfiles()$fullname[1]
library(angstroms)

## first time slice
avar <- romsdata3d(f, varname = "temp", slice = 1L)

## but each layer is relative to Cs_r/h
## so this provides a 3D brick with the depth at every cell
h3d <- romshcoords(f, S = "Cs_r", depth = "h")

## use calc to find the index of the layer for a depth
depth <- -500  ## note this is absolute, negative matters
layerindex <- calc(h3d, function(x) which.min(abs(x - depth)))

## now we are ready for raster::stackSelect
a <- stackSelect(avar, layerindex)
## don't have the capacity to think deeply enough about this atm
## but we could get the layer below and above ...
#b <- stackSelect(avar, layerindex + 1)
extract_depth_var <- function(x, varname, depth = 0, timeslice = 1L) {
  avar <- romsdata3d(x, varname = "temp", slice = timeslice)
  h3d <- romshcoords(x, S = "Cs_r", depth = "h")
  print(sprintf("calculating layer index at depth %01f", round(depth)))
  ## note how we reverse the pixel sequence, because NetCDF
  ## indexing counts up from negatory (maybe we should flip the h3d ...)
  layerindex <- calc(h3d, function(pixel) which.min(abs(rev(pixel) - depth)))
  print("extracting matching layer pixels ...")
  stackSelect(avar, layerindex)
}

## can vectorize depth in this function to avoid rebuilding h3d
## but need to also do that for time ...
regdepth <- brick(lapply(c(0, -500, -1000, -1500, -2000, -2500), function(a) extract_depth_var(f, "temp", a)))

[mdsumner]

Version 3.

We have to worry about using Cs_r or Cs_w, here Cs_r is appropriate for the mass properties, until we get the interval properly figured out

f <- raadtools::cpolarfiles()$fullname[1]
library(angstroms)

## but each layer is relative to Cs_r/h
## so this provides a 3D brick with the depth at every cell
h3d <- romshcoords(f, S = "Cs_r", depth = "h")

depth_index <- function(depth_array, depth) {
  a <- values(depth_array)
  #ind <- apply(a, 1:2, function(x) which.min(abs(x - depth)))
  ## a group_by will probably be faster here
  ind <- apply(a, 1, function(x) which.min(abs(x - depth)))
  setValues(depth_array[[1]], ind)
}
## this is a faster (memory-bound) stackSelect
layer_select <- function(x, index) {
  longmat <- values(x)
  setValues(index, longmat[cbind(seq_len(ncell(x)), values(index))])
}

## note this depth is absolute, negative matters
index <- depth_index(h3d, -500)

temp <- layer_select(romsdata3d(f, varname = "temp", slice = 1L), 
                     index)
salt <- layer_select(romsdata3d(f, varname = "salt", slice = 1L),
                     index)

[mdsumner]

Needs testing again in light of new romshcoords/romsdepth implementation, thanks to @raymondben

[mdsumner]

These functions work well, now need to combine the 2D coords (quadmesh::mesh_plot style) with the h-coordinates from the coupling work. Too painful to do that from raster, so wrap the romsh with slicing as per xz/yz etc.

library(raadtools)
cf <- cpolarfiles()
library(angstroms)
m <- roms_yz(cf$fullname[18], slice = c(800, 1), varname = "salt")
quadmesh::mesh_plot(m, asp = "")
plot(m, asp = "", col = viridis::viridis(100))
contour(m, add = TRUE)
n <- roms_xz(cf$fullname[18], slice = c(300, 1), varname = "salt")
n <- crop(n, extent(595, 941, 1, 31))
zl <- c(33.5, 35.18)

par(mfrow = c(2, 2))
plot(crop(romsdata2d(cf$fullname[18], "salt", slice = c(31, 1)), ex), zlim = zl, asp = "", col = viridis::viridis(100), main = "surface")
contour(crop(romsdata2d(cf$fullname[18], "salt", slice = c(31, 1)), ex), add = TRUE)

abline(v = 800)
abline(lty = 2, h = 300)
plot(crop(romsdata2d(cf$fullname[18], "salt", slice = c(1, 1)), ex), zlim = zl, asp = "", col = viridis::viridis(100), main = "bottom")
contour(crop(romsdata2d(cf$fullname[18], "salt", slice = c(1, 1)), ex), add = TRUE)
abline(v = 800)
abline(lty = 2, h = 300)
plot(m, asp = "", col = viridis::viridis(100), main = "S-N profile", zlim = zl)
contour(m, add = TRUE)
abline(v = 300, lty = 2)
plot(n, asp = "", col = viridis::viridis(100), main = "E-W profile", zlim = zl)
abline(v = 800, lty = 2)
contour(n, add = TRUE)