<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=utf-8"><meta name=Generator content="Microsoft Word 14 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:Helvetica;
panose-1:2 11 6 4 2 2 2 2 2 4;}
@font-face
{font-family:Helvetica;
panose-1:2 11 6 4 2 2 2 2 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p
{mso-style-priority:99;
mso-margin-top-alt:auto;
margin-right:0cm;
mso-margin-bottom-alt:auto;
margin-left:0cm;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
span.EmailStyle18
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:72.0pt 72.0pt 72.0pt 72.0pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-GB link=blue vlink=purple><div class=WordSection1><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Why rasterize polygons at all in this case, if you can run much simpler queries directly on the polygons? <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Rasterizing polygons is actually a waste of disk space. All rasterized pixels will supposedly have the same (single?) attribute as the polygon itself.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>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<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>(e.g. where desc = “Very High”, etc.). Reprojection works on polygons as well, of course.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>GL<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><div><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Guido Lemoine<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>European Commission, Joint Research Centre<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Institute for Environment and Sustainability, Monitoring Agricultural Resources Unit<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>TP 266 – 21027 Ispra, Italy, T. +39 0332 786239, F. +39 0332 785162<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>E-mail: <a href="mailto:guido.lemoine@jrc.ec.europa.eu"><span style='color:blue'>guido.lemoine@jrc.ec.europa.eu</span></a> URL: <a href="http://ec.europa.eu/jrc"><span style='color:blue'>http://ec.europa.eu/jrc</span></a><o:p></o:p></span></p></div><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><div><div style='border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm'><p class=MsoNormal><b><span lang=EN-US style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>From:</span></b><span lang=EN-US style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> postgis-users-bounces@lists.osgeo.org [mailto:postgis-users-bounces@lists.osgeo.org] <b>On Behalf Of </b>Jason Mathis<br><b>Sent:</b> 14 May 2014 00:30<br><b>To:</b> PostGIS Users Discussion; Pierre Racine<br><b>Subject:</b> Re: [postgis-users] getting raster values for point from st_value<o:p></o:p></span></p></div></div><p class=MsoNormal><o:p> </o:p></p><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'>Hi Pierre,<o:p></o:p></span></p></div><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'><o:p> </o:p></span></p></div><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'>Thanks so much for super quick response. You will have to forgive me I am fairly new to postgis. <o:p></o:p></span></p></div><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'><o:p> </o:p></span></p></div><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'>Creating the table is no problem but how would I check if it overlaps with my points? <o:p></o:p></span></p></div><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'><o:p> </o:p></span></p></div><div id="bloop_customfont"><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'>thanks!<o:p></o:p></span></p></div><div id="bloop_sign_1400020101137347072"><div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'><o:p> </o:p></span></p></div></div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'><o:p> </o:p></span></p><p><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif";color:black'>On May 13, 2014 at 4:08:40 PM, Pierre Racine (<a href="mailto:pierre.racine@sbf.ulaval.ca">pierre.racine@sbf.ulaval.ca</a>) wrote:<o:p></o:p></span></p><blockquote style='margin-top:5.0pt;margin-bottom:5.0pt'><div><div><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Helvetica","sans-serif"'>Create a new table with the extents of the raster tiles and check if it 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: <a href="mailto:postgis-users-bounces@lists.osgeo.org">postgis-users-bounces@lists.osgeo.org</a> [mailto:postgis-users- <br>> <a href="mailto:bounces@lists.osgeo.org">bounces@lists.osgeo.org</a>] On Behalf Of Jason Mathis <br>> Sent: Tuesday, May 13, 2014 5:56 PM <br>> To: <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a> <br>> Subject: [postgis-users] getting raster values for point from st_value <br>> <br>> Hi all, <br>> <br>> I have recently loaded a good amount of rasterized shape files into a db. I <br>> mainly used gdal_rasterize and raster2pgsql. The shape files were 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 30 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 | psql -d <br>> raster_db <br>> <br>> I had two sets of shape files one set had about 10 burn values and another <br>> had about 50. The issue is getting these values out of the raster. I believe I <br>> should be using something from this page: <br>> <a href="http://postgis.net/docs/RT_ST_Value.html">http://postgis.net/docs/RT_ST_Value.html</a>. <br>> <br>> I can visualize this file in QGIS and it looks good, the histogram 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 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 this: <br>> <br>> SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -111.762866, <br>> 34.874309),5070),4326))) <br>> FROM tbl_raster <br>> WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point( - <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, proj4text, srtext) <br>> values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 <br>> +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 <br>> +units=m +no_defs ', 'PROJCS["NAD83 / Conus <br>> Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHER <br>> OID["GRS <br>> 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0, <br>> 0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORIT <br>> Y["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG <br>> ","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_ <br>> Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_par <br>> allel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude <br>> _of_center",- <br>> 96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT[" <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 using a query <br>> with a point value. I am thinking it has to do with the srid? Since I 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 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>_______________________________________________<br>postgis-users mailing list<br><a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br><a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><o:p></o:p></span></p></div></div></blockquote><p class=MsoNormal><o:p> </o:p></p><p><span style='font-size:10.0pt;font-family:"Arial","sans-serif"'>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.</span><o:p></o:p></p></div></body></html>