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

[Paragon Corporation]

> Ok; I realise now that I assumed the special case "rectangular regularly
tiled raster coverage" (defined as arrangement d. in
> http://trac.osgeo.org/postgis/wiki/WKTRaster/Documentation01 ).

> But still given I have imported one single raster "rectangular regularly
tiled raster coverage"
> * Don't you think this and my related queries above could be a typical and
use case?

Its typical that you would have a regular raster coverage, but I'm still not
clear why knowing the original grid column/row from the original raster is
relevant for any use case.
To me its like trying to remember the original row a geometry was in in the
ESRI shape file I imported. 

If you can state a specific example that would be great :)  

I'm pretty green as far as raster analysis goes so please don't take my lack
of understanding  as an insutl.

> * This probably deserves simplified queries where the user does'nt have to
think about speed and implementation?
> * After doing a raster2pgsql (.py) with tiling option "-k 20x20" does the
fact get lost, that this originally was a single raster file?

The memory of it being a single raster is kept in the raster_columns table.
You can use the additional file switchin import  which would include the
original file name as a column in the raster table.  This of course is most
Useful if you are importing a whole folder of rasters presumably
hierarchically ordered.

 Perhaps Pierre, Jorge, or Mat have more thoughts on the matter.

> (this fact could certainly give hints to PostGIS to speed up queries).
How does it give hints?  Right now each raster tile has a spatial bounding
box around it.  So that is wha tit uses to speed up queries.  The only more
in depth hint I can think of is when getting into the actual
Pixels of a tile (by generating some histogram of the distribution of pixel
ranges to make pixel value queries much much faster), then again -- knowing
that a tile is part of a bigger set I still don't see the usefulness of as
far as speed is concerned.

The assumption I am making of course is that your raster coverage is so big
that you would never ask for the whole thing. And if you were planning to do
it that way, then you would never chop up your raster file  in the first
place
So you would only have one record in your table where the columns and rows
map directly to your original file.



> So I have to revise e.g. my query 4 (and your solution):
> Query 4: Given one real projected coordinate (e.g. 489500/102500) return
the raster objects (or tiles) and the ("closest") raster grid point(s) and
> value(s) therein.

> Thanks for your solution proposal:

> SELECT rid, ST_Value(rast, ST_GeomFromText('POINT(489500 102500)', 21781))
>  FROM ch1000_tif_tiled
>  WHERE ST_Intersects(rast, ST_GeomFromText('POINT(489500 102500)', 21781)
);
1: 181;
2: 201;0

> The result would be ok for the general case (not knowing how many
different rasters are in the table).

> But the resuls seems weird (having two rows) given that I imported only
and only one raster as tiles.

Oh you are getting 2 records back.  This should only happen if your point
happens to lie at the edge of 2 pixels.

I see your point though that if your raster were in a single record, you
would get one record back and the record would in fact be arbitrary of the
2.
I think Pierre had discussed this a while back what to do in this scenario.


Really in my mind, the correct thing to do would be to return an array.  But
correctness is less useful in many cases than just arbitrarily picking one.
Especially
Since this SHOULD BE a rare occurrence that a point lies at the edge of two
pixels.

Though the fact they are sitting in two different rows means the point is a
the edge of a tile and a pixel which would seem even rarer unless you
specifically
Picked up a point in that location.  Is that the case?
If not then this is a bug.

> Regarding metadata I played around with the examples in chapter 13 using
an untiled raster 'relief_jpg' with 3 bands:

> -- Listing 13.1 Getting general summary info about rasters in a kauai
table SELECT count(*) As num_rasters, ST_Height(rast) As height,
    ST_Width(rast) As width, ST_SRID(rast) As srid,
    ST_NumBands(rast) As num_bands,
    ST_BandPixelType(rast,1) As btype
FROM relief_jpg
GROUP BY ST_Height(rast) ,
    ST_Width(rast), ST_SRID(rast),
    ST_NumBands(rast),
    ST_BandPixelType(rast,1);
-- 1;584;891;21781;3;"8BUI"

> -- Solution using metadata
SELECT num_rasters, blocksize_x as height, blocksize_y as width, srid,
array_length(pixel_types, 1) AS num_bands, pixel_types as btypes
  FROM
    raster_columns as r,
    (SELECT count(*) AS num_rasters FROM relief_jpg) dummy
  WHERE r_table_name = 'relief_jpg';
-- 1;;;21781;3;"{8BUI,8BUI,8BUI}"

>  So, two things are unclear to me:
> * btypes should be an array in your solution since there are 3 bands,
aren't they?
Yes they are which the text representation of an array in PostgreSQL is
{8BUI,8BUI,8BUI}

> * in my solution there is a mix of table name as element and as text
string.
> * Obviously raster2pgsql did not set blocksize_x/blocksize_y

Don't quite understand what you mean by mix of table name as element and as
text string.
>From your query it looks like you only have one raster record in you
relief_jpg.
I think the raster2pgsql routine only sets the blocksize_x and y if you
choose to tile the raster.  I suppose it might make 
Sense to just have it set to the width and height of the raster if the
raster is not tiled.  Then again the fact its not set 
Gives you the clue that its tiled verses not tiled, which seems kind of
useful.

Thanks for reviewing the chapter,
Regina

http://www.postgis.us