gdal_mdim_getView.md

cd gdal/autotest

python3

from osgeo import gdal
ds = gdal.OpenEx('gdrivers/data/netcdf/byte_no_cf.nc', gdal.OF_MULTIDIM_RASTER)
ds.GetRootGroup().OpenMDArray("Band1").GetShape()
# (20, 20)


## create a new dataset with view spec
nds = gdal.MultiDimTranslate("/vsimem/array_view.zarr", ds, format="Zarr", arraySpecs=['name=Band1,view=[::2,::4]'])
## no change from source
nds.GetRootGroup().OpenMDArray("Band1").GetShape()
# (20, 20)

## I expect that to work at command line with 
## gdalmdimtranslate gdrivers/data/netcdf/byte_no_cf.nc out.zarr -of Zarr -array 'name=Band1,view=[::2,::4]'

## calling GetView directly works
ds.GetRootGroup().OpenMDArray("Band1").GetView("[::2,::4]").GetShape()
# (10, 5)

demonstrate successful use of sf utils vs gdal_read_mdim to get a subset with the GDAL mdim arraySpec (this only works the same in both in GDAL >= 3.8.0 and some backport patch of 3.7)

OSGeo/gdal#8297

reported on gdal-dev: https://lists.osgeo.org/pipermail/gdal-dev/2023-August/057591.html

  library(sf)
#> Linking to GEOS 3.12.0dev, GDAL 3.8.0dev-a4c7265179, PROJ 9.1.1; sf_use_s2() is
#> TRUE
dsn <- "/vsicurl/https://github.com/OSGeo/gdal/raw/master/autotest/gdrivers/data/netcdf/byte_no_cf.nc"

## this does not work in GDAL currently, it's fixed in pull/8297 (we get only the 20x20 original array)
gdal_utils("mdimtranslate", dsn, 
           tf <- tempfile(fileext = ".nc"), options = c("-array",  "name=Band1,view=[0:5:1, 0:7:2]"))


sf::gdal_read_mdim(tf, proxy = FALSE)$array_list$Band1
#>      [,1] [,2] [,3] [,4]
#> [1,]  181  255  123  132
#> [2,]  156  132   99  165
#> [3,]  156  140  189  156
#> [4,]  156  156  140  107
#> [5,]  173  140  165  107
#> attr(,"units")
#> [1] ""

## same as

cnt <- c(5, 4)
stp <- c(1, 2)
sf::gdal_read_mdim(dsn, 
                  offset = c(0, 0), count = cnt, step = stp)$array_list$Band1
#>      [,1] [,2] [,3] [,4]
#> [1,]  181  255  123  132
#> [2,]  156  132   99  165
#> [3,]  156  140  189  156
#> [4,]  156  156  140  107
#> [5,]  173  140  165  107
#> attr(,"units")
#> [1] ""

^{Created on 2023-08-31 with reprex v2.0.2}

prior to this fix we see

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.8.0dev-57740a5d1e, PROJ 8.2.0; sf_use_s2() is
#> TRUE
dsn <- "/vsicurl/https://github.com/OSGeo/gdal/raw/master/autotest/gdrivers/data/netcdf/byte_no_cf.nc"

## this does not work in GDAL currently, it's fixed in pull/8297 (we get only the 20x20 original array)
gdal_utils("mdimtranslate", dsn, 
           tf <- tempfile(fileext = ".nc"), options = c("-array",  "name=Band1,view=[0:5:1, 0:7:2]"))


sf::gdal_read_mdim(tf, proxy = FALSE)$array_list$Band1
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
#>  [1,]  181  173  156  189  165   NA   NA   NA   NA    NA    NA    NA    NA
#>  [2,]  156  255  140  140  156   NA   NA   NA   NA    NA    NA    NA    NA
#>  [3,]  156  132  123  165  107   NA   NA   NA   NA    NA    NA    NA    NA
#>  [4,]  156  140   99  132  107   NA   NA   NA   NA    NA    NA    NA    NA
#>  [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#>  [6,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#>  [7,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#>  [8,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#>  [9,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [10,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [11,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [12,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [13,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [14,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [15,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [16,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [17,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [18,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [19,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#> [20,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA    NA
#>       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
#>  [1,]    NA    NA    NA    NA    NA    NA    NA
#>  [2,]    NA    NA    NA    NA    NA    NA    NA
#>  [3,]    NA    NA    NA    NA    NA    NA    NA
#>  [4,]    NA    NA    NA    NA    NA    NA    NA
#>  [5,]    NA    NA    NA    NA    NA    NA    NA
#>  [6,]    NA    NA    NA    NA    NA    NA    NA
#>  [7,]    NA    NA    NA    NA    NA    NA    NA
#>  [8,]    NA    NA    NA    NA    NA    NA    NA
#>  [9,]    NA    NA    NA    NA    NA    NA    NA
#> [10,]    NA    NA    NA    NA    NA    NA    NA
#> [11,]    NA    NA    NA    NA    NA    NA    NA
#> [12,]    NA    NA    NA    NA    NA    NA    NA
#> [13,]    NA    NA    NA    NA    NA    NA    NA
#> [14,]    NA    NA    NA    NA    NA    NA    NA
#> [15,]    NA    NA    NA    NA    NA    NA    NA
#> [16,]    NA    NA    NA    NA    NA    NA    NA
#> [17,]    NA    NA    NA    NA    NA    NA    NA
#> [18,]    NA    NA    NA    NA    NA    NA    NA
#> [19,]    NA    NA    NA    NA    NA    NA    NA
#> [20,]    NA    NA    NA    NA    NA    NA    NA
#> attr(,"units")
#> [1] ""

## same as

cnt <- c(5, 4)
stp <- c(1, 2)
sf::gdal_read_mdim(dsn, 
                   offset = c(0, 0), count = cnt, step = stp)$array_list$Band1
#>      [,1] [,2] [,3] [,4]
#> [1,]  181  255  123  132
#> [2,]  156  132   99  165
#> [3,]  156  140  189  156
#> [4,]  156  156  140  107
#> [5,]  173  140  165  107
#> attr(,"units")
#> [1] ""

^{Created on 2023-09-01 with reprex v2.0.2}