gdal mdim mosaic: accept GeoTIFF as input

[mdsumner]

Feature description

base line support would generate a new dimension and line up the first band from each file in the multidim

next level add support for nominating new dimension name and type of indexing variable, specifying input values for the dimension or or start/step values, a common use-case could provide support for templating date time string from the file/uri name

There are potential complications:

• a huge diversity of creation options and types of dimensions and indexing variables
• meaning of intervals in the new dimension/s
• what to do with multiband?
• maybe this is a new app 'gdal mdim stack ...'

but, this is a very common scenario that implemented by band-proxy in many downstream packages. I'd argue it represents a key missing link from other array schemes in this space, fwiw I'm actively exploring it with workarounds in R and python and what impact a simple or more complex implementation would have on the current GDAL soucrce.

Additional context

GDAL mdim VRT supports input from GeoTIFF, but only (afaics) by manual construction on the VRT XML.

[mdsumner]

I think it makes more sense to add a new utility gdal mdim concat (or stack). Proposed interface:

gdal mdim concat [OPTIONS] <INPUT>... -o <OUTPUT.vrt>

  --input, -i          Input raster files (glob, @filelist, or explicit paths/URIs)
  --output, -o         Output file (.vrt or materialized)
  --output-format, -f  Output format (default: VRT)
  --array-name               Name for the output array (default: "data")
  --band               Source band number (default: 1)

Validation options (matching gdalbuildvrt/gdal raster mosaic):

  --strict              Fail on characteristic mismatch (default: skip with warning)
  --non-strict          Skip mismatched files with warning (default)
  --allow-crs-difference    Allow differing CRS

Dimension options:

  --dim-name           Name for the new dimension (required)
  --dim-values         Explicit numeric values (comma-separated)
  --dim-start          Start value for regular axis
  --dim-increment      Increment for regular axis
  --dim-type           TEMPORAL, VERTICAL, PARAMETRIC, etc.
  --dim-direction      NORTH, SOUTH, UP, DOWN, FUTURE, PAST, etc.
  --dim-unit           Unit string (e.g. "days since 1970-01-01", "meters")

Dimension value extraction:

  --dim-pattern        Regex with one capture group to extract value from filename
  --parse-format       strptime format string for date/time parsing (implies temporal)
  --time-origin        Origin date for temporal offset computation
  --time-unit          seconds | minutes | hours | days (default: days)

Thoughts? I'm working on this.

[rouault]

Thoughts? I'm working on this.

