Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active February 13, 2025 05:01
Show Gist options
  • Save mdsumner/adf6a823019ba1d68aa69b4cf10213c6 to your computer and use it in GitHub Desktop.
Save mdsumner/adf6a823019ba1d68aa69b4cf10213c6 to your computer and use it in GitHub Desktop.

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)

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment