[gdal-dev] Multidimensional raster support in GDAL
even.rouault at spatialys.com
Wed Oct 25 04:01:42 PDT 2017
(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
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:
GetGeoTransform(double adfGT) // read side
SetGeoTransform(double adfGT) // write side
RasterIO(int nXOff, int nYOff, int nXSize, int nYSize, ... int nBufXSize, int
nBufYSize, ...) + IRasterIO()
Create( int nXSize, int nYSize, ... )
GetBlockSize(int* pnBlockXSize, int *pnBlockYSize)
IReadBlock(int nBlockXOff, int nBlockYOff, ...)
RasterIO() + IRasterIO()
All GDAL algorithms and utilities assume 2D for subsetting, redimensionning, etc...
In a pure ND model, this would probably be something like:
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 + adfGT
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 /
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
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:
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
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
- 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.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the gdal-dev