[gdal-dev] Multidimensional raster support in GDAL

Ari Jolma ari.jolma at gmail.com
Wed Oct 25 04:42:54 PDT 2017


I'd like to first know / decide what would we mean by a 
"multidimensional raster"?

What we now have as a raster is a dataset with one or more bands. The 
bands represent the data dimension and thus there is one of those 
(2D+1). Would we like to have more data dimensions? What Even sketches 
below is ND, but would we really want to have 2D + ND?

Also now there is no unit for the data dimension and the value is 
strictly increasing positive integer, would we like to have an explicit 
unit and arbitrary value (double, date, time, ...)?

Ari


Even Rouault kirjoitti 25.10.2017 klo 14:01:
>
> Hi
>
> (top posting to clearly mark the start of a new thread)
>
> > I too think that multidimensional raster support would be useful.
>
> > Besides by drastically redesigning data structures, could we get there
>
> > incrementally?
>
> One difficulty is that there are 154 raster drivers that use the 
> current data structures.
>
> 2D is assumed in a number of base data structures. Non exhaustive list:
>
> class GDALDataset:
>
> int nRasterXSize
>
> int nRasterYSize
>
> int GetRasterXSize()
>
> int GetRasterYSize()
>
> GetGeoTransform(double adfGT[6]) // read side
>
> SetGeoTransform(double adfGT[6]) // write side
>
> RasterIO(int nXOff, int nYOff, int nXSize, int nYSize, ... int 
> nBufXSize, int nBufYSize, ...) + IRasterIO()
>
> [...]
>
> class GDADriver
>
> Create( int nXSize, int nYSize, ... )
>
> class GDALRasterBand
>
> int nRasterXSize
>
> int nRasterYSize
>
> int nBlockXSize
>
> int nBlockYSize
>
> int nBlocksPerRow
>
> int nBlocksPerColumn
>
> int GetXSize()
>
> int GetYSize()
>
> GetBlockSize(int* pnBlockXSize, int *pnBlockYSize)
>
> IReadBlock(int nBlockXOff, int nBlockYOff, ...)
>
> RasterIO() + IRasterIO()
>
> [...]
>
> class GDALRasterBlock
>
> int nXOff;
>
> int nYOff;
>
> int nXSize;
>
> int nYSize;
>
> [...]
>
> All GDAL algorithms and utilities assume 2D for subsetting, 
> redimensionning, etc...
>
> In a pure ND model, this would probably be something like:
>
> class GDALDataset
>
> int nAxisCount
>
> char* apszAxisName[nAxisCount] // "X", "Y", "Z", "T", ...
>
> char* apszAxisUnit[nAxisCount] // "deg", "m", "secs since 
> 1970-01-01T00:00:00Z", ...
>
> int panAxisSize[nAxisCount] // generalizes nRasterXSize + nRasterYSize
>
> double padfOrigin[nAxisCount] // generalizes adfGT[0] + adfGT[3]
>
> double papadfAxisVectors[nAxisCount][nAxisCount] // generalizes other 
> coefficients of the geotransform
>
> IMultiDimRasterIO( const int* panWinOffsets, const size_t* 
> panWinSizes, ...., const size_t* panBufSizes )
>
> (some members could actually be only present in the driver class 
> extending GDALDataset. GDALDataset would only have getter and setter 
> similarly to the above GetGeoTransform / SetGeoTransform )
>
> And with that model, we only support regular grids and skewed and 
> rotated grids, but with constant sampling space along each axis. Which 
> might not be enough for the N > 2 'T' or 'Z' dimensions for which I 
> can imagine sampling to be less frequently regular than on the 2D plane.
>
> In netCDF, you can have support for irregular sampling since a 
> variable is indexed by dimensions, and dimensions are generally 
> associated with a 1D variable that describe the coordinate value for 
> each sample point along the axis/dimension.
>
> So a more general model would be:
>
> class GDALDataset
>
> int nAxisCount
>
> char* apszAxisName[nAxisCount]
>
> char* apszAxisUnit[nAxisCount]
>
> int panAxisSize[nAxisCount]
>
> double papadfAxisCoordinates[nAxisCount][] where the size of the 2nd 
> dimension of this array is the value of panAxisSize[i]
>
> (That would still not support fully irregular grids where basically 
> each intersection of the grid should have a coordinate tuple, but we 
> probably don't need to go to that complexity.)
>
> Or perhaps put all axis related stuff in a dedicated class, and with a 
> flag to indicate regular vs irregular spacing, so as to simplify some 
> processing in the regular spacing case
>
> class GDALGridAxis
>
> char* pszName
>
> char* pszUnit
>
> int nSize
>
> bool bRegularSpacing;
>
> double dfOrigin; // if bRegularSpacing == true
>
> double adfVector[nAxisCount]; // if bRegularSpacing == true
>
> double adfAxisCoordinates[nSize]; // if bRegularSpacing == false
>
> In a transition, we'd want:
>
> * all existing 2D functionnalities and public API to be preserved. And 
> some rather mechanical way of converting existing driver code to the 
> new internal API (helpers for the 2D case) since 2D only drivers are 
> and will remain 95% of existing drivers.
>
> * addition of a restricted set of ND functionnalities. Among the 
> potential restrictions for a first stage: read-only support, nearest 
> neighbour resampling. Minimum functionnality needed:
>
> - ND read support in netCDF driver
>
> - A base GDALRasterBand::IMultiDimRasterIO() implementation, which 
> requires the GDALRasterBlock / cache mechanism to support ND
>
> - We'd want some support for gdal_translate to be able to do coverage 
> subsetting and slicing (you'd need to slice down to 2D if no drivers 
> support ND output).
>
> - As the VRT format is the fundamental mechanics for any non trivial 
> gdal_translate operation (any switch besides -of and -co goes through 
> VRT internally), it would need to be updated to support ND.
>
> - C API and Python bindings should be updated.
>
> Even
>
> -- 
>
> Spatialys - Geospatial professional services
>
> http://www.spatialys.com
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171025/b8c7e225/attachment-0001.html>


More information about the gdal-dev mailing list