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

[Stefan Keller]

> the result of a ST_Intersection(raster, geom) call is a set of geomval
> results.

Ok; geomval seems to be a predefined struct which contains first the
geometry, second the associated value.
And ST_Intersection() currently polygonizes the raster to vector.
But how is the query specified (if yet) in order to get a raster as
the result of calling ST_Intersection(), i.e. rasterizing first the
vector geom component before doing the overlay with raster?

Yours, S.


2010/11/25 Jorge Arévalo <jorge.arevalo at deimos-space.com>:
> 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
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>