Changes between Version 88 and Version 89 of WKTRasterTutorial01

Timestamp:

Jul 2, 2010, 6:05:31 AM (14 years ago)

Author:

pracine

Comment:

—

WKTRasterTutorial01

@@ -151 +151 @@
 }}}
 
-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.
+"rast::geometry" cast the raster envelope (actually the convex hull) to a geometry thus making OpenJump to display 187200 polygons.
+
+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. Right-clic on the added layer and select "Zoom To Layer" to adjust the display to their extent and see them.
 
 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;
 }}}
-        
+
+ST_Union assemble all the raster envelope into one unique polygon and ST_Buffer make sure it become a simple polygon.
+
 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)...
 
@@ -171 +175 @@
     WHERE rid=3278;
 }}}
-        
-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).
+
+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. 
+
+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).
 
 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
     SELECT id, 
-        (st_intersection(rast, the_geom)).geom AS the_geom, 
-        (st_intersection(rast, the_geom)).val
+        (ST_Intersection(rast, the_geom)).geom AS the_geom, 
+        (ST_Intersection(rast, the_geom)).val
     FROM cariboupoint_buffers_wgs, 
-        srtm_tiled
-    WHERE st_intersects(rast, the_geom);
+         srtm_tiled
+    WHERE ST_Intersects(rast, the_geom);
 }}}
         
@@ -337 +343 @@
           FROM srtm_tiled,
                cariboupoint_buffers_wgs
-          WHERE st_intersects(rast, the_geom)
+          WHERE ST_Intersects(rast, the_geom)
          ) foo;
 }}}
-        
-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.
+
+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.
+
+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.
 
 [[Image(WKTRasterViewOfIntersectionResult.gif, align=right, 50%)]]
@@ -351 +359 @@
 }}}
 
-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.
+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.
 
 
@@ -361 +369 @@
     CREATE TABLE result01 AS
     SELECT id, 
-           sum(st_area(ST_Transform(the_geom, 32198)) * val) / 
-           sum(st_area(ST_Transform(the_geom, 32198))) AS meanelev
+           sum(ST_Area(ST_Transform(the_geom, 32198)) * val) / 
+           sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev
     FROM caribou_srtm_inter
     GROUP BY id