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

Thu May 15 00:24:03 PDT 2014

Hey,
if perfo is important,
you may want to use the version with many points at the same time :
http://postgis.net/docs/RT_ST_SetValues.html.

Cheers,
Rémi-C

2014-05-15 1:05 GMT+02:00 Jason Mathis <jmathis at redzonesoftware.com>:

> 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.
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20140515/201a6da0/attachment.html>