zipzarr.md

https://eopf-public.s3.sbg.perf.cloud.ovh.net/eoproducts/S02MSIL1C_20230629T063559_0000_A064_T3A5.zarr.zip

In GDAL we have to inner-quote the string like this (a subtle bug that should be fixed IMO):

"ZARR:\"/vsizip//vsicurl/https:/eopf-public.s3.sbg.perf.cloud.ovh.net/eoproducts/S02MSIL1C_20230629T063559_0000_A064_T3A5.zarr.zip\""

Here's a few examples using GDALWarp() api fun from R. Not very fast because there's no overviews in the store.

We're targetting just one dataset in the store '/quality/l1c_quicklook/r10m/tci' and treating that as a classic 2D raster. GDAL's not getting the crs EPSG or wkt2 from the metadata, so I'm not sure if that's lack of compliance by the store or a problem in GDAL.

dsn <- "vrt://ZARR:\"/vsizip//vsicurl/https://eopf-public.s3.sbg.perf.cloud.ovh.net/eoproducts/S02MSIL1C_20230629T063559_0000_A064_T3A5.zarr.zip\":/quality/l1c_quicklook/r10m/tci?a_srs=EPSG:32636"

library(vapour)  ## use -ts arg of GDALWarp() and let it work out the y size
a <- gdal_raster_data(dsn, target_dim = c(1024, 0), bands = 1:3)
ximage::ximage(a, asp = 1)


library(terra)  ## use res*10 
r <- rast(dsn)
rr <- rast(r)
res(rr) <- c(100, 100)
plotRGB(project(r, rr, by_util = TRUE))



library(gdalraster)  ## use -ts arg of GDALwarp()
warp(dsn, "/vsimem/tci.tif", t_srs = "", cl_arg = c("-ts", 1024, 0))

plot_raster(read_ds(new(GDALRaster, "/vsimem/tci.tif")))

In Python I would do this, to read a decimated array the same as above

from osgeo import gdal
gdal.UseExceptions()
dsn = "ZARR:\"/vsizip//vsicurl/https://eopf-public.s3.sbg.perf.cloud.ovh.net/eoproducts/S02MSIL1C_20230629T063559_0000_A064_T3A5.zarr.zip\""


dataset = gdal.OpenEx(dsn, gdal.OF_MULTIDIM_RASTER)
rg = dataset.GetRootGroup()
rg.GetGroupNames()

quality = rg.OpenGroup("quality")
quality.GetGroupNames()
l1c_quicklook = quality.OpenGroup("l1c_quicklook")
l1c_quicklook.GetGroupNames()
r10m = l1c_quicklook.OpenGroup("r10m")
r10m.GetMDArrayNames()

tci = r10m.OpenMDArray("tci")
[dim.GetSize() for dim in tci.GetDimensions()]

a = tci.ReadAsArray(array_start_idx= [0, 0, 0], count=[3, 1098, 1098], array_step = [1, 10, 10])

but I don't have capacity to plot that, so wrapping osgeo.gdal in R I do

dsn = "ZARR:\"/vsizip//vsicurl/https://eopf-public.s3.sbg.perf.cloud.ovh.net/eoproducts/S02MSIL1C_20230629T063559_0000_A064_T3A5.zarr.zip\""
library(pymdim) ## remotes::install_github("mdsumner/pymdim") very very WIP
ds = pymdim:::multidimraster(dsn)
rg = ds$GetRootGroup()
tci <- rg$OpenMDArrayFromFullname("/quality/l1c_quicklook/r10m/tci")
tci
ar <- tci$ReadAsArray(array_start_idx= c(0L, 0L, 0L), count=c(3L, 1098L, 1098L), array_step = c(1L, 10L, 10L))
dim(ar)
ximage::ximage(aperm(ar, c(2, 3, 1)), asp = 1)

(saving the georeferencing from the multidim model for another day, for now)

here with rioxarray, obviously not general for the multidim case but there is no xarray support for multidim (yet)

 import xarray
dsn  = 'ZARR:"/vsizip//vsicurl/https://eopf-public.s3.sbg.perf.cloud.ovh.net/eoproducts/S02MSIL1C_20230629T063559_0000_A064_T3A5.zarr.zip":/quality/l1c_quicklook/r10m/tci'
xarray.open_dataset(dsn, engine = "rasterio")

<xarray.Dataset> Size: 1GB
Dimensions:      (band: 3, x: 10980, y: 10980)
Coordinates:
  * band         (band) int64 24B 1 2 3
  * x            (x) float64 88kB 3e+05 3e+05 3e+05 ... 4.098e+05 4.098e+05
  * y            (y) float64 88kB 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06
    spatial_ref  int64 8B ...
Data variables:
    band_data    (band, y, x) float32 1GB ...