[postgis-users] PostGIS - Simple Raster Point Query

[Pierre Racine]

ST_Point does not set the SRID you have to do it explicitly:

select st_value(rast,ST_SetSRID(ST_Point(1821416, 5649720), 2193)) from nzdem3

If your raster is big and you want this kind of query to be fast you have to tile your raster when you load it using the -t raster2pgsql option and then use ST_Intersects(rast, yourgeom) in the WHERE clause.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-
> bounces at lists.osgeo.org] On Behalf Of Tim-Hinnerk Heuer
> Sent: Saturday, December 08, 2012 6:29 AM
> To: postgis-users at lists.osgeo.org
> Subject: [postgis-users] PostGIS - Simple Raster Point Query
> 
> Hi,
> 
> I'm new to PostGIS, at least in the backend, so please excuse my naive question,
> but I couldn't figure this out Googling.
> 
> So, I've got it all working (PostGIS 2.0, imported the raster file and got the right
> projection information etc).
> SELECT PostGIS_full_version();
>                                                                 postgis_full_version
> --------------------------------------------------------------------------------------------------
> ---------------------------------------------------
>  POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March
> 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" TOPOLOGY
> RASTER
> (1 row)
> 
> 
> Now I come to querying but cannot seem to get there. I would like to make a
> simple point query for the height in a point on the raster that we have.
> 
> I tried this:
> select st_value(rast,ST_Point(1821416, 5649720)) from nzdem3;
> where ST_Point(1821416, 5649720) is in NZTM (EPSG:2193). However, it doesn't
> work and just complains it is out of range. So, I did some research and found
> that it wants the pixel coordinate. Now, I could calculate it by getting the Geo
> Transformation with GDAL etc, but what is the point of PostGIS if I have to use
> GDAL instead? The point is that I would like to use PostGIS, because it seems like
> it's much more powerful. And we have PostGIS running anyway and it seems like
> the sensible choice.
> 
> I also tried this:
> select DISTINCT ST_World2RasterCoordX(rast,ST_Point(1821416, 5649720))
> from nzdem3;
> 
> but that gives me a pixel offset for each row in the database which really
> confuses me. How am I supposed to query if there are 56 different values for x
> alone:
> st_world2rastercoordx
> -----------------------
>                  -6038
>                  10762
>                  40762
>                  -4838
>                  31162
>                  33562
>                  38362
>                    -38
>                  45562
>                  14362
>                  19162
>                   7162
>                  11962
>                  21562
>                 -13238
>                  20362
>                  46762
>                   1162
>                   8362
>                 -14438
>                  27562
>                  15562
>                  26362
>                 -15638
>                  28762
>                  -9638
>                  -2438
>                 -12038
>                  49162
>                   3562
>                  37162
>                  39562
>                  22762
>                   9562
>                  35962
>                  -3638
>                 -16838
>                   4762
>                  47962
>                   5962
>                  34762
>                   2362
>                  41962
>                  25162
>                  17962
>                  23962
>                  -1238
>                  16762
>                  -7238
>                  32362
>                  13162
>                  43162
>                  -8438
>                 -10838
>                  44362
>                  29962
> (56 rows)
> 
> I see that PostGIS has a lot of potential also with raster queries and ultimately, I
> would like to be able to intersect it with a line string or even a multi poligon. But
> for now I would be happy if I could get a point with the projection's geographic
> coordinates.
> 
> Please help!
> 
> Thanks,
> Tim
> 
> Tim-Hinnerk Heuer
> 
> Twitter: @geekdenz
> 
> Blog: http://www.thheuer.com <http://www.thheuer.com/>