Open mdsumner opened 6 years ago
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)))
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)
Needs testing again in light of new romshcoords/romsdepth implementation, thanks to @raymondben
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)
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
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.
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