| Version 39 (modified by , 13 years ago) ( diff ) |
|---|
GDAL Driver for PostGIS Raster Working Specifications
Quick Links
Current status of the driver (February 2011)
The driver is:
- Able to read in-db evenly blocked rasters (all blocks with same size)
- Able to read in-db one-row-rasters:
- If the table really has more than one row: using -where clause in connection string
- If the table has more than one row: the table must have been marked as "regularly blocked table", with -k in loader
- Able to generate two kind of raster object based on two modes:
- ONE_RASTER_PER_ROW ('mode = 1' in connection string, or nothing): The default mode. Each table row is considered as an independent raster. If the table required has more than one row, and no -where clause has been specified in connection string, all the table rows will be considered as reported as Subdatasets. Unless you specify the other working mode
- ONE_RASTER_PER_TABLE ('mode = 2' in connection string): Each table is considered as a raster coverage, and each row is a raster tile.
- Too slow (Reads the entire table metadata for constructing the GDALDataset object, and needs one server round per GDALRasterBand::IReadBlock call)
The driver is not:
- Able to read out-db rasters (developed, but not tested, and with known bugs)
- Able to create new rasters
- Able to manage all the PostGIS Raster arrangements
- Able to provide a color interpretation for bands
Design principles
Topic: The basis
The main class of a GDAL driver is GDALDataset: A set of associated raster bands. So, 1 GDALDataset must be able to contain:
- An untiled image stored in a raster table's row.
- A tiled image stored in a raster table (regular or irregular, rectangular or not, with or without missing tiles, with or without overlapping between tiles)
- A raster object coverage from the rasterization of a vector coverage stored in a raster table (regular or irregular, rectangular or not, with or without missing tiles, with or without overlapping between tiles)
In the first case, 1 GDALDataset = 1 PostGIS Raster object. In the other two cases, 1 GDALDataset = Several PostGIS Raster objects. For this reason, the GDAL PostGIS Raster driver has 2 working modes: ONE_RASTER_PER_TABLE, ONE_RASTER_PER_ROW.
However, currently the driver only deals with continuous tiled raster layers, when all the raster tiles are the same size, snap to the same grid and do not overlap (the ideal case).
Open question: Are 2 working modes enough to manage all the raster arrangements?
Pierre: I think yes. We have to distinguish "want we want to produce" from "what we have to deal with". The two modes answer "want we want to produce" and the different table arrangement are "what we have to deal with".
From a GDAL user point of view I know there is a bunch of raster rows in the DB and there is only two things I want to do: extract those rasters rows one by one creating one raster per row or treat them all as a single raster and blend them all together. Furthermore, I want to be able to SELECT those rows using a WHERE statement. If I want a single raster from the db, I have to build my WHERE clause accordingly. There is no need for an extra mode for this. Beside, I don't want to know or have to know what is the raster table arrangement. I expect the driver is able to deal with them all.
Then, the driver has to deal with all the possible arrangement of those selected rows in both mode (overlap, gaps, missing tiles, etc…). You tried to enumerate the posssible arrangement above but I think there is only two cases: the tiles are regularly tiled or they are not, whatever the number of tile there is (1 or more). To me the irregular case is a generalization of the first one.
Jorge: I think we have 3 cases: untiled raster, regularly tiled raster and irregularly tiled raster.
If, and only if, you can optimize the regularly tiled case, then you write is as an exception. The problem is to make sure the table is REALLY regularly tiled without relying on the user knowledge. Just the introduction of the -a option to raster2pgsql.py allowing to append tiles to an existing table make the "regularly blocked" flag untrustable. If really we want to maintain this flag we will have to create something like a ST_ValidateRegularBlocking aggregate function.
Jorge: fully agree. The only way to ensure a raster is regularly tiled is a checking function. To be used carefully.
Pierre: Then if we can not rely on the raster_columns flag and if a ST_ValidateRegularBlocking() would be too slow, we have to treat "regularly tiled" and "irregularly tiled" as a one unique case hoping that the "regular" one will be faster because it involves less processing when merging the tiles together.
Topic: Constructing the GDALDataset object
To construct a GDALDataset object, the driver must:
- Open the dataset (create db connection)
- Read some data about the dataset (metadata): srid, georeference information, projection information, raster data size, band information (number of bands, pixel size, color interpretation, if present), any other driver-specific dataset related information (i.e.: in our case, schema and table name)
- Construct the structure for raster bands, with instances of GDALRasterBand class. You need to provide some basic information: data type (pixel size), block size (GDAL contains a concept of the natural block size of rasters so that applications can organize data access efficiently for some file formats) and color interpretation (if any).
The metadata must be read from the raster table, using SQL functions like ST_Extent (used for raster data extent), ST_Metadata (used for general raster metadata) or functions like ST_SRID, ST_Width, ST_Height, etc. When your GDALDataset matches only one raster row (a raster tile) this is not a problem. But when your GDALDataset matches a whole raster table (ONE_RASTER_PER_TABLE mode), you have 2 options:
- Call the functions over the whole table and filter the result (i.e.: select distinct st_srid(rast) from raster_table, select distinct st_metadata(rast) from raster table). It can be a really slow operation, but you can check if all tiles are like expected (for example: if they are the same size, if they share the same srid, if they overlap or not, etc)
Pierre comment: I think the driver should not try to detect bad raster arrangement with SQL queries. It should just get what it needs from the DB and burn rasters tiles as they come. If it detected that the arrangement is regular then burn them as regular. If they are not, then burn them accordingly.
Jorge comment: What does accordingly mean here? My bet: if the user wants ONE_RASTER_PER_ROW, no problem. Burn one raster file for each tile. If the user wants ONE_RASTER_PER_TABLE and the tiles are not regular, the driver may warn the user and abort or force ONE_RASTER_PER_ROW mode (warning the user first). Any other options?
Pierre's response: Accordingly to their geolocation. 1) create an empty raster buffer the size of the whole query area 2)you query all the raster rows (tiles) you need 3) write pixel values from those tiles to the correct location (deduced from the georeference) in the buffer one after the other. That they are regularly tiled or not does not matter. Just write the last raster last, overwriting existing underlapping values.
Jorge: Ok. Understood.
- Call the functions limiting the output to one result. Fast operation, but may be incorrect
Pierre comment: What might be incorrect? There should not be different srid, pixeltype, or pixelsize in the same table. We have have to warn that we do not support this bad arrangement yet.
Jorge comment: Related with the previous comment, we could simply warn the user about using that band arrangement and maybe force the ONE_RASTER_PER_ROW arrangement instead.
Currently, the driver takes the first (and slow) option. That caused performance problems (see ticket #497)
Open question: How to fetch the information needed to construct the GDALDataset? Pay attention to the fact that you are not asking for raster data yet. You only need metadata, for constructing the basic GDALDataset object.
Pierre: So just knowing the raster extent, the pixelsize and the pixeltype is not sufficient? You could do this with a quick query SELECT ST_Extent(rast), min(ST_BandPixelSize(rast, band)), min(ST_BandPixelType(rast, band)) assuming all tiles have the same pixelsize and pixeltype.
Jorge: Yes, it's sufficient. This question is deprecated. Sorry.
Topic: Reading/Writing raster data
Once constructed the basic structure (GDALDataset object and related GDALRasterBand objects), you need to choose the strategy for raster data reading/writing:
- Natural block oriented r/w: The driver reads/writes data in equal sized blocks. The potentially more efficient way of r/w data. Really, the natural block size for this dataset is chosen during GDALRasterBand creation. So, it's driver's responsibility to provide the desired value for block size. To use this method, your driver must provide an implementation of IReadBlock.
Pierre question: How is the natural block size choosen in our case?
Jorge: By ST_Metadata (width and height).In case of non-regular blocked rasters the function raises an error.
- Region oriented r/w: The driver reads/writes arbitrary regions of data. It's a potentially less efficient method, because you have to take care of data type translation if the data type of the buffer is different than that of the GDALRasterBand. You also must takes care of image decimation / replication if the buffer size (nBufXSize x nBufYSize) is different than the size of the region being accessed (nXSize x nYSize). To use this method, your driver must provide an implementation of IRasterIO.
Pierre question: In which case "the data type of the buffer is different than that of the GDALRasterBand"?
Pierre question: In which case "the buffer size (nBufXSize x nBufYSize) is different than the size of the region being accessed (nXSize x nYSize)"?
Jorge: Ok, these are not use cases of our driver. The point is: if you implement IRasterIO, you must deal with the fact that maybe the caller is using arbitrary data type and buffer size to call your implementation (these are function's arguments). We could simply detect this situation by comparing the data type and buffer size provided as arguments with the known data type and raster size, and raise an error if needed. But I think if the caller uses an integer buffer and we have a raster with floating point values, for example, we still should provide the raster data and warn the user about possible unaccurate values, because of truncating oepration.
Pierre question: Do we agree that "Region oriented" is a general case of "Natural block oriented" and that IReadBlock can be implemented as a wrapper around IRasterIO?
Jorge: Agree. But taking into account the previous comment.
Clearly, there's no best method for reading/writing data in our case. In the ideal case of regulary blocked rasters, with no overlapping and same grid for all tiles, the block oriented r/w is the more appropiate strategy. But in the rest of the cases, a more general r/w method must be provided.
Currently, the natural block oriented r/w method is the one implemented for the driver. This is a limitation for 2 reasons:
- Obviously, it only fits one raster arrangement
- Each ReadBlock call forces a new server round, constructing a Box and getting the raster row that contains it. This can be really slow, in case of huge raster coverages (question raised in ticket #497 too).
Pierre question: Why does each IReadBlock force a new server round and IRasterIO not?
Open question: How to get the needed metadata in case of ONE_RASTER_PER_TABLE arrangement. As argued in ticket #497, executing ST_Extent or ST_Metadata without limits over a big table can be a really heavy process.
Pierre question: How long take a ST_Extent query on a 1 000 000 tiles indexed table?
Jorge: Maybe acceptable time. But 1 000 000 server calls to get the raster data portion surrounded by the box constructed with 4 coords (each IReadBlock round) is too slow, I think. Each IReadBlock call may imply a simple select st_band(rast, nband) from table where… in case of regularly blocked rasters (1 natural GDAL block = 1 raster table row), but the general case (to be implemented in IRasterIO) is different. The caller may want a region covered by 2 raster tiles. Ok, right now we simply support this particular case (regularly blocked raster), but for the more general case, we'd need a st_intersection version returning raster data, IMHO.
Pierre: But how can we avoid the fact that IReadBlock fetch one tile at a time? There seems to be no way around it. It seems to me that in our case IRasterIO is generally more efficient in fetching raster data because it can fetch a larger ectent than a simple tile. Right? How does GDAL decide when to use IReadBlock or IRasterIO? Can't we force it to use IRasterIO all (most) of the time?
You could use a future ST_Union() (not ST_Intersection()) to merge many tiles into one filling the need of IRasterIO but you can also implement your more simple "burn-the-last-fetched-tile-at-the-proper-location" algorythm.
Jorge: The I/O system in GDAL works like this (thanks to Even Rouault, who provided me this explanation some time ago): The GDALRasterBand object has an array of pointers to blocks (papoBlocks). Initially, all the pointers are set to NULL. When a I/O operation over a band is issued, the default implementation in GDALRasterBand::IRasterIO() will divide the requested source window into a list of blocks matching the block size of the band. It will call the GetLockedBlockRef() method of GDALRasterBand to see if there's already a matching block stored in papoBlocks. If yes, the caller will fetch the pixel buffer from the cached block. If not, it will call the IReadBlock() method on the band (a particular driver's implementation. Our driver provides one) and will store the resulting buffer in a new GDALRasterBlock object.
So, GDAL performs its I/O operation in a block-oriented mode. It may be natural to think that a GDAL block = PostGIS Raster tile, and you only need to provide a IReadBlock implementation, but this is a limited point of view. Only works with regular blocking arrangement, like now.
We could get rid of IReadBlock and implement our own IRasterIO version, as you suggested. This IRasterIO implementation would perform the query to get the data. Again, we could call IReadBlock only in case the requested region size matches the tile size, but it would only work if all the tiles have the same size. GDAL block system is not really intended to deal with different size's blocks.
Another important thing is IRasterIO doesn't make any assumption about data buffer size or data type. So, the function must check if the buffer size is enough to store the requested raster data (in case of data read), and if not, use the proper overview to read the data (in case overviews exist), or resample the raster data to fit into the buffer. You must check the data type too. If the raster data type and the buffer data type are different, you must perform the data type translation.
There are some useful functions to deal with these problems, like GDALCopyWords or GDALBandGetBestOverviewLevel. But we're, in some way, breaking the GDAL I/O philosophy.
Open question: What should be the general r/w algorithm? jorgearevalo: I think the strategy read as much data as you can should be the right one, to minimize server rounds. This is: construct a query that, using ST_Intersects, fetches as much rows as possible. This query would be executed in IRasterIO method. But I don't know how to choose the geographic limits for the query (how much data is as much data as you can?)
Pierre question: What do you mean by "I don't know how to choose the geographic limits"? Can't you deduce them from the parameter passed to the function and the metadata of the rasters?
Jorge: Yes, deprecated question. Forget it.
Implementation
Pierre: Can't it not be as simple as the pseudo code below? In the best case the required blocks fits what is in the table and everything is optimized. If not it is slower. We don't have to know in advance whether the table is regularly tiled or not.
GDALRasterBand::IRasterIO(required tile metadata) {
Deduce a world coordinate rectangle which four corners are the center of the upper left, lower left, lower right, upper right pixels of the area requested.
Fetch, in a unique SQL query, every rasters intersecting this area along with the upper left X & Y, width, height, scalex, scaleY, skewX and skewY
If there is only one row fetched and the metadata of this raster fits the required tile metadata (meaning we are querying based on the natural block size)
just copy the values as a block (memcopy) (is this possible? Should we have/isn't there a ST_Bytea() SQL function?)
If not iterate on the required tile
copy pixel values one by one from each raster fetched (that means if there is overlaps, only the pixels of the last raster fetched is burned. We can enhace this later by providing a BLENDING_MODE with the driver)
}
GDALRasterBand::ReadBlock(block x & y){
Deduce IRasterIO required tile metadata
call IRasterIO(required tile metadata)
}
Jorge: I think we should get rid of IReadBlock, because the argument above. About the IRasterIO psc, yes, I think is mostly correct.
