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

[Jorge Arévalo]

On Fri, Nov 26, 2010 at 11:08 AM, Stefan Keller <sfkeller at gmail.com> wrote:
>> 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.
>

Currently, there's no such function. We now are able to intersect a
raster and a vector, and obtain a vector as output. We do it by
polygonizing the raster first, and the intersecting the polygonized
raster with the vector.

The intersection of a raster and a vector to get a raster as output is
planned, but not implemented yet.

Best regards,


>
> 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
>>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>



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