Changes between Version 88 and Version 89 of WKTRasterTutorial01
 Timestamp:
 Jul 2, 2010, 6:05:31 AM (11 years ago)
Legend:
 Unmodified
 Added
 Removed
 Modified

WKTRasterTutorial01
v88 v89 151 151 }}} 152 152 153 As the rasters in the "srtm_tiled" table are in a different spatial reference system than the point you might not see then overlap with the caribou points. Rightclic on the added layer and select "Zoom To Layer" to adjust the display to their extent and see them. 153 "rast::geometry" cast the raster envelope (actually the convex hull) to a geometry thus making OpenJump to display 187200 polygons. 154 155 As the rasters in the "srtm_tiled" table are in a different spatial reference system than the point you might not see them overlap with the caribou points. Rightclic on the added layer and select "Zoom To Layer" to adjust the display to their extent and see them. 154 156 155 157 You can also view the complete area covered by the raster coverage as one unique geometry by doing this query: … … 159 161 FROM srtm_tiled; 160 162 }}} 161 163 164 ST_Union assemble all the raster envelope into one unique polygon and ST_Buffer make sure it become a simple polygon. 165 162 166 As you can see, unlike most raster coverage loaded into a GIS, the area covered by this single table raster coverage is not strictly rectangular and the two big missing square areas are not filled with nodata values. Here you have only one table for the whole coverage which is build from many raster files loaded in a single step. You don't need 13 tables to store 13 rasters. We did this with one band SRTM files in TIF format, but it could have been 3 bands BMP files or JPEG and the number of files and the total size of the coverage does not really matter (PostgreSQL has a limit of 32 terabytes)... 163 167 … … 171 175 WHERE rid=3278; 172 176 }}} 173 174 Vectorizing all the table would have been way too long and the huge amount of geometry would have been impossible to load in OpenJUMP (and most GIS software). 177 178 ST_DumpAsPolygons returns polygons representing raster area sharing a single value and the associated value as a complex "geomval" object. The first item of the SELECT is the geometry and the second the associated value. 179 180 Vectorizing all the table would have been way too long and the huge number of geometry would have been impossible to load in OpenJUMP (and most GIS software). 175 181 176 182 First notice: some pixels are grouped together into complex polygons. ST_DumpAsPolygons() not only blindly dump each pixels, it actually vectorize the raster grouping areas with shared values together. In a raster having a low range of values, this would result in a significantly lower number of polygons than the number of pixels. A higher number of different values (like continuous floating point value rasters) would result in many, many little polygons, most of them square corresponding to only one pixel. A future ST_Reclass() unction should help reducing the number of value in a raster and give more analysis flexibility. For now you can always reclassify your raster BEFORE loading them into PostGIS WKT Raster. … … 319 325 CREATE TABLE caribou_srtm_inter AS 320 326 SELECT id, 321 ( st_intersection(rast, the_geom)).geom AS the_geom,322 ( st_intersection(rast, the_geom)).val327 (ST_Intersection(rast, the_geom)).geom AS the_geom, 328 (ST_Intersection(rast, the_geom)).val 323 329 FROM cariboupoint_buffers_wgs, 324 srtm_tiled325 WHERE st_intersects(rast, the_geom);330 srtm_tiled 331 WHERE ST_Intersects(rast, the_geom); 326 332 }}} 327 333 … … 337 343 FROM srtm_tiled, 338 344 cariboupoint_buffers_wgs 339 WHERE st_intersects(rast, the_geom)345 WHERE ST_Intersects(rast, the_geom) 340 346 ) foo; 341 347 }}} 342 343 This version takes only 7.4 minutes. Almost half the time spent by the preceding one... This little increase in complexity certainly worths the gain in performance. The size of the tiles has also a major effect here. If we restart the whole process with 200 pixels x 200 pixels tiles the intersection operation takes 47 minutes. 100 pixels x 100 pixels: 17 minutes... 30 pixels x 30 pixels: 5 minutes. On the other side, more tiles there is to write in the sql file, longer the loading process. You choose... Smaller tiles is generally a better choice. 348 349 As its geometry/geometry sister, ST_Intersects(geometry, raster, band) returns TRUE if the withvalue area of a raster or a raster tile intersects a geometry and ST_Intersection(geometry, raster, band) returns the geometry/value set of geometries representing the intersection between the geometry and each polygonized group of pixel sharing a same value from the raster and its associated value. It is important to note here that ST_Intersects and ST_Intersection take the nodata value into account. That means if a polygon fall entirely in a nodata value area of a raster, ST_Intersects return FALSE and ST_Intersection returns an EMPTY GEOMETRY. This is a very important raster analysis feature. 350 351 The second version of the query takes only 7.4 minutes to complete. Almost half the time spent by the first one... This little increase in complexity certainly worths the gain in performance. The size of the tiles has also a major effect here. If we restart the whole process with 200 pixels x 200 pixels tiles the intersection operation takes 47 minutes. 100 pixels x 100 pixels: 17 minutes... 30 pixels x 30 pixels: 5 minutes. On the other side, more tiles there is to write in the sql file, longer the loading process. You choose... Smaller tiles is generally a better choice. 344 352 345 353 [[Image(WKTRasterViewOfIntersectionResult.gif, align=right, 50%)]] … … 351 359 }}} 352 360 353 One approach to make this operation without PostGIS WKT Raster would have been to first convert our 13 rasters to 13 shapefiles, import them in PostGIS and do a traditional intersection query between geometries only. The problem is that converting only one of those SRTM file is very difficult... The vector version of only one of those files exceed the limit of 2 GB imposed on shapefiles and the GML equivalent file is more than 10 GB... The WKT Raster approach goes like this: "since it is mostly impossible to convert those raster to vector, load them all into PostGIS as they are and only the smallest part of them will be vectorized to do the requested operation", this using the very efficient GiST index, the versatility of GDAL and the power of the postGIS vector engine.361 One approach to make this operation without PostGIS WKT Raster would have been to first convert our 13 rasters to 13 shapefiles, import them in PostGIS and do a traditional intersection query between geometries only. The problem is that converting only one of those SRTM file is very difficult... The vector version of only one of those files exceed the limit of 2 GB imposed on shapefiles and the GML equivalent file is more than 10 GB... (130 GB for the whole coverage). The WKT Raster approach goes like this: "since it is mostly impossible to convert those rasters to vector, load them all into PostGIS as they are and WKT Raster will vectorized only the what is needed to do the requested operation", this using the very efficient GiST index, the versatility of GDAL and the power of the postGIS vector engine. 354 362 355 363 … … 361 369 CREATE TABLE result01 AS 362 370 SELECT id, 363 sum( st_area(ST_Transform(the_geom, 32198)) * val) /364 sum( st_area(ST_Transform(the_geom, 32198))) AS meanelev371 sum(ST_Area(ST_Transform(the_geom, 32198)) * val) / 372 sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev 365 373 FROM caribou_srtm_inter 366 374 GROUP BY id