[gdal-dev] Multidimensional raster support in GDAL

[Even Rouault]

On mardi 21 novembre 2017 12:23:31 CET Ari Jolma wrote:
> I'd like to get back to this topic a bit since I learned something from
> the standards that I believe was not discussed in the earlier discussion
> 
> https://lists.osgeo.org/pipermail/gdal-dev/2017-October/047464.html
> 
> In WCS 1 the range of the data - remember the main conceptual division
> to the domain (the X,Y, etc dimensions) and the range (the data if you
> will) - consists of fields, which may consist of bands. That is, a field
> may have an array value (an array of bands). In WCS 2, this separation
> is, firstly, outsourced to SWE since the datatype of range is
> SWE::DataRecord, and, secondly, left as further work in the range
> subsetting extension. The current range subsetting extension makes no
> difference between fields and bands, they are the same and there is no
> structure.

Looking at the WCS 1.1 spec, the RangeSubset KVP syntax shows indeed a lot of complexity. 

RangeSubset = FieldSubset *( “;” FieldSubset )
FieldSubset = FieldId [ “:” Interpolation ] [ “[” AxisSubsets “]” ] AxisSubsets = AxisSubset *( “,” 
AxisSubset )
AxisSubset = AxisId “[” Keys “]”
Keys = Key *( “,” Key )

Not sure if implemetations in practice go beyond having a single field made of a single axis 
(at least this is what MapServer supports).

> 
> The GDAL WCS driver has had some experimental support for time as field

I'd say more than the time support is a kind of ad-hoc dimensional parameter.

> and different time steps as bands.

My reading of the code and doc of the current driver is that time steps are treated as 
subdatasets, and not bands.
See https://trac.osgeo.org/gdal/ticket/3449

> I'm not sure if there are many people
> using it since I believe the functionality is currently broken when used
> together with MapServer. The brokenness seems to be so serious that when
> making a gdal_translate call where the source is a WCS that has a field
> with more than one band (and GDAL recognizes it as such) and the call
> does not involve scaling, the call fails.

Is it really an issue specific of time or something more general about multiband support ?

I've just spent some time looking at mapserver code & doc, and there's some confusion 
regarding multiband support . In particular the WCS 2.0 code has a mode where it uses the 
WCS 1.x mapfile metadata item "wcs_rangeset_axes" in a way not compatible with the WCS 
1.x implementation, see
https://github.com/mapserver/mapserver/pull/5523

However for WCS 1.0, I managed to make band subsetting work with the following conf in 
the mapfile (no time component)

   "wcs_bandcount" "9"
   "wcs_rangeset_axes" "bands"
   "wcs_rangeset_name" "Landsat 5 TM Bands"
   "wcs_rangeset_label" "Bands"
   "wcs_rangeset_description" "Bands for Landsat 5 TM"

DescribeCoverage reports:

    <rangeSet>
      <RangeSet>
        <description>Bands for Landsat 5 TM</description>
        <name>Landsat 5 TM Bands</name>
        <label>Bands</label>
        <axisDescription>
          <AxisDescription>
            <name>bands</name>
            <label>Bands/Channels/Samples</label>
            <values>
              <singleValue>1</singleValue>
              <singleValue>2</singleValue>
              <singleValue>3</singleValue>
              <singleValue>4</singleValue>
              <singleValue>5</singleValue>
              <singleValue>6</singleValue>
              <singleValue>7</singleValue>
              <singleValue>8</singleValue>
              <singleValue>9</singleValue>
            </values>
          </AxisDescription>
        </axisDescription>
        <nullValues>
          <singleValue>0</singleValue>
        </nullValues>
      </RangeSet>
    </rangeSet>

And the following GetCoverage request succeeds and returns 3 bands
SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=multi&FORMAT=GEOT
IFF_8&BBOX=15,48,16,49&bands=9,5,1&CRS=EPSG:4326&WIDTH=5&HEIGHT=5

> The reason seems to be that
> gdal_translate accesses the data band by band.

gdal_translate (actually GDALDatasetCopyWholeRaster()) uses the INTERLEAVE=PIXEL/
BAND metadata item of the IMAGE_STRUCTURE domain as a hint whether to read band by 
band or all bands together

So if you add the following, it should request all bands :
poDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE")

> 
> I guess my point is here is that this is another difference between the
> GDAL data model and the general coverage data model. The simple approach
> to the data model mismatch is to allow GDAL datasets to be made from
> coverage sources with queries. Also one thing I realized I need for the
> new WCS driver is something which distinguishes different GDAL datasets
> made from a single coverage.

I'm not sure what you have in mind,

> 
> Best regards,
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171121/58b3b0ae/attachment-0001.html>