To better understand, it would help if you could give a pseudo gdal info of 2 input datasets , arguments used, and info on the result (don't need to be exact, but to give the idea of the use cases)

[mdsumner]

two bland geotiff inputs

Driver: GTiff/GeoTIFF
Files: /tmp/slice_1.tif
Size is 4, 3
Origin = (0.000000000000000,3.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,   3.0000000) 
Lower Left  (   0.0000000,   0.0000000) 
Upper Right (   4.0000000,   3.0000000) 
Lower Right (   4.0000000,   0.0000000) 
Center      (   2.0000000,   1.5000000) 
Band 1 Block=4x3 Type=Int32, ColorInterp=Gray
  Description = lyr.1
  Min=0.000 Max=11.000 
  Minimum=0.000, Maximum=11.000, Mean=-9999.000, StdDev=-9999.000
  NoData Value=-2147483648
  Metadata:
    STATISTICS_MINIMUM=0
    STATISTICS_MAXIMUM=11
    STATISTICS_MEAN=-9999
    STATISTICS_STDDEV=-9999


Driver: GTiff/GeoTIFF
Files: /tmp/slice_2.tif
Size is 4, 3
Origin = (0.000000000000000,3.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,   3.0000000) 
Lower Left  (   0.0000000,   0.0000000) 
Upper Right (   4.0000000,   3.0000000) 
Lower Right (   4.0000000,   0.0000000) 
Center      (   2.0000000,   1.5000000) 
Band 1 Block=4x3 Type=UInt32, ColorInterp=Gray
  Description = lyr.1
  Min=11.000 Max=22.000 
  Minimum=11.000, Maximum=22.000, Mean=-9999.000, StdDev=-9999.000
  NoData Value=4294967295
  Metadata:
    STATISTICS_MINIMUM=11
    STATISTICS_MAXIMUM=22
    STATISTICS_MEAN=-9999
    STATISTICS_STDDEV=-9999

combine them as depth slices

gdal mdim concat \
 --dim-name "depth" --dim-values 10,20  --dim-type VERTICAL --dim-direction DOWN \
 /tmp/slice_1.tif /tmp/slice_2.tif  depth.vrt

<VRTDataset>
  <Group name="/">
    <Dimension name="depth" size="2" type="VERTICAL" direction="DOWN" indexingVariable="depth"/>
    <Dimension name="y" size="3" type="HORIZONTAL_Y" direction="NORTH" indexingVariable="y"/>
    <Dimension name="x" size="4" type="HORIZONTAL_X" direction="EAST" indexingVariable="x"/>
    <Array name="depth">
      <DataType>Float64</DataType>
      <DimensionRef ref="depth"/>
      <RegularlySpacedValues start="10" increment="10"/>
    </Array>
    <Array name="y">
      <DataType>Float64</DataType>
      <DimensionRef ref="y"/>
      <RegularlySpacedValues start="2.5" increment="-1"/>
    </Array>
    <Array name="x">
      <DataType>Float64</DataType>
      <DimensionRef ref="x"/>
      <RegularlySpacedValues start="0.5" increment="1"/>
    </Array>
    <Array name="slice">
      <DataType>Int32</DataType>
      <DimensionRef ref="depth"/>
      <DimensionRef ref="y"/>
      <DimensionRef ref="x"/>
      <BlockSize>1,3,4</BlockSize>
      <NoDataValue>-2147483648</NoDataValue>
      <Source>
        <SourceFilename>/tmp/slice_1.tif</SourceFilename>
        <SourceBand>1</SourceBand>
        <SourceTranspose>-1,0,1</SourceTranspose>
        <DestSlab offset="0,0,0"/>
      </Source>
      <Source>
        <SourceFilename>/tmp/slice_2.tif</SourceFilename>
        <SourceBand>1</SourceBand>
        <SourceTranspose>-1,0,1</SourceTranspose>
        <DestSlab offset="1,0,0"/>
      </Source>
    </Array>
  </Group>
</VRTDataset>

combine equivalent date/time versions of those as temporal slices (with dates in the filename), with handy templating from file name for the indexing var

(probably we wouldn't want y to be inverted like this just my workaround is very literal for now)

gdal mdim concat \
 --dim-name "time" --dim-pattern "[0-9]{8}" --parse-format "%Y%m%d" --time-unit days \
/tmp/slice_20200101.tif /tmp/slice_20200102.tif time.vrt

<VRTDataset>
  <Group name="/">
    <Dimension name="time" size="2" type="TEMPORAL" indexingVariable="time"/>
    <Dimension name="y" size="3" type="HORIZONTAL_Y" direction="NORTH" indexingVariable="y"/>
    <Dimension name="x" size="4" type="HORIZONTAL_X" direction="EAST" indexingVariable="x"/>
    <Array name="time">
      <DataType>Float64</DataType>
      <DimensionRef ref="time"/>
      <RegularlySpacedValues start="18262" increment="1"/>
      <Unit>days since 1970-01-01</Unit>
    </Array>
    <Array name="y">
      <DataType>Float64</DataType>
      <DimensionRef ref="y"/>
      <RegularlySpacedValues start="2.5" increment="-1"/>
    </Array>
    <Array name="x">
      <DataType>Float64</DataType>
      <DimensionRef ref="x"/>
      <RegularlySpacedValues start="0.5" increment="1"/>
    </Array>
    <Array name="data">
      <DataType>Int32</DataType>
      <DimensionRef ref="time"/>
      <DimensionRef ref="y"/>
      <DimensionRef ref="x"/>
      <BlockSize>1,3,4</BlockSize>
      <NoDataValue>-2147483648</NoDataValue>
      <Source>
        <SourceFilename>/tmp/slice_20200101.tif</SourceFilename>
        <SourceBand>1</SourceBand>
        <SourceTranspose>-1,0,1</SourceTranspose>
        <DestSlab offset="0,0,0"/>
      </Source>
      <Source>
        <SourceFilename>/tmp/slice_20200102.tif</SourceFilename>
        <SourceBand>1</SourceBand>
        <SourceTranspose>-1,0,1</SourceTranspose>
        <DestSlab offset="1,0,0"/>
      </Source>
    </Array>
  </Group>
</VRTDataset>

I had also considered a vrt:// syntax for elevating a raster to an mdim slice (the opposite of ?transpose=arr:1,0, which then would make sense for mdim mosaic I think. There's a deep well of options for covering what anyone might want, but I think we're missing at the moment a formal way to cast a list of tiffs into a temporal multidim (which is realistically the only case I really care about).

[mdsumner]

A motivating set of data is on source.coop where we make COGs of GHRSST:

https://source.coop/ausantarctic/ghrsst-mur-v2

using a dev-R package I published an mdim VRT version of the 'analysed_sst' tifs from that

from osgeo import gdal

gdal.UseExceptions()

ds = gdal.OpenEx(
    "/vsicurl/https://github.com/hypertidy/vrtstack/raw/refs/heads/main/inst/examples/ghrsst-tif-mdim.vrt",
    gdal.OF_MULTIDIM_RASTER
)

rg = ds.GetRootGroup()
arr = rg.OpenMDArray("analysed_sst")

for dim in arr.GetDimensions():
    print(dim.GetFullName())
    print(dim.GetSize())

#/time
#8719
#/y
#17999
#/x
#36000
rg.OpenMDArray("time").ReadAsArray()
#array([11839., 11840., 11841., ..., 20555., 20556., 20557.])