hypertidy/anglr

re-orientation

Opened this issue · 1 comments

mdsumner commented

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.

image