[postgis-users] getting raster values for point from st_value

[Jason Mathis]

Hi all, 

I have recently loaded a good amount of rasterized shape files into a db. I mainly used gdal_rasterize and raster2pgsql. The shape files were in NAD83 so I used the below commands to rasterized and load. 

gdal_rasterize -at -a desc -where "desc='Very High'" -burn 2 -tr 30 30 -ot byte -a_nodata -99 shape_file.shp raster.tiff
raster2pgsql -d -s 5070 -t 100x100 -F -I -C -Y raster.tiff tbl_raster | psql -d raster_db

I had two sets of shape files one set had about 10 burn values and another had about 50. The issue is getting these values out of the raster. I believe I should be using something from this page: http://postgis.net/docs/RT_ST_Value.html. 

I can visualize this file in QGIS and it looks good, the histogram shows values. I can also use the postgis functions, st_summarystats or st_histogram  to see some data but I need to be able to query a given point and return the raster value. 

I have this:

  SELECT ST_Value(rast,(ST_Point( -111.750185, 34.886948)))
  FROM  tbl_raster
  WHERE st_intersects(rast, (ST_Point(-111.750185, 34.886948)))

Then i thought oh wait should I set the srid?:

  SELECT ST_Value(rast,(ST_SetSRID(ST_Point( -111.750185, 34.886948),5070)))
  FROM  tbl_raster 
  WHERE st_intersects(rast, (ST_SetSRID(ST_Point(-111.750185, 34.886948),5070))) 

Still I get nothing, then I thought I need to transform so I got this:

  SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -111.762866, 34.874309),5070),4326)))
  FROM   tbl_raster 
  WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point( -111.762866, 34.874309),5070),4326)))

But in the end nothing. 

I did need to add the SRID to the spatial_ref_sys table: 

  INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m     +no_defs ', 'PROJCS["NAD83 / Conus Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","5070"]]’);

But ultimately I have nothing I can’t seem to get any data out using a query with a point value. I am thinking it has to do with the srid? Since I can visually inspect it and see data. 

Anyone have any hot ideas as to why I am not seeing data?

Thanks!


-- 
 

This transmission contains confidential and privileged information intended 
solely for the party identified above. If you receive this message in 
error, you must not use it or convey it to others. Please destroy it 
immediately and contact the sender at (303) 386-3955 or by return e-mail to 
the sender.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20140513/5db037de/attachment.html>