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

Wed May 14 16:05:27 PDT 2014

So I figured this out. It was my lack of understanding the functions, projections, and coordinates. Seems all the data is correct but my queries were not. My final query that works is:

  SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_MakePoint(-111.574508,35.296170),4326),5070)))
  FROM  tbl_raster 
  WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_MakePoint(-111.574508,35.296170),4326),5070))) 

I needed to first set the srid in 4326 AND then transform to 5070. Pretty simple solution! 

Thanks,
jason 

On May 14, 2014 at 1:46:04 PM, Jason Mathis (jmathis at redzonesoftware.com) wrote:

Hi Guido,

The example I gave in my question is just the start of how we want to analyze the data. Giving a point and getting back the value I thought would be fairly easy to do, so I started there. I believe some other calculations maybe more complicated and better suited to raster. One example, given a line we would need to calculate all the values along that line to a given point. I can’t go into too much detail because I fairly new to this world, and are being told to do things from much smarter people. 

Also this data was originally raster data. Its actually LandFire data if you are familiar. We partnered with another company who vectorized the data probably for similar reasons you stated. They also reprojected the data before giving it to us. Sooooo hence the process of reproject and rasterize back to the original formats. 

Hope that answers your question. We got about 11GB of shape files after the rasterize and loading it was only a little over 1GB of data in the db; so space was not too much of a concern. 

thanks for all the help hopefully I can figure this out:) 

-jason 

On May 14, 2014 at 12:29:46 AM, Guido Lemoine (guido.lemoine at jrc.ec.europa.eu) wrote:

Why rasterize polygons at all in this case, if you can run much simpler queries directly on the polygons?

Rasterizing polygons is actually a waste of disk space. All rasterized pixels will supposedly have the same (single?) attribute as the polygon itself.

Try using shp2pgsql first to check if you can directly work on the polygons. The attribute(s) will appear as (a) separate column(s), so you can include them in your predicates

(e.g. where desc = “Very High”, etc.). Reprojection works on polygons as well, of course.

GL

Guido Lemoine

European Commission, Joint Research Centre

Institute for Environment and Sustainability, Monitoring Agricultural Resources Unit

TP 266 – 21027 Ispra, Italy, T. +39 0332 786239, F. +39 0332 785162

E-mail: guido.lemoine at jrc.ec.europa.eu  URL: http://ec.europa.eu/jrc

From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-bounces at lists.osgeo.org] On Behalf Of Jason Mathis
Sent: 14 May 2014 00:30
To: PostGIS Users Discussion; Pierre Racine
Subject: Re: [postgis-users] getting raster values for point from st_value

Hi Pierre,

Thanks so much for super quick response. You will have to forgive me I am fairly new to postgis. 

Creating the table is no problem but how would I check if it overlaps with my points? 

thanks!

On May 13, 2014 at 4:08:40 PM, Pierre Racine (pierre.racine at sbf.ulaval.ca) wrote:

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.

_______________________________________________
postgis-users mailing list
postgis-users at lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

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.

_______________________________________________
postgis-users mailing list
postgis-users at lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
-- 

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/20140514/be5a9f13/attachment.html>