[postgis-devel] How to handle tiled rasters (ST_Value, grid_pt=>coord, coord=>grid_pt etc)?

[Paragon Corporation]

 Stefan,

Hope the following answers your questions.
> First, I've successfully imported into PostGIS 2.0.0 a tiled raster into
table ch1000_tif_tiled, like this:
> raster2pgsql.py" -r ch1000.ti
 >  f -t ch1000_tif_tiled -s 21781 -k 20x20 -I -o ch1000_tif_tiled.sql
> Query 1 - Get no. of bands of a tiled raster:
> SELECT ST_numbands(rast) FROM ch1000_tif_tiled;
  1

>   240 times!
> This is not what I expected; I'd like one row, since I inputted one raster
file (now tiled though). And "SELECT MAX(ST_numbands(rast)) does'nt make it
better. => ?
all the meta data for a raster tile is part of the definition of each tile.
The reason being is that there really isn't anything stopping me from
throwing a bunch of disparate tiles into the same table
much like you can do with geometries.  By convention you don't just like by
convention you don't throw geometries of different types in the same table.

So there are two ways to get at what you want

SELECT DISTINCT ST_NumBands(rast) FROM ch1000_tif_tiled;  (slower if you
have a lot of records but more reliable)

Or if your table was imported via the tools or raster column added via
AddRasterColumn, it would be registered in raster_columns table -- in which
case you can just query that table much like you would query
geometry_columns.

Though I see we don't have a num bands in that table, so you would check for
the array length of the pixel_types field and that would be the number of
bands

SELECT array_upper(pixel_types,1) As num_bands 
  FROM raster_columns WHERE r_table_name = 'ch1000_tif_tiled' AND
r_table_schema = 'public';

----------------------------------------------------------------------------
----------------------------------------------------------------------------
----------------------------------------------

Query 2 - Get value at raster "grid" point 10/20 (default band=1).
> SELECT ST_Value(rast,10,20) FROM ch1000_tif_tiled;
> This query again returns 240 rows (values)! What I'd expect is one value
and one grid point.

The 10,20 represents the 10th column of the 20th row of a given raster tile,
not of the whole coverage (table) .  Again you have to think about each
raster tile as being distinct entities if you query in that fashion.
Think of what would happen if you had an uneven tile coverage and you
inserted a tile in the middle that wasn't there before -- what would you
expect to get when you type in 10,20?


Query 3 - Get projected coordinate and value given a raster "grid"
point at 10/20 (after a note at ST_Value):
> SELECT ST_AsText(ST_SetSRID(ST_Point(ST_UpperLeftX(rast) +
ST_PixelSizeX(rast)*10, ST_UpperLeftY(rast) + ST_PixelSizeY(rast)*20),
ST_SRID(rast))), ST_Value(rast,10,20) 
>    FROM ch1000_tif_tiled;
>  "POINT(489500 102500)";0
>  "POINT(489500 122500)";468
>  "POINT(489500 142500)";1273
>  "POINT(489500 162500)";948
>  "POINT(489500 182500)";0
> => I'd expext only one POINT and it's value.
Same answer as for query 2.

> There are of course more functions, like ST_Envelope(rast) which gives
back 24 envelopes..., etc.
> => Should I really have to care about tiled raster doing MAX, ST_Union,
GROUP BY, etc.?
Yes.  Keep in mind the first objective of PostGIS Raster is to make raster
operations fast.  When you are doing that, you will be working with small
subsections of the raster - a TILE, not the whole raster coverage.
You will want to select subsections to group by etc, not the whole coverage.


Query 4: Given one real projected coordinate (e.g. 489500/102500) return the
corresponding ("closest") grid point and its value.
=> ?  What do you mean by grid point ? Again no concept when you are talking
about a tiled raster at least in PostGIS Raster unles you multiply by the
number of rows in the table etc.

The value would be

SELECT ST_Value(rast, ST_GeomFromText('POINT(489500 102500)', 21781)) 
 FROM ch1000_tif_tiled
WHERE ST_Intersects(rast,  ST_GeomFromText('POINT(489500 102500)', 21781) );



> P.S. BTW, table "raster_overviews" was emtpy after import (as described
above). What's its purpose?
That is expected. I haven't used it myself, but my understanding is
raster_overviews would be filled in if you requested raster overviews in
your import using the -l switch.
It is useful if you want different resolutions of your raster stored.  The
raster_overviews would then list parallel raster tables that hold different
resolutions of the original table.  Kind of like a pyramid.

Mateusz is currently making revisions to that functionality, so he might
have more to say about it.

http://www.postgis.org/pipermail/postgis-devel/2010-December/010836.html


Hope that  helps,
Regina
http://www.postgis.us