42 | | * An untiled image stored in a raster table's row. |
43 | | * 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) |
44 | | * 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) |
| 38 | * ~~ An untiled image stored in a raster table's row ~~ |
| 39 | * ~~ 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) ~~ |
| 40 | * ~~ 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) ~~ |
| 41 | '''UPDATE''': As Pierre suggested, there're only 2 arrangements |
| 42 | * Regulary tiled raster |
| 43 | * Irregulary tiled raster |
71 | | * 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) |
72 | | * Construct the structure for raster bands, with instances of [http://www.gdal.org/classGDALRasterBand.html 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). |
| 74 | * Determine, in a 1st, very fast query to the db, by looking in the raster_overview view, what lower resolution table are available for the requested raster table |
| 75 | * Determine, in a 2nd, fast enough query to the db, the extent and the maximum number of bands of the requested raster be aggregating the extents of all the rasters. This takes about 1 second on 360000 tiles even if there is no index. |
| 76 | * Determine, in a 3rd, very fast query to the db, the pixel size & rotation, the band types and the nodata value for each band of ONLY ONE raster (LIMIT 1). The driver should assume those values will be the same for every other rasters in the table. If when fetching the other tiles, it realizes one does not, we must say that we do not support this arrangement. (I'm still a bit perplex about the nodata value though.) |
74 | | 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: |
75 | | * 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) |
76 | | |
77 | | '''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. |
78 | | |
79 | | '''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? |
80 | | |
81 | | 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. |
82 | | |
83 | | Jorge: Ok. Understood. |
84 | | |
85 | | * Call the functions limiting the output to one result. Fast operation, but may be incorrect |
86 | | |
87 | | '''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. |
88 | | |
89 | | '''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. |
90 | | |
91 | | Currently, the driver takes the first (and slow) option. That caused performance problems (see ticket #497) |
92 | | |
93 | | '''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. |
94 | | |
95 | | 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. |
96 | | |
97 | | Jorge: Yes, it's sufficient. This question is deprecated. Sorry. |
| 78 | '''Open Question''': If in the first query we find a lower resolution table, does the rest of the work must be performed with this lower resolution table? At least these 3 queries, until we want to read the actual raster data to burn it into the buffer. The queries should be faster in an overview table, but the pixel size will not be the same using an overview table instead the normal resolution table. And you don't read from overviews unless you want to implement decimation because your buffer size is different from your raster size. Am I right? |
104 | | * ''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 [http://www.gdal.org/classGDALRasterBand.html#09e1d83971ddff0b43deffd54ef25eef IReadBlock]. |
105 | | |
106 | | '''Pierre question:''' How is the natural block size choosen in our case? |
107 | | |
108 | | '''Jorge''': By ST_Metadata (width and height).In case of non-regular blocked rasters the function raises an error. |
109 | | |
110 | | * 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 [http://www.gdal.org/classGDALRasterBand.html#5497e8d29e743ee9177202cb3f61c3c7 IRasterIO]. |
111 | | |
112 | | We choose the region oriented approach. Please read the Pierre suggested pseudocode and Jorge's comments: |
113 | | |
114 | | '''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. |
| 85 | More specific: |
136 | | Basic GDAL I/O follow this schema: GDALDataSet::RasterIO --> GDALDataset::IRasterIO --> (for each band) 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)[[BR]] |
137 | | |
138 | | 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.[[BR]] |
139 | | |
140 | | 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.[[BR]] |
141 | | |
142 | | 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: |
143 | | * 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) |
144 | | * 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)[[BR]] |
145 | | |
146 | | 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: |
147 | | * 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. |
148 | | * 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?''' '''Response (jorgearevalo): Yes, IWriteBlock persists the blocks to database''' |
149 | | |
150 | | 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). |
151 | | |
152 | | |
153 | | About the IReadBlock method: |
| 107 | About the IReadBlock method (to be implemented in the rasterband class): |