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

Pierre Racine Pierre.Racine at sbf.ulaval.ca
Wed May 14 10:09:34 PDT 2014 

Make sure your layers align properly before importing them in the DB and don't forget to answer Guido's question...

> -----Original Message-----
> From: Jason Mathis [mailto:jmathis at redzonesoftware.com]
> Sent: Wednesday, May 14, 2014 12:35 PM
> To: Pierre Racine; PostGIS Users Discussion
> Subject: RE: [postgis-users] getting raster values for point from st_value
> 
> Its NAD83/Conus Albers
> http://epsg.io/5070
> 
> I needed to reproject the shape file as well. It was originally in an Albers
> projection but before handed over to me was projected into wgs84. The
> raster and shape files did not line up properly until we reprojected the
> shape file first and then did the conversion to raster.
> 
> I had some issues with the reprojection as well and someone on the gdal
> mailing list suggested using the epsg:5070 for the reprojection. So the
> whole process really went
> 
> 1. reproject shapfile
> 2. create raster files
> 3. Import files in the database
> 
> I just reprojected one file using 4269 and it fits in the extents. Although the
> point is still hanging out there.
> 
> I created the point like:
> 
>   SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-111.750185,
> 34.886948),5070),5070)::geometry(Point,5070)
> 
> Which should be in arizona.
> 
> 
> 
> On May 14, 2014 at 9:52:31 AM, Pierre Racine (pierre.racine at sbf.ulaval.ca)
> wrote:
> 
> 
> 	What is this -s 5070 that you put in your raster2pgsql command?
> Should not it be 4269 (NAD 83)?
> 
> 	> -----Original Message-----
> 	> From: Jason Mathis [mailto:jmathis at redzonesoftware.com]
> 	> Sent: Wednesday, May 14, 2014 11:47 AM
> 	> To: Pierre Racine; PostGIS Users Discussion
> 	> Subject: RE: [postgis-users] getting raster values for point from
> st_value
> 	>
> 	> Well I was visualizing the tiff file not the postgis raster table. Its
> seems
> 	> somewhere along the lines the qgis raster plugin dropped off.
> Although I am
> 	> able to see the raster tables in the db manager and can drag/drop
> the table,
> 	> horribly slow but works.
> 	>
> 	>
> 	> The raster extents or just the raster table do not overlay with the
> point,
> 	> which explains the problem. I still don’t understand what I did
> wrong.
> 	> Attaching a screen shot.
> 	>
> 	> thanks,
> 	>
> 	>
> 	>
> 	> On May 14, 2014 at 7:01:27 AM, Pierre Racine
> (pierre.racine at sbf.ulaval.ca)
> 	> wrote:
> 	>
> 	>
> 	> You said you could visualize this table in QGIS? Or were you
> 	> speaking about a shapefile? You can visualize the rastextets table
> and the
> 	> points table the same way.
> 	>
> 	> Pierre
> 	>
> 	> > -----Original Message-----
> 	> > From: Jason Mathis [mailto:jmathis at redzonesoftware.com]
> 	> > Sent: Tuesday, May 13, 2014 6:30 PM
> 	> > 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.
> 	>
> 	>
> 	>
> 	>
> 	> 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.
> 
> 
> 
> 
> 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.

More information about the postgis-users mailing list