<div dir="ltr"><div>Hi Pierre,</div><div>That did result in the correct response: </div><div><br></div><div><font face="courier new, monospace">NWBC=# select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1, geom, true)))</font></div>
<div><font face="courier new, monospace">NWBC-# FROM aet_1,(select * from zben_allwshds where bid = 61) foo</font></div><div><font face="courier new, monospace">NWBC-# WHERE ST_Intersects(rast, geom)</font></div><div><font face="courier new, monospace">NWBC-# GROUP BY bid</font></div>
<div><font face="courier new, monospace">NWBC-# 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><br></div><div><br></div><div>Did I have an old version of the function?</div><div>Thanks!</div><div>Hailey</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Aug 21, 2013 at 2:34 PM, Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca" target="_blank">Pierre.Racine@sbf.ulaval.ca</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Heiley,<br>
<br>
What if you execute this and then retry your query:<br>
<br>
CREATE OR REPLACE FUNCTION raster_summarystatsstate(ss summarystats, rast raster, nband int DEFAULT 1, exclude_nodata_value boolean DEFAULT TRUE, sample_percent double precision DEFAULT 1)<br>
    RETURNS summarystats<br>
    AS $$<br>
    DECLARE<br>
        newstats summarystats;<br>
        ret summarystats;<br>
    BEGIN<br>
        IF rast IS NULL OR ST_HasNoBand(rast) OR ST_IsEmpty(rast) THEN<br>
            RETURN ss;<br>
        END IF;<br>
        newstats := _ST_SummaryStats(rast, nband, exclude_nodata_value, sample_percent);<br>
        IF $1 IS NULL THEN<br>
            ret := (newstats.count,<br>
                    newstats.sum,<br>
                    null,<br>
                    null,<br>
                    newstats.min,<br>
                    newstats.max)::summarystats;<br>
        ELSE<br>
            ret := (COALESCE(ss.count,0) + COALESCE(newstats.count, 0),<br>
                    COALESCE(ss.sum,0) + COALESCE(newstats.sum, 0),<br>
                    null,<br>
                    null,<br>
                    least(ss.min, newstats.min),<br>
                    greatest(ss.max, newstats.max))::summarystats;<br>
        END IF;<br>
        RETURN ret;<br>
    END;<br>
    $$<br>
    LANGUAGE 'plpgsql';<br>
