[postgis-users] Queries and (helper) functions for viewing and analysing PostGIS Raster data?

[Jorge Arévalo]

Hello Stefan,

On Thu, Nov 25, 2010 at 6:25 PM, Stefan Keller <sfkeller at gmail.com> wrote:
> I'm still trying to understand PostGIS Raster especially regarding
> analysis and viewing PostGIS raster data. Let's begin with the latter.
>
> Viewing:
> => What are the pros & cons to do this (and which is preferred?):
> ST_DumpAsPolygons or ST_PixelAsPolygons?
>
> Here is an example for viewing PostGIS Raster e.g. in OpenJUMP I found
> in the Tutorial
> (http://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01 )
>
> SELECT
>  ST_AsBinary((ST_DumpAsPolygons(rast)).geom),
>  (ST_DumpAsPolygons(rast)).val
> FROM srtm_tiled
> WHERE rid=3278;
>
> This seems to polygonize raster to vector. Unfortunately it's not
> explained in the tutorial what the constraint clause "rid=3278" means:
> it's obviously a single, distinct tile(?). From the book "PostGIS in
> Action": "ST_DumpAsPolygons returns a set of single polygon, pixvalue
> pairs for a given raster band and relies on the GDAL library. In many
> cases this function will be easier and faster to use than going down
> to the level of the pixel using ST_Value.". A solution for doing this
> I found here as a pl/pgsql function called ST_PixelAsPolygons
> (http://trac.osgeo.org/postgis/wiki/WKTRasterUsefulFunctions).
>

On the origins, we simply wanted a function to polygonize a raster, as
a base for a seamless vector-raster intersection. For doing this
purpose, we implemented a C function that uses the GDAL Polygonize
algorithm (http://trac.osgeo.org/gdal/browser/trunk/gdal/alg/polygonize.cpp).
This core function is the one called by ST_DumpAsPolygons.

The ST_PixelAsPolygons function is a function that polygonizes a
raster too, but it has 2 differences with ST_DumpAsPolygons:
- ST_PixelAsPolygons is a pure PL/pgSQL implementation, and
ST_DumpAsPolygons relies on a core C function.
- The algorithm is different in both functions. ST_DumpAsPolygons try
to collect all neighboring pixels with the same value in one polygon.
ST_PixelAsPolygons transforms each pixel in a geometry.

So, ST_PixelAsPolygons was a "temporary" function, until having a
working and more efficient implementation of raster polygonization:
ST_DumpAsPolygons. I'll use ST_DumpAsPolygons (I'd like to add it some
improvements, but I don't know when)

About the clause "rid=3278", it was only for practical reasons: It
takes a long time to polygonize big raster coverages, and can generate
a big amount of polygons. So, Pierre polygonized only one raster tile,
to show the result in the tutorial. One raster tile, in PostGIS Raster
context, means "1 raster table row". And remember there's no real
difference between "raster" and "tile" in PostGIS Raster. Each row of
a raster table can be treated as a single raster object, because it
has its own georeference information, even when this row belongs to a
bigger raster coverage (a raster table).


>
> Analysis: Regarding analysis I'd like to begin with the observation,
> that all examples I've found so far are polygonizing rasters to
> vector. From the tutorial doing overlay:
>
> CREATE TABLE caribou_srtm_inter AS
> SELECT id,
>  (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);
>
> => What is the intended exact result type of this overlay operation
> query (ST_Intersects)? POINT?
> => Which query would generate another raster layer (instead of a point
> vector like in the above example)?
>

The intersection function is based on ST_DumpAsPolygons function. So,
the result of a ST_Intersection(raster, geom) call is a set of geomval
results. Each geometry represents all the neighboring pixels with the
same value. The ST_Intersects function simply returns TRUE/FALSE

About generating new raster layers, we're working on it. Just now on
ST_MapAlgebra function. To output a PostGIS Raster as a different type
of raster, you can use GDAL PostGIS Raster driver, available on GDAL
1.8.0SVN (http://trac.osgeo.org/gdal/wiki/frmts_wtkraster.html)

Actually, you have a version of the driver since GDAL 1.7.x, but I
committed a new version some weeks ago, with memory leaks and bugs
fixed. I know some people are testing the driver, and finding bugs,
and I'd like to fix them all, but I don't have enough time :-( for
everything.


> Yours, S.
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

Best regards,

-- 
Jorge Arévalo
Internet & Mobilty Division, DEIMOS
jorge.arevalo at deimos-space.com
http://mobility.grupodeimos.com/
http://gis4free.wordpress.com