2

I've successfully loaded a geoTif into PostGIS (Postgresql 9.3.9) using the following command on my Ubuntu 14.04 server:-

sudo -u postgres raster2pgsql -s 4326 -d -I -C -M -R -l 4 /path/to/my_raster.tif -F -t 10000x10000 my_raster | sudo -u postgres psql -d my_database

It's a 1 banded grayscale raster (lat/lon) with cell values representing some environmental factor on the ground. Hence, I simply want to get the cell value at a given lat/lon.

Following some examples I found I came up with the following SQL:-

SELECT
  ST_Value( rast, ST_SetSRID(ST_MakePoint(0,51.5),4326) ) AS rast_val
FROM 
  my_raster
WHERE
  ST_Value( rast, ST_SetSRID(ST_MakePoint(0,51.5),4326) ) IS NOT NULL;

This seems to get the right value for me when I cross reference the same location loading the raster in QGIS.

My question is basically a sanity check if the above SQL is correct as I am not 100% sure of the WHERE clause.

2
  • It looks correct. You don't really need the WHERE clause, though I imagine the query parser will convert the query into something that doesn't actually call ST_Value twice. You could probably handle that more elegantly using exclude_nodata_value true/false that is built into the function. Commented Mar 4, 2016 at 15:13
  • Thanks, but "exclude_nodata_value" does not seem make any difference either it is set to TRUE of FALSE, unless I am doing it the wrong way: 'SELECT ST_Value( rast, ST_SetSRID(ST_MakePoint(0,51.5),4326), TRUE ) AS rast_val FROM my_raster' or 'SELECT ST_Value( rast, ST_SetSRID(ST_MakePoint(0,51.5),4326), FALSE ) AS rast_val FROM my_raster'
    – Jason
    Commented Mar 7, 2016 at 11:09

1 Answer 1

2

A faster way is to cut your raster into small tiles while loading it with the -t option (e.g. 100x100). Make sure also to build an index on the tiles with the -I option.

Then you can query the pixel values with a query like this taking advantage of the indexed tiles:

SELECT ST_Value(rast, ST_SetSRID(ST_MakePoint(0,51.5),4326)) AS val
FROM my_raster
WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(0,51.5),4326));

Should be way faster....

1
  • Thanks for the tips. Turning on \timing, I can see your query is at least twice as fast as my original query.
    – Jason
    Commented Mar 8, 2016 at 10:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.