<div class="im HOEnZb"><br>
> -----Original Message-----<br>
> From: Hailey Eckstrand [mailto:<a href="mailto:haileyeckstrand@gmail.com">haileyeckstrand@gmail.com</a>]<br>
> Sent: Wednesday, August 21, 2013 5:05 PM<br>
> To: PostGIS Users Discussion<br>
</div><div class="im HOEnZb">> Cc: <a href="mailto:dustymugs@gmail.com">dustymugs@gmail.com</a>; Pierre Racine<br>
> Subject: Re: [postgis-users] raster intersect with polygon resulting in false<br>
> NULL value<br>
><br>
</div><div class="HOEnZb"><div class="h5">> Hi Pierre,<br>
> I just downloaded & installed this morning:<br>
> <a href="http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_su" target="_blank">http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_su</a><br>
> mmarystatsagg.sql<br>
><br>
><br>
> But here is the print out:<br>
> <a href="http://hastebin.com/tubelisoda.pas" target="_blank">http://hastebin.com/tubelisoda.pas</a><br>
><br>
><br>
> Hailey<br>
><br>
><br>
> On Wed, Aug 21, 2013 at 1:45 PM, Pierre Racine<br>
> <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>> wrote:<br>
><br>
><br>
>       Hailey,<br>
><br>
>       Could you post here the definition of ST_SummaryStatsAgg() you<br>
> are using. I'm pretty sure I fixed this some time ago. I'd like to compare<br>
> what you have with what I have.<br>
><br>
>       Pierre<br>
><br>
><br>
>       > -----Original Message-----<br>
>       > From: <a href="mailto:postgis-users-bounces@lists.osgeo.org">postgis-users-bounces@lists.osgeo.org</a> [mailto:<a href="mailto:postgis-">postgis-</a><br>
> users-<br>
>       > <a href="mailto:bounces@lists.osgeo.org">bounces@lists.osgeo.org</a>] On Behalf Of Hailey Eckstrand<br>
>       > Sent: Wednesday, August 21, 2013 4:18 PM<br>
>       > To: PostGIS Users Discussion; <a href="mailto:dustymugs@gmail.com">dustymugs@gmail.com</a><br>
>       > Subject: Re: [postgis-users] raster intersect with polygon resulting<br>
> in false<br>
>       > NULL value<br>
>       ><br>
>       > It looks like you were right bborie..<br>
>       ><br>
>       > WITH foo AS (<br>
>       > SELECT * FROM zben_allwshds WHERE bid = 61<br>
>       > )<br>
>       > SELECT<br>
>       > a.rid,<br>
>       > ST_SummaryStats(ST_Clip(a.rast,1, f.geom, true))<br>
>       > FROM aet_1 a<br>
>       > JOIN foo f<br>
>       > ON ST_Intersects(a.rast, f.geom);<br>
>       > NOTICE:  The two rasters provided have no intersection.<br>
> Returning no band<br>
>       > raster<br>
>       > CONTEXT:  PL/pgSQL function "st_clip" line 39 at assignment<br>
>       > NOTICE:  Could not find raster band of index 1 when setting pixel<br>
> value.<br>
>       > Nodata value not set. Returning original raster<br>
>       > CONTEXT:  PL/pgSQL function "st_clip" line 41 at assignment<br>
>       > NOTICE:  Invalid band index (must use 1-based). Returning NULL<br>
>       > CONTEXT:  SQL function "st_summarystats" statement 1<br>
>       >  rid | st_summarystats<br>
>       > -----+------------------<br>
>       >  122 | (79,237,3,0,3,3)<br>
>       >  146 |<br>
>       > (2 rows)<br>
>       ><br>
>       ><br>
>       ><br>
>       > On Wed, Aug 21, 2013 at 1:05 PM, Bborie Park<br>
> <<a href="mailto:dustymugs@gmail.com">dustymugs@gmail.com</a>><br>
>       > wrote:<br>
>       ><br>
>       ><br>
>       >       Hailey,<br>
>       ><br>
>       >       Can you run a query such as the following and post the<br>
> results? I<br>
>       > had to do a double-take for the ST_SummaryStatsAgg() function<br>
> as I'm not<br>
>       > familiar with that function.<br>
>       ><br>
>       >       WITH foo AS (<br>
>       >       SELECT * FROM zben_allwshds WHERE bid = 61<br>
>       >       )<br>
>       >       SELECT<br>
>       >       a.rid,<br>
>       >       ST_SummaryStats(ST_Clip(a.rast,1, f.geom, true))<br>
>       >       FROM aet_1 a<br>
>       >       JOIN foo f<br>
>       >       ON ST_Intersects(a.rast, f.geom)<br>
>       ><br>
>       >       I'm suspecting that ST_SummaryStatsAgg() is having an issue<br>
> where<br>
>       > one tile's summary stats has NULL values being combined with<br>
> the<br>
>       > aggregate summary stats.<br>
>       ><br>
>       ><br>
>       >       -bborie<br>
>       ><br>
>       ><br>
>       ><br>
>       >       On Wed, Aug 21, 2013 at 12:26 PM, Hailey Eckstrand<br>
>       > <<a href="mailto:haileyeckstrand@gmail.com">haileyeckstrand@gmail.com</a>> wrote:<br>
>       ><br>
>       ><br>
>       >               Hello,<br>
>       >               I am trying perform aggregate statistics on a raster vector<br>
>       > overlay with a polygon of 192 sub polys. Many of the resulting<br>
> 192 mean<br>
>       > values are correct.. however, I am encountering that one vector<br>
> overlay<br>
>       > mean is NULL which I know to be incorrect. The polygon and<br>
> raster are in<br>
>       > the same projection and when I look at them in QGIS they are<br>
> overlaid<br>
>       > correctly and I can see that the polygon is over top of ~30 raster<br>
> cells that<br>
>       > are all the value 3 (it is not over top of any NA values). One thing<br>
> that is sort<br>
>       > of fishy is that QGIS is reading the raster statistics and no data<br>
> value<br>
>       > incorrectly (values included below) compared to the base ascii<br>
> raster<br>
>       > gdalinfo which is included below.<br>
>       ><br>
>       >               My query:<br>
>       ><br>
>       >               select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1, geom,<br>
>       > true))).mean<br>
>       >               FROM aet_1,(select * from zben_allwshds where bid = 61)<br>
>       > foo<br>
>       >               WHERE ST_Intersects(rast, geom)<br>
>       >               GROUP BY bid<br>
>       >               order by bid;<br>
>       >               NOTICE:  The two rasters provided have no intersection.<br>
>       > Returning no band raster<br>
>       >               CONTEXT:  PL/pgSQL function "st_clip" line 39 at<br>
>       > assignment<br>
>       >               NOTICE:  Could not find raster band of index 1 when<br>
> setting<br>
>       > pixel value. Nodata value not set. Returning original raster<br>
>       >               CONTEXT:  PL/pgSQL function "st_clip" line 41 at<br>
>       > assignment<br>
>       >               NOTICE:  Invalid band index (must use 1-based). Returning<br>
>       > NULL<br>
>       >               CONTEXT:  PL/pgSQL function<br>
> "raster_summarystatsstate"<br>
>       > line 9 at assignment<br>
>       >               SQL function "raster_summarystatsstate" statement 1<br>
>       >                bid | mean<br>
>       >               -----+------<br>
>       >                 61 |<br>
>       ><br>
>       ><br>
>       >               Query which returns all the Summary Stats (order of stats<br>
> is<br>
>       > count | sum | mean | stddev | min | max)<br>
>       ><br>
>       >               select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1, geom,<br>
>       > true)))<br>
>       >               FROM aet_1,(select * from zben_allwshds where bid = 61)<br>
>       > foo<br>
>       >               WHERE ST_Intersects(rast, geom)<br>
>       >               GROUP BY bid<br>
>       >               order by bid;<br>
>       >               NOTICE:  The two rasters provided have no intersection.<br>
>       > Returning no band raster<br>
>       >               CONTEXT:  PL/pgSQL function "st_clip" line 39 at<br>
>       > assignment<br>
>       >               NOTICE:  Could not find raster band of index 1 when<br>
> setting<br>
>       > pixel value. Nodata value not set. Returning original raster<br>
>       >               CONTEXT:  PL/pgSQL function "st_clip" line 41 at<br>
>       > assignment<br>
>       >               NOTICE:  Invalid band index (must use 1-based). Returning<br>
>       > NULL<br>
>       >               CONTEXT:  PL/pgSQL function<br>
> "raster_summarystatsstate"<br>
>       > line 9 at assignment<br>
>       >               SQL function "raster_summarystatsstate" statement 1<br>
>       >                bid | st_summarystatsagg<br>
>       >               -----+--------------------<br>
>       >                 61 | (,237,,,3,3)<br>
>       ><br>
>       ><br>
>       ><br>
>       >               The raster was loaded into postgresql with the<br>
> raster2pgsql<br>
>       > tool:<br>
>       ><br>
>       >               raster2pgsql -s 3005 -I -C -M aet_1.asc -F -t 100x100<br>
>       > public.aet_1 | psql -d NWBC<br>
>       ><br>
>       ><br>
>       >               gdalinfo of the loaded raster:<br>
>       ><br>
>       >               gdalinfo aet_1.asc<br>
>       >               Driver: AAIGrid/Arc/Info ASCII Grid<br>
>       >               Files: aet_1.asc<br>
>       >                      aet_1.asc.aux.xml<br>
>       >               Size is 2331, 2895<br>
>       >               Coordinate System is `'<br>
>       >               Origin =<br>
>       > (173461.904000000009546,1929781.912000000011176)<br>
>       >               Pixel Size = (400.000000000000000,-<br>
>       > 400.000000000000000)<br>
>       >               Corner Coordinates:<br>
>       >               Upper Left  (  173461.904, 1929781.912)<br>
>       >               Lower Left  (  173461.904,  771781.912)<br>
>       >               Upper Right ( 1105861.904, 1929781.912)<br>
>       >               Lower Right ( 1105861.904,  771781.912)<br>
>       >               Center      (  639661.904, 1350781.912)<br>
>       >               Band 1 Block=2331x1 Type=Float32,<br>
> ColorInterp=Undefined<br>
>       >                 Min=0.000 Max=9.000<br>
>       >                 Minimum=0.000, Maximum=9.000, Mean=1.955,<br>
>       > StdDev=2.375<br>
>       >                 NoData Value=-9999<br>
>       >                 Metadata:<br>
>       >                   STATISTICS_MAXIMUM=9<br>
>       >                   STATISTICS_MEAN=1.9546895049572<br>
>       >                   STATISTICS_MINIMUM=0<br>
>       >                   STATISTICS_STDDEV=2.3747620319231<br>
>       ><br>
>       >               Incorrect QGIS Raster Statistics:<br>
>       ><br>
>       ><br>
>       >               STATISTICS_MAXIMUM=1.2341209168929e+033<br>
>       ><br>
>       >               STATISTICS_MEAN=5.9293249344197e+031<br>
>       ><br>
>       >               STATISTICS_MINIMUM=-9999<br>
>       ><br>
>       >               STATISTICS_STDDEV=2.4574481861697e+032<br>
>       ><br>
>       >               NODATA VALUE= -32768<br>
>       ><br>
>       ><br>
>       ><br>
>       >               PostGIS raster info:<br>
>       ><br>
>       >                 rid serial NOT NULL,<br>
>       >                 rast raster,<br>
>       >                 filename text,<br>
>       >                 CONSTRAINT aet_1_pkey PRIMARY KEY (rid),<br>
>       >                 CONSTRAINT enforce_height_rast CHECK (st_height(rast)<br>
> =<br>
>       > 100),<br>
>       >                 CONSTRAINT enforce_max_extent_rast CHECK<br>
>       > (st_coveredby(st_convexhull(rast),<br>
>       ><br>
> '0103000020BD0B000001000000050000001D5A643BAF2C0541FCA9F1D2<br>
>       ><br>
> EB7D27411D5A643BAF2C0541FED478E935723D41448B6CE7954B3141FE<br>
>       ><br>
> D478E935723D41448B6CE7954B3141FCA9F1D2EB7D27411D5A643BAF2C<br>
>       > 0541FCA9F1D2EB7D2741'::geometry)),<br>
>       >                 CONSTRAINT enforce_nodata_values_rast CHECK<br>
>       > (_raster_constraint_nodata_values(rast)::numeric(16,10)[] = '{-<br>
>       > 9999}'::numeric(16,10)[]),<br>
>       >                 CONSTRAINT enforce_num_bands_rast CHECK<br>
>       > (st_numbands(rast) = 1),<br>
>       >                 CONSTRAINT enforce_out_db_rast CHECK<br>
>       > (_raster_constraint_out_db(rast) = '{f}'::boolean[]),<br>
>       >                 CONSTRAINT enforce_pixel_types_rast CHECK<br>
>       > (_raster_constraint_pixel_types(rast) = '{32BF}'::text[]),<br>
>       >                 CONSTRAINT enforce_same_alignment_rast CHECK<br>
>       > (st_samealignment(rast,<br>
>       ><br>
> '0100000000000000000000794000000000000079C01D5A643BAF2C0541<br>
>       ><br>
> FED478E935723D4100000000000000000000000000000000BD0B000001<br>
>       > 000100'::raster)),<br>
>       >                 CONSTRAINT enforce_scalex_rast CHECK<br>
>       > (st_scalex(rast)::numeric(16,10) = 400::numeric(16,10)),<br>
>       >                 CONSTRAINT enforce_scaley_rast CHECK<br>
>       > (st_scaley(rast)::numeric(16,10) = (-400)::numeric(16,10)),<br>
>       >                 CONSTRAINT enforce_srid_rast CHECK (st_srid(rast) =<br>
>       > 3005),<br>
>       >                 CONSTRAINT enforce_width_rast CHECK (st_width(rast)<br>
> =<br>
>       > 100)<br>
>       ><br>
>       ><br>
>       ><br>
>       >               Install info:<br>
>       ><br>
>       >               PostgreSQL Version 9.1<br>
>       ><br>
>       ><br>
>       >               postgis_full_version:<br>
>       >                POSTGIS="2.0.1 r9979" GEOS="3.3.8-CAPI-1.7.8"<br>
>       > PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released<br>
> 2012/10/08"<br>
>       > LIBXML="2.7.8" TOPOLOGY RASTER<br>
>       >               (1 row)<br>
>       ><br>
>       ><br>
>       >               QGIS version:<br>
>       >               1.8.0-Lisboa<br>
>       >               Compiled against GDAL/OGR<br>
>       >               1.9.2<br>
>       >               Running against GDAL/OGR<br>
>       >               1.9.2<br>
>       >               GEOS Version<br>
>       >               3.3.8<br>
>       >               PostgreSQL Client Version<br>
>       >               8.3.10<br>
>       ><br>
>       ><br>
>       >               Sorry for all the extra information, I just don't know what<br>
>       > people need to know to help. Thanks in advance!<br>
>       ><br>
>       ><br>
>       >               Hailey<br>
>       ><br>
>       ><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-" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-</a><br>
>       > users<br>
>       ><br>
>       ><br>
>       ><br>
>       ><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" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
>       ><br>
>       ><br>
>       ><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" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>