Changes between Version 88 and Version 89 of WKTRasterTutorial01


Ignore:
Timestamp:
Jul 2, 2010, 6:05:31 AM (11 years ago)
Author:
pracine
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • WKTRasterTutorial01

    v88 v89  
    151151}}}
    152152
    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. Right-clic 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
     155As 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. Right-clic on the added layer and select "Zoom To Layer" to adjust the display to their extent and see them.
    154156
    155157You can also view the complete area covered by the raster coverage as one unique geometry by doing this query:
     
    159161    FROM srtm_tiled;
    160162}}}
    161        
     163
     164ST_Union assemble all the raster envelope into one unique polygon and ST_Buffer make sure it become a simple polygon.
     165
    162166As 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)...
    163167
     
    171175    WHERE rid=3278;
    172176}}}
    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
     178ST_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
     180Vectorizing 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).
    175181
    176182First 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.
     
    319325    CREATE TABLE caribou_srtm_inter AS
    320326    SELECT id,
    321         (st_intersection(rast, the_geom)).geom AS the_geom,
    322         (st_intersection(rast, the_geom)).val
     327        (ST_Intersection(rast, the_geom)).geom AS the_geom,
     328        (ST_Intersection(rast, the_geom)).val
    323329    FROM cariboupoint_buffers_wgs,
    324         srtm_tiled
    325     WHERE st_intersects(rast, the_geom);
     330         srtm_tiled
     331    WHERE ST_Intersects(rast, the_geom);
    326332}}}
    327333       
     
    337343          FROM srtm_tiled,
    338344               cariboupoint_buffers_wgs
    339           WHERE st_intersects(rast, the_geom)
     345          WHERE ST_Intersects(rast, the_geom)
    340346         ) foo;
    341347}}}
    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
     349As 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
     351The 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.
    344352
    345353[[Image(WKTRasterViewOfIntersectionResult.gif, align=right, 50%)]]
     
    351359}}}
    352360
    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.
     361One 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.
    354362
    355363
     
    361369    CREATE TABLE result01 AS
    362370    SELECT id,
    363            sum(st_area(ST_Transform(the_geom, 32198)) * val) /
    364            sum(st_area(ST_Transform(the_geom, 32198))) AS meanelev
     371           sum(ST_Area(ST_Transform(the_geom, 32198)) * val) /
     372           sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev
    365373    FROM caribou_srtm_inter
    366374    GROUP BY id