AustralianAntarcticDivision / angstroms

Tools for model grid data (with raster in R)
http://australianantarcticdivision.github.io/angstroms/
4 stars 1 forks source link

general mesh extraction (WIP) #18

Open mdsumner opened 4 years ago

mdsumner commented 4 years ago
f <- "/rdsi/PUBLIC/raad/data/sose.ucsd.edu/SO6/ITER106/bsose_i106_2008to2012_monthly_Theta.nc"
library(angstroms)
library(quadmesh)

coords <- function(x, idx, ...) {
    stopifnot(is.numeric(idx))
    stopifnot(length(idx) == 2L)
    tnc <- tidync::tidync(x)##ncmeta::nc_dims(x)
    cd <- tidync::hyper_dims(tnc)$name[idx]
    X <- tnc$transforms[[cd[1]]][[cd[1]]]
    Y <- tnc$transforms[[cd[2]]][[cd[2]]]
    list(X, Y)
}

coords_xy <- function(x, ...) {
    coords(x, c(1, 2))
}
coords_xz <- function(x, ...) {
    coords(x, c(1, 3))
}
coords_xt <- function(x, ...) {
    coords(x, c(1, 4))
}
coords_yz <- function(x, ...) {
    coords(x, c(2, 3))
}
coords_yt <- function(x, ...) {
    coords(x, c(2, 4))
}
coords_zt <- function(x, ...) {
    coords(x, c(4, 3))
}

xy <- flip(roms_xy(f, varname = "THETA", slice = c(1, 1)), "y")
xz <- flip(roms_xz(f, varname = "THETA", slice = c(500, 1)), "y")
yz <- flip(roms_yz(f, varname = "THETA", slice = c(200, 1)), "y")
yt <- flip(roms_yt(f, varname = "THETA", slice = c(200, 1)), "y")
zt <- flip(roms_zt(f, varname = "THETA", slice = c(200, 1)), "y")

xx <- coords_xz(f)
cds <- brick(setValues(xz, rep(xx[[1]], nrow(xz))), 
             setValues(xz, rep(rev(xx[[2]]), each = ncol(xz))))

xx <- coords_yz(f)
cds <- brick(setValues(xz, rep(xx[[1]], nrow(xz))), 
             setValues(xz, rep(rev(xx[[2]]), each = ncol(xz))))

xx <- coords_yz(f)
cds <- brick(setValues(yz, rep(xx[[1]], nrow(yz))), 
             setValues(yz, rep(rev(xx[[2]]), each = ncol(yz))))

mesh_plot(flip(yz, "y"), coords = cds)