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

[Pierre Racine]

Create a new table with the extents of the raster tiles and check if it overlaps with your points:

CREATE TABLE rastextets AS
SELECT rast::geometry
FROM tbl_raster

Pierre

> -----Original Message-----
> From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-
> bounces at lists.osgeo.org] On Behalf Of Jason Mathis
> Sent: Tuesday, May 13, 2014 5:56 PM
> To: postgis-users at lists.osgeo.org
> Subject: [postgis-users] getting raster values for point from st_value
> 
> 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",SPHER
> OID["GRS
> 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,
> 0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORIT
> Y["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_par
> allel_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],AUT
> HORITY["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.