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

Sat Dec 11 09:50:23 PST 2010

Many thanks for the answers!

>> 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.

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?
* 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?
(this fact could certainly give hints to PostGIS to speed up queries).

>> 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.

With "grid point" I mean the same as in query 2 (raster "grid" point 10/20).

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.

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?
* 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

Yours, S.

2010/12/11 Paragon Corporation <lr at pcorp.us>:
> Hmm slight correction  I spoke too soon about the whole grid thing.  I
> think you have to take the whole extent of the table then divide by the
> width and height of the extent by the ST_Width and ST_Height to get numbers
> of tiles in each original grid row.
> Sounds very messy and more than I can think of at the moment.
>
> Again not sure what the utility of that is even when you are rendering on a
> map since you would be rendering tiles that intersect with your area of
> interest using ST_Intersects.
>
>
>
> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net
> [mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of Paragon
> Corporation
> Sent: Saturday, December 11, 2010 11:51 AM
> To: 'PostGIS Users Discussion'; 'Stefan Keller'
> Cc: 'PostGIS Development Discussion'
> Subject: Re: [postgis-users] How to handle tiled rasters
> (ST_Value,grid_pt=>coord, coord=>grid_pt etc)?
>
> Stefan,
>
> Slight addition to what I said in response to this question
>
>
>> 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 rid,  ST_Value(rast, 1, g.pt) As pixvalband1,
> ST_World2RasterX(rast,pt) As x, ST_World2RasterY(rast,pt) As y FROM
> ch1000_tif_tiles INNER JOIN
>        (SELECT ST_GeomFromText('POINT(489500 102500)', 21781)  As pt) As g
>   ON ST_Intersects(rast,  g.pt );
>
> Note the x and y are for that given tile.  If you want the position across
> the whole grid and you have regular gridding and you imported using
> raster2pgsql, then you would use something like
>
> (rid - 1)*ST_Width(rast) + ST_World2RasterX(rast,pt) As x,  (rid -
> 1)*ST_Height(rast) + ST_World2RasterY(rast,pt) As y
>
> Relevant links here
> http://www.postgis.org/documentation/manual-svn/RT_ST_World2RasterCoordX.htm
> l
> http://www.postgis.org/documentation/manual-svn/RT_ST_World2RasterCoordY.htm
> l
> http://www.postgis.org/documentation/manual-svn/RT_ST_Width.html
> http://www.postgis.org/documentation/manual-svn/RT_ST_Height.html
>
> I'm curious though -- which use cases do you have where you need the grid x,
> y from the original raster?
>
>
> Hope that  helps,
> Regina
> http://www.postgis.us
>
>
>
> _______________________________________________
> 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
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>