[gdal-dev] Multidimensional raster support in GDAL

[Even Rouault]

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171025/d87fc2b7/attachment-0001.html>