<html><head><style>body{font-family:Helvetica,Arial;font-size:13px}</style></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;"><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">Its NAD83/Conus Albers</div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;"><a href="http://epsg.io/5070">http://epsg.io/5070</a></div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;"><br></div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">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. </div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;"><br></div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">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</div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;"><br></div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">1. reproject shapfile</div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">2. create raster files</div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">3. Import files in the database</div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;"><br></div><div id="bloop_customfont" style="font-family:Helvetica,Arial;font-size:13px; color: rgba(0,0,0,1.0); margin: 0px; line-height: auto;">I just reprojected one file using 4269 and it fits in the extents. Although the point is still hanging out there. </div> <div id="bloop_sign_1400083152300542976" class="bloop_sign"><div style="font-family:helvetica,arial;font-size:13px"><br></div></div> I created the point like:<div><br></div><div> SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-111.750185, 34.886948),5070),5070)::geometry(Point,5070)</div><div><br></div><div>Which should be in arizona.</div><div><br></div><div><br></div><div><p style="color:#000;">On May 14, 2014 at 9:52:31 AM, Pierre Racine (<a href="mailto:pierre.racine@sbf.ulaval.ca">pierre.racine@sbf.ulaval.ca</a>) wrote:</p> <blockquote type="cite" class="clean_bq"><span><div><div></div><div>What is this -s 5070 that you put in your raster2pgsql command? Should not it be 4269 (NAD 83)?
<br>
<br>> -----Original Message-----
<br>> From: Jason Mathis [mailto:jmathis@redzonesoftware.com]
<br>> Sent: Wednesday, May 14, 2014 11:47 AM
<br>> To: Pierre Racine; PostGIS Users Discussion
<br>> Subject: RE: [postgis-users] getting raster values for point from st_value
<br>>
<br>> Well I was visualizing the tiff file not the postgis raster table. Its seems
<br>> somewhere along the lines the qgis raster plugin dropped off. Although I am
<br>> able to see the raster tables in the db manager and can drag/drop the table,
<br>> horribly slow but works.
<br>>
<br>>
<br>> The raster extents or just the raster table do not overlay with the point,
<br>> which explains the problem. I still don’t understand what I did wrong.
<br>> Attaching a screen shot.
<br>>
<br>> thanks,
<br>>
<br>>
<br>>
<br>> On May 14, 2014 at 7:01:27 AM, Pierre Racine (pierre.racine@sbf.ulaval.ca)
<br>> wrote:
<br>>
<br>>
<br>> You said you could visualize this table in QGIS? Or were you
<br>> speaking about a shapefile? You can visualize the rastextets table and the
<br>> points table the same way.
<br>>
<br>> Pierre
<br>>
<br>> > -----Original Message-----
<br>> > From: Jason Mathis [mailto:jmathis@redzonesoftware.com]
<br>> > Sent: Tuesday, May 13, 2014 6:30 PM
<br>> > To: PostGIS Users Discussion; Pierre Racine
<br>> > Subject: Re: [postgis-users] getting raster values for point from
<br>> st_value
<br>> >
<br>> > Hi Pierre,
<br>> >
<br>> > Thanks so much for super quick response. You will have to forgive
<br>> me I am
<br>> > fairly new to postgis.
<br>> >
<br>> > Creating the table is no problem but how would I check if it
<br>> overlaps with
<br>> > my points?
<br>> >
<br>> > thanks!
<br>> >
<br>> >
<br>> >
<br>> > On May 13, 2014 at 4:08:40 PM, Pierre Racine
<br>> (pierre.racine@sbf.ulaval.ca)
<br>> > wrote:
<br>> >
<br>> >
<br>> > Create a new table with the extents of the raster tiles and check if
<br>> it
<br>> > overlaps with your points:
<br>> >
<br>> > CREATE TABLE rastextets AS
<br>> > SELECT rast::geometry
<br>> > FROM tbl_raster
<br>> >
<br>> > Pierre
<br>> >
<br>> > > -----Original Message-----
<br>> > > From: postgis-users-bounces@lists.osgeo.org [mailto:postgis-
<br>> > users-
<br>> > > bounces@lists.osgeo.org] On Behalf Of Jason Mathis
<br>> > > Sent: Tuesday, May 13, 2014 5:56 PM
<br>> > > To: postgis-users@lists.osgeo.org
<br>> > > Subject: [postgis-users] getting raster values for point from
<br>> > st_value
<br>> > >
<br>> > > Hi all,
<br>> > >
<br>> > > I have recently loaded a good amount of rasterized shape files
<br>> into
<br>> > a db. I
<br>> > > mainly used gdal_rasterize and raster2pgsql. The shape files
<br>> were
<br>> > in NAD83
<br>> > > so I used the below commands to rasterized and load.
<br>> > >
<br>> > > gdal_rasterize -at -a desc -where "desc='Very High'" -burn 2 -tr
<br>> 30
<br>> > 30 -ot
<br>> > > byte -a_nodata -99 shape_file.shp raster.tiff
<br>> > > raster2pgsql -d -s 5070 -t 100x100 -F -I -C -Y raster.tiff tbl_raster
<br>> |
<br>> > psql -d
<br>> > > raster_db
<br>> > >
<br>> > > I had two sets of shape files one set had about 10 burn values
<br>> and
<br>> > another
<br>> > > had about 50. The issue is getting these values out of the raster.
<br>> I
<br>> > believe I
<br>> > > should be using something from this page:
<br>> > > http://postgis.net/docs/RT_ST_Value.html.
<br>> > >
<br>> > > I can visualize this file in QGIS and it looks good, the histogram
<br>> > shows
<br>> > > values. I can also use the postgis functions, st_summarystats or
<br>> > > st_histogram to see some data but I need to be able to query a
<br>> > given point
<br>> > > and return the raster value.
<br>> > >
<br>> > > I have this:
<br>> > >
<br>> > > SELECT ST_Value(rast,(ST_Point( -111.750185, 34.886948)))
<br>> > > FROM tbl_raster
<br>> > > WHERE st_intersects(rast, (ST_Point(-111.750185, 34.886948)))
<br>> > >
<br>> > > Then i thought oh wait should I set the srid?:
<br>> > >
<br>> > > SELECT ST_Value(rast,(ST_SetSRID(ST_Point( -111.750185,
<br>> > > 34.886948),5070)))
<br>> > > FROM tbl_raster
<br>> > > WHERE st_intersects(rast, (ST_SetSRID(ST_Point(-111.750185,
<br>> > > 34.886948),5070)))
<br>> > >
<br>> > > Still I get nothing, then I thought I need to transform so I got
<br>> this:
<br>> > >
<br>> > > SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -
<br>> > 111.762866,
<br>> > > 34.874309),5070),4326)))
<br>> > > FROM tbl_raster
<br>> > > WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point(
<br>> -
<br>> > > 111.762866, 34.874309),5070),4326)))
<br>> > >
<br>> > > But in the end nothing.
<br>> > >
<br>> > > I did need to add the SRID to the spatial_ref_sys table:
<br>> > >
<br>> > > INSERT into spatial_ref_sys (srid, auth_name, auth_srid,
<br>> > proj4text, srtext)
<br>> > > values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5
<br>> > +lat_0=23
<br>> > > +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80
<br>> > +towgs84=0,0,0,0,0,0,0
<br>> > > +units=m +no_defs ', 'PROJCS["NAD83 / Conus
<br>> > >
<br>> >
<br>> Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHER
<br>> > > OID["GRS
<br>> > >
<br>> >
<br>> 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,
<br>> > >
<br>> >
<br>> 0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORIT
<br>> > >
<br>> >
<br>> Y["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG
<br>> > >
<br>> >
<br>> ","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_
<br>> > >
<br>> >
<br>> Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_par
<br>> > >
<br>> >
<br>> allel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude
<br>> > > _of_center",-
<br>> > >
<br>> >
<br>> 96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["
<br>> > >
<br>> >
<br>> metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUT
<br>> > > HORITY["EPSG","5070"]]’);
<br>> > >
<br>> > > But ultimately I have nothing I can’t seem to get any data out
<br>> > using a query
<br>> > > with a point value. I am thinking it has to do with the srid? Since
<br>> I
<br>> > can
<br>> > > visually inspect it and see data.
<br>> > >
<br>> > > Anyone have any hot ideas as to why I am not seeing data?
<br>> > >
<br>> > > Thanks!
<br>> > >
<br>> > >
<br>> > >
<br>> > > This transmission contains confidential and privileged
<br>> information
<br>> > intended
<br>> > > solely for the party identified above. If you receive this message
<br>> in
<br>> > error,
<br>> > > you must not use it or convey it to others. Please destroy it
<br>> > immediately
<br>> > > and contact the sender at (303) 386-3955 or by return e-mail to
<br>> > the sender.
<br>> >
<br>> > _______________________________________________
<br>> > postgis-users mailing list
<br>> > postgis-users@lists.osgeo.org
<br>> > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
<br>> >
<br>> >
<br>> > This transmission contains confidential and privileged information
<br>> intended
<br>> > solely for the party identified above. If you receive this message in
<br>> error,
<br>> > you must not use it or convey it to others. Please destroy it
<br>> immediately
<br>> > and contact the sender at (303) 386-3955 or by return e-mail to
<br>> the sender.
<br>>
<br>>
<br>>
<br>>
<br>> This transmission contains confidential and privileged information intended
<br>> solely for the party identified above. If you receive this message in error,
<br>> you must not use it or convey it to others. Please destroy it immediately
<br>> and contact the sender at (303) 386-3955 or by return e-mail to the sender.
<br>
<br></div></div></span></blockquote></div></body></html>
<br>
<p><font face="Arial" size="2">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.</font></p>