mdsumner
commented
1 month ago next increment, here write write to Zarr and then open that with GDAL multdim, read at /10 and plot
library(reticulate)
xarray <- import("xarray")
pandas <- import("pandas")
ftodate <- function(x) pandas$DatetimeIndex(list(
format(as.Date(stringr::str_extract(basename(x), "[0-9]{8}"), "%Y%m%d"))
))
da <- function(x) {
ds <- xarray$open_dataset(x)
ds$expand_dims(time = ftodate(x))
}
files <- raadtools::read_par(returnfiles = TRUE)
ds <- xarray$concat(lapply(files$fullname[1:10], da), dim = "time")
ds$to_zarr(store = "/tmp/par.zarr")
gdal <- import("osgeo.gdal")
gdal$UseExceptions()
ds <- gdal$OpenEx("/tmp/par.zarr", gdal$OF_MULTIDIM_RASTER)
rg <- ds$GetRootGroup()
rg$GetMDArrayNames()
# [1] "lat" "lon" "time" "palette" "par"
par <- rg$OpenMDArray("par")
par$GetBlockSize()
par$GetDimensionCount()
lapply(par$GetDimensions(), function(.x) .x$GetSize())
b <- par$GetView("[:,:,1]")
b$GetBlockSize()
d <- b$GetDimensions()
lapply(d, function(.x) .x$GetSize())
bytes <- b$Read()
length(bytes)
##10 * 4320 * 2
#m <- matrix(readBin(bytes, "integer", n = 4320 * 10, size = 2), 10)
#ximage::ximage(m)
b <- par$GetView("[0:1,::10,::10]")
b$GetBlockSize()
d <- b$GetDimensions()
shape <- lapply(d, function(.x) .x$GetSize())
sbytes <- b$Read()
length(sbytes)
m <- matrix(readBin(sbytes, "integer", n = length(sbytes)/2, size = 2), shape[[2]] , byrow = TRUE)
ximage::ximage(m, c(-180, 180, -90, 90))
here's the basics of how we could return an xarray object (sometimes you'll want the date from the fileset, sometimes it will be in the object already)
Don't yet know how to specify what the time is (interval or start or middle or whatever)
I think this could avoid the actual calls to open datasets by using kerchunk ...