<div dir="ltr">Gdal is honoring the no-data value encoded in the raster file, but postgis is not. I see this problem in a lot of code out there.<div><br></div><div>THK</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">
On Wed, Aug 21, 2013 at 2:26 PM, Hailey Eckstrand <span dir="ltr"><<a href="mailto:haileyeckstrand@gmail.com" target="_blank">haileyeckstrand@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">Hello,<div>I am trying perform aggregate statistics on a raster vector overlay with a polygon of 192 sub polys. Many of the resulting 192 mean values are correct.. however, I am encountering that one vector overlay mean is NULL which I know to be incorrect. The polygon and raster are in the same projection and when I look at them in QGIS they are overlaid correctly and I can see that the polygon is over top of ~30 raster cells that are all the value 3 (it is not over top of any NA values). One thing that is sort of fishy is that QGIS is reading the raster statistics and no data value incorrectly (values included below) compared to the base ascii raster gdalinfo which is included below.</div>
<div><br></div><div>My query:</div><div><br></div><div><div><font face="courier new, monospace">select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1, geom, true))).mean</font></div><div><font face="courier new, monospace">FROM aet_1,(select * from zben_allwshds where bid = 61) foo</font></div>
<div><font face="courier new, monospace">WHERE ST_Intersects(rast, geom)</font></div><div><font face="courier new, monospace">GROUP BY bid</font></div><div><font face="courier new, monospace">order by bid;</font></div></div>
<div><div><font face="courier new, monospace">NOTICE: The two rasters provided have no intersection. Returning no band raster</font></div><div><font face="courier new, monospace">CONTEXT: PL/pgSQL function "st_clip" line 39 at assignment</font></div>
<div><font face="courier new, monospace">NOTICE: Could not find raster band of index 1 when setting pixel value. Nodata value not set. Returning original raster</font></div><div><font face="courier new, monospace">CONTEXT: PL/pgSQL function "st_clip" line 41 at assignment</font></div>
<div><font face="courier new, monospace">NOTICE: Invalid band index (must use 1-based). Returning NULL</font></div><div><font face="courier new, monospace">CONTEXT: PL/pgSQL function "raster_summarystatsstate" line 9 at assignment</font></div>
<div><font face="courier new, monospace">SQL function "raster_summarystatsstate" statement 1</font></div><div><font face="courier new, monospace"> bid | mean</font></div><div><font face="courier new, monospace">-----+------</font></div>
<div><font face="courier new, monospace"> 61 |</font></div></div><div><br></div><div><br></div><div>Query which returns all the Summary Stats (order of stats is count | sum | mean | stddev | min | max)</div><div><br></div>
<div><div><font face="courier new, monospace">select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1, geom, true)))</font></div><div><font face="courier new, monospace">FROM aet_1,(select * from zben_allwshds where bid = 61) foo</font></div>
<div><font face="courier new, monospace">WHERE ST_Intersects(rast, geom)</font></div><div><font face="courier new, monospace">GROUP BY bid</font></div><div><font face="courier new, monospace">order by bid;</font></div><div>
<font face="courier new, monospace">NOTICE: The two rasters provided have no intersection. Returning no band raster</font></div><div><font face="courier new, monospace">CONTEXT: PL/pgSQL function "st_clip" line 39 at assignment</font></div>
<div><font face="courier new, monospace">NOTICE: Could not find raster band of index 1 when setting pixel value. Nodata value not set. Returning original raster</font></div><div><font face="courier new, monospace">CONTEXT: PL/pgSQL function "st_clip" line 41 at assignment</font></div>
<div><font face="courier new, monospace">NOTICE: Invalid band index (must use 1-based). Returning NULL</font></div><div><font face="courier new, monospace">CONTEXT: PL/pgSQL function "raster_summarystatsstate" line 9 at assignment</font></div>
<div><font face="courier new, monospace">SQL function "raster_summarystatsstate" statement 1</font></div><div><font face="courier new, monospace"> bid | st_summarystatsagg</font></div><div><font face="courier new, monospace">-----+--------------------</font></div>
<div><font face="courier new, monospace"> 61 | (,237,,,3,3)</font></div></div><div><br></div><div><br></div><div><br></div><div>The raster was loaded into postgresql with the raster2pgsql tool:</div><div><br></div><div>
<font face="courier new, monospace">raster2pgsql -s 3005 -I -C -M aet_1.asc -F -t 100x100 public.aet_1 | psql -d NWBC</font><br>
</div><div><br></div><div>gdalinfo of the loaded raster:</div><div><br></div><div><div><font face="courier new, monospace">gdalinfo aet_1.asc</font></div><div><font face="courier new, monospace">Driver: AAIGrid/Arc/Info ASCII Grid</font></div>
<div><font face="courier new, monospace">Files: aet_1.asc</font></div><div><font face="courier new, monospace"> aet_1.asc.aux.xml</font></div><div><font face="courier new, monospace">Size is 2331, 2895</font></div>
<div>
<font face="courier new, monospace">Coordinate System is `'</font></div><div><font face="courier new, monospace">Origin = (173461.904000000009546,1929781.912000000011176)</font></div><div><font face="courier new, monospace">Pixel Size = (400.000000000000000,-400.000000000000000)</font></div>
<div><font face="courier new, monospace">Corner Coordinates:</font></div><div><font face="courier new, monospace">Upper Left ( 173461.904, 1929781.912)</font></div><div><font face="courier new, monospace">Lower Left ( 173461.904, 771781.912)</font></div>
<div><font face="courier new, monospace">Upper Right ( 1105861.904, 1929781.912)</font></div><div><font face="courier new, monospace">Lower Right ( 1105861.904, 771781.912)</font></div><div><font face="courier new, monospace">Center ( 639661.904, 1350781.912)</font></div>
<div><font face="courier new, monospace">Band 1 Block=2331x1 Type=Float32, ColorInterp=Undefined</font></div><div><font face="courier new, monospace"> Min=0.000 Max=9.000</font></div><div><font face="courier new, monospace"> Minimum=0.000, Maximum=9.000, Mean=1.955, StdDev=2.375</font></div>
<div><font face="courier new, monospace"> NoData Value=-9999</font></div><div><font face="courier new, monospace"> Metadata:</font></div><div><font face="courier new, monospace"> STATISTICS_MAXIMUM=9</font></div><div>
<font face="courier new, monospace"> STATISTICS_MEAN=1.9546895049572</font></div><div><font face="courier new, monospace"> STATISTICS_MINIMUM=0</font></div><div><font face="courier new, monospace"> STATISTICS_STDDEV=2.3747620319231</font></div>
</div><div><br></div><div>Incorrect QGIS Raster Statistics:</div><div><br></div><div><p style="margin:0px">STATISTICS_MAXIMUM=1.2341209168929e+033</p>
<p style="margin:12px 0px">STATISTICS_MEAN=5.9293249344197e+031</p>
<p style="margin:12px 0px">STATISTICS_MINIMUM=-9999</p>
<p style="margin:12px 0px">STATISTICS_STDDEV=2.4574481861697e+032</p><p style="margin:12px 0px">NODATA VALUE= -32768 <br></p></div><div><br></div><div>PostGIS raster info:</div><div><br></div><div><div><font face="courier new, monospace"> rid serial NOT NULL,</font></div>
<div><font face="courier new, monospace"> rast raster,</font></div><div><font face="courier new, monospace"> filename text,</font></div><div><font face="courier new, monospace"> CONSTRAINT aet_1_pkey PRIMARY KEY (rid),</font></div>
<div><font face="courier new, monospace"> CONSTRAINT enforce_height_rast CHECK (st_height(rast) = 100),</font></div><div><font face="courier new, monospace"> CONSTRAINT enforce_max_extent_rast CHECK (st_coveredby(st_convexhull(rast), '0103000020BD0B000001000000050000001D5A643BAF2C0541FCA9F1D2EB7D27411D5A643BAF2C0541FED478E935723D41448B6CE7954B3141FED478E935723D41448B6CE7954B3141FCA9F1D2EB7D27411D5A643BAF2C0541FCA9F1D2EB7D2741'::geometry)),</font></div>
<div><font face="courier new, monospace"> CONSTRAINT enforce_nodata_values_rast CHECK (_raster_constraint_nodata_values(rast)::numeric(16,10)[] = '{-9999}'::numeric(16,10)[]),</font></div><div><font face="courier new, monospace"> CONSTRAINT enforce_num_bands_rast CHECK (st_numbands(rast) = 1),</font></div>
<div><font face="courier new, monospace"> CONSTRAINT enforce_out_db_rast CHECK (_raster_constraint_out_db(rast) = '{f}'::boolean[]),</font></div><div><font face="courier new, monospace"> CONSTRAINT enforce_pixel_types_rast CHECK (_raster_constraint_pixel_types(rast) = '{32BF}'::text[]),</font></div>
<div><font face="courier new, monospace"> CONSTRAINT enforce_same_alignment_rast CHECK (st_samealignment(rast, '0100000000000000000000794000000000000079C01D5A643BAF2C0541FED478E935723D4100000000000000000000000000000000BD0B000001000100'::raster)),</font></div>
<div><font face="courier new, monospace"> CONSTRAINT enforce_scalex_rast CHECK (st_scalex(rast)::numeric(16,10) = 400::numeric(16,10)),</font></div><div><font face="courier new, monospace"> CONSTRAINT enforce_scaley_rast CHECK (st_scaley(rast)::numeric(16,10) = (-400)::numeric(16,10)),</font></div>
<div><font face="courier new, monospace"> CONSTRAINT enforce_srid_rast CHECK (st_srid(rast) = 3005),</font></div><div><font face="courier new, monospace"> CONSTRAINT enforce_width_rast CHECK (st_width(rast) = 100)</font></div>
</div><div><br></div><div><br></div><div><br></div><div>Install info:</div><div><br></div><div><font face="courier new, monospace">PostgreSQL Version 9.1</font></div><div><font face="courier new, monospace"><br></font></div>
<div><div><font face="courier new, monospace">postgis_full_version:</font></div><div><font face="courier new, monospace"> POSTGIS="2.0.1 r9979" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.8" TOPOLOGY RASTER</font></div>
<div><font face="courier new, monospace">(1 row)</font></div></div><div><font face="courier new, monospace"><br></font></div><div><div><font face="courier new, monospace">QGIS version:</font></div><div><font face="courier new, monospace">1.8.0-Lisboa</font></div>
<div><font face="courier new, monospace">Compiled against GDAL/OGR</font></div><div><font face="courier new, monospace">1.9.2</font></div><div><font face="courier new, monospace">Running against GDAL/OGR</font></div><div>
<font face="courier new, monospace">1.9.2</font></div><div><font face="courier new, monospace">GEOS Version</font></div><div><font face="courier new, monospace">3.3.8</font></div><div><font face="courier new, monospace">PostgreSQL Client Version</font></div>
<div><font face="courier new, monospace">8.3.10</font></div><div><br></div></div><div><br></div><div>Sorry for all the extra information, I just don't know what people need to know to help. Thanks in advance!</div><span class="HOEnZb"><font color="#888888"><div>
<br></div><div>Hailey</div><div><br></div></font></span></div>
<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" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
<br></blockquote></div><br><br clear="all"><div><br></div>-- <br><div dir="ltr"><a href="http://www.keittlab.org/" target="_blank">http://www.keittlab.org/</a></div>
</div>