| Version 40 (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.
We choose the region oriented approach. Please read the Pierre suggested pseudocode and Jorge's comments:
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)
}
Jorge: Basically correct, but I need to add some notes here:
Basic GDAL I/O follow this schema: GDALDataSet::RasterIO --> GDALRasterBand::RasterIO --> Divide requested window in blocks, matching GDAL block size, specified in band constructor --> Look for block data in cache --> (if data not in cache) ReadBlock --> IReadBlock (your driver's implementation)
A lot of formats raster formats use tiles, and GDAL is prepared for that, with the block concept. So, if a raster coverage is dividied in tiles, GDAL defines natural block size as the block size that is most efficient for accessing the format, understanding block as raster tile. For this reason, most drivers implement its own version of IReadBlock method. They delegate the I/O system in GDAL core, and only provide blocks of data when required. GDAL, besides, provides a block cache.
The above scheme has a problem: you need a fixed block size. All the blocks of a given coverage have the same width and height. And the regularly blocked scheme is only one of the possible data arrangements PostGIS Raster provides.
To allow all the possible data arrangements, the PostGIS Raster driver must go deeper on GDAL I/O system. It must implement its own version of IRasterIO method, at Dataset level. For two main reasons:
- IRasterIO is a GDALDataset method that allows I/O of an arbitrary size region of the data. This is useful in cases where the blocked I/O isn't the optimal way or even it isn't possible. This is the PostGIS Raster's case (withouht a fixed block size, there's no blocked I/O)
- Getting all the data in a given region (no matter how many raster rows are that) implies only one SQL query. And minimize the server round is interesting to PostGIS Raster (more speed, better performance)
So, basically the psc for IRasterIO is right, but as GDAL is prepared for a blocked I/O with same size blocks and PostGIS Raster can't ensure all the coverage tiles have the same size, we must take a decission:
- Get rid of the GDAL block concept and simply get all the requested data in a query from IRasterIO method. That implies getting rid of GDAL block cache too.
- Adapt our messy arrangement to GDAL block system, by dividing the data fetched from database in equally sized blocks, and manually adding them to the cache. An important thing here is what happens when some raster data, stored in the cache, is modified. Using the GDAL cache system, you can simply replace the data in the cache and mark the block as dirty (using GDALRasterBlock::MarkDirty method), but... Q: Are the dirty blocks written to disk (database, in our case) by IWriteBlock method? If not, how?
After a conversation between Pierre Racine and Jorge Arevalo, we agree on the approach of getting all data in IRasterIO, dividing it into equally sized blocks and storing them in the cache. Jorge Arevalo is working on it (September 2011).
About the IReadBlock method:
GDALRasterBand::!ReadBlock(block x & y){
Deduce IRasterIO required tile metadata
call IRasterIO(required tile metadata)
}
We should get rid of it, for the reasons explained above. We aren't working with equally sized rasters, in general. And even in that case, getting all the requested data in one query with IRasterIO method is a better approach (faster).
