following on from https://gist.github.com/mdsumner/4ff89897e95d834073a2f49bf648707f we're working on bindings for the GDAL multidimensional api in R
## this code is very WIP and only in-dev
library(gdalraster) ## a draft branch on mdsumner/gdalraster@multidimnew
dsn <- "/vsigs/gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3"
Sys.setenv("GS_NO_SIGN_REQUEST" = "YES")
ds <- new(GDALMultiDimRaster, dsn)
arr <- ds$getArrayNames()
print(str(arr))
# chr [1:277] "latitude" "level" "longitude" "time" ...
(dims <- ds$getDimensionNames("sea_surface_temperature"))
sizes <- ds$getDimensionSizes("sea_surface_temperature")
## using GDMDArray->GetView
val <- matrix(ds$getViewValues("sea_surface_temperature", "[350616,:,:]"), sizes[2], byrow = TRUE)
dim(val)
## [1] "time" "latitude" "longitude"
range(val, na.rm = T)
# 270.5703 304.8348
lon <- ds$getCoordinateValues("longitude")
lat <- ds$getCoordinateValues("latitude")
ex <- c(range(lon), range(lat)) ## close enough, for now
ds$close()
ximage::ximage(val, ex, asp = 1, col = hcl.colors(128))
maps::map("world2", add = TRUE)