Open mdsumner opened 4 years ago
Here's a different approach, trying to do what angstroms does with NetCDF files with an in-memory brick
time0.tif is derived with
f <- "...data/sose.ucsd.edu/SO6/ITER122/bsose_i122_2013to2017_monthly_Chl.nc"
## lvar says "treat time as level, read all Zs"
br <- raster::brick(f, lvar = 4, stopIfNotEqualSpaced = FALSE)
## use index extent (ignore longlat for the moment)
br <- setExtent(br, extent(0, 2160, 0, 588))
writeRaster(crop(br, extent(800, 1100, 300, 510)), "time0.tif")
## slice a brick
br <- raster::brick("time0.tif")
library(raster)
## slice X =
get_z <- function(x) {
z <- getZ(x)
if (is.null(z)) {
z <- seq(0, -nlayers(x))
}
z
}
slice_x <- function(x, x0 = 1) {
xx <- crop(x, extent(x, 1, nrow(x), x0, x0))
dm <- dim(xx)[dim(xx) > 1]
out <- setValues(raster(nrows = dm[2], ncols = dm[1]), as.vector(values(xx)))
ex <- extent(ymin(xx), ymax(xx), range(get_z(xx)))
setExtent(out, ex)
}
slice_y <- function(x, y0 = 1) {
xx <- crop(x, extent(x, y0, y0, 1, ncol(x)))
dm <- dim(xx)[dim(xx) > 1]
out <- setValues(raster(nrows = dm[2], ncols = dm[1]), as.vector(values(xx)))
ex <- extent(xmin(xx), xmax(xx), range(get_z(xx)))
setExtent(out, ex)
}
mesh_vertical_xy <- function(x, z0 = 0) {
q <- anglr::as.mesh3d(x * 0, triangles = FALSE)
q$vb[3, ] <- z0
q$material$color <- colourvalues::colour_values(values(x))
q$meshColor <- "faces"
q
}
mesh_vertical_xz <- function(x, y0 = 0) {
q <- anglr::as.mesh3d(x * 0, triangles = FALSE)
q$vb[3, ] <- q$vb[2, ]
q$vb[2, ] <- y0
q$material$color <- colourvalues::colour_values(values(x))
q$meshColor <- "faces"
q
}
mesh_vertical_yz <- function(x, x0 = 0) {
q <- anglr::as.mesh3d(x * 0, triangles = FALSE)
q$vb[3, ] <- q$vb[2, ]
q$vb[2, ] <- q$vb[1, ]
q$vb[1, ] <- x0
q$material$color <- colourvalues::colour_values(values(x))
q$meshColor <- "faces"
q
}
library(rgl)
z0 <- 10
y0 <- 100
plot3d(mesh_vertical_xz(slice_y(br, y0), y0 = yFromRow(br, y0)), add = FALSE)
plot3d(mesh_vertical_xy(br[[z0]], get_z(br)[z0]), add = TRUE, alpha = 0.3)
x0 <- 230
plot3d(mesh_vertical_yz(flip(slice_x(br, x0), "x"), x0 = xFromCol(br, x0)), add = TRUE)
The axes here are just x-index, y-index, and z-index (negative), but all works with native coords so it's just a lookup.
We need helpers to drag axes around like this
https://github.com/hypertidy/anglr/wiki/Examples#stacked-walls-plot