[postgis-users] raster intersect with polygon resulting in false NULL value
Hailey Eckstrand
haileyeckstrand at gmail.com
Wed Aug 21 14:43:38 PDT 2013
Hi Pierre,
That did result in the correct response:
NWBC=# select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1, geom, true)))
NWBC-# FROM aet_1,(select * from zben_allwshds where bid = 61) foo
NWBC-# WHERE ST_Intersects(rast, geom)
NWBC-# GROUP BY bid
NWBC-# order by bid;
NOTICE: The two rasters provided have no intersection. Returning no band
raster
CONTEXT: PL/pgSQL function "st_clip" line 39 at assignment
NOTICE: Could not find raster band of index 1 when setting pixel value.
Nodata value not set. Returning original raster
CONTEXT: PL/pgSQL function "st_clip" line 41 at assignment
Did I have an old version of the function?
Thanks!
Hailey
On Wed, Aug 21, 2013 at 2:34 PM, Pierre Racine
<Pierre.Racine at sbf.ulaval.ca>wrote:
> Heiley,
>
> What if you execute this and then retry your query:
>
> 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)
> RETURNS summarystats
> AS $$
> DECLARE
> newstats summarystats;
> ret summarystats;
> BEGIN
> IF rast IS NULL OR ST_HasNoBand(rast) OR ST_IsEmpty(rast) THEN
> RETURN ss;
> END IF;
> newstats := _ST_SummaryStats(rast, nband, exclude_nodata_value,
> sample_percent);
> IF $1 IS NULL THEN
> ret := (newstats.count,
> newstats.sum,
> null,
> null,
> newstats.min,
> newstats.max)::summarystats;
> ELSE
> ret := (COALESCE(ss.count,0) + COALESCE(newstats.count, 0),
> COALESCE(ss.sum,0) + COALESCE(newstats.sum, 0),
> null,
> null,
> least(ss.min, newstats.min),
> greatest(ss.max, newstats.max))::summarystats;
> END IF;
> RETURN ret;
> END;
> $$
> LANGUAGE 'plpgsql';
>
> > -----Original Message-----
> > From: Hailey Eckstrand [mailto:haileyeckstrand at gmail.com]
> > Sent: Wednesday, August 21, 2013 5:05 PM
> > To: PostGIS Users Discussion
> > Cc: dustymugs at gmail.com; Pierre Racine
> > Subject: Re: [postgis-users] raster intersect with polygon resulting in
> false
> > NULL value
> >
> > Hi Pierre,
> > I just downloaded & installed this morning:
> > http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_su
> > mmarystatsagg.sql
> >
> >
> > But here is the print out:
> > http://hastebin.com/tubelisoda.pas
> >
> >
> > Hailey
> >
> >
> > On Wed, Aug 21, 2013 at 1:45 PM, Pierre Racine
> > <Pierre.Racine at sbf.ulaval.ca> wrote:
> >
> >
> > Hailey,
> >
> > Could you post here the definition of ST_SummaryStatsAgg() you
> > are using. I'm pretty sure I fixed this some time ago. I'd like to
> compare
> > what you have with what I have.
> >
> > Pierre
> >
> >
> > > -----Original Message-----
> > > From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-
> > users-
> > > bounces at lists.osgeo.org] On Behalf Of Hailey Eckstrand
> > > Sent: Wednesday, August 21, 2013 4:18 PM
> > > To: PostGIS Users Discussion; dustymugs at gmail.com
> > > Subject: Re: [postgis-users] raster intersect with polygon
> resulting
> > in false
> > > NULL value
> > >
> > > It looks like you were right bborie..
> > >
> > > WITH foo AS (
> > > SELECT * FROM zben_allwshds WHERE bid = 61
> > > )
> > > SELECT
> > > a.rid,
> > > ST_SummaryStats(ST_Clip(a.rast,1, f.geom, true))
> > > FROM aet_1 a
> > > JOIN foo f
> > > ON ST_Intersects(a.rast, f.geom);
> > > NOTICE: The two rasters provided have no intersection.
> > Returning no band
> > > raster
> > > CONTEXT: PL/pgSQL function "st_clip" line 39 at assignment
> > > NOTICE: Could not find raster band of index 1 when setting pixel
> > value.
> > > Nodata value not set. Returning original raster
> > > CONTEXT: PL/pgSQL function "st_clip" line 41 at assignment
> > > NOTICE: Invalid band index (must use 1-based). Returning NULL
> > > CONTEXT: SQL function "st_summarystats" statement 1
> > > rid | st_summarystats
> > > -----+------------------
> > > 122 | (79,237,3,0,3,3)
> > > 146 |
> > > (2 rows)
> > >
> > >
> > >
> > > On Wed, Aug 21, 2013 at 1:05 PM, Bborie Park
> > <dustymugs at gmail.com>
> > > wrote:
> > >
> > >
> > > Hailey,
> > >
> > > Can you run a query such as the following and post the
> > results? I
> > > had to do a double-take for the ST_SummaryStatsAgg() function
> > as I'm not
> > > familiar with that function.
> > >
> > > WITH foo AS (
> > > SELECT * FROM zben_allwshds WHERE bid = 61
> > > )
> > > SELECT
> > > a.rid,
> > > ST_SummaryStats(ST_Clip(a.rast,1, f.geom, true))
> > > FROM aet_1 a
> > > JOIN foo f
> > > ON ST_Intersects(a.rast, f.geom)
> > >
> > > I'm suspecting that ST_SummaryStatsAgg() is having an issue
> > where
> > > one tile's summary stats has NULL values being combined with
> > the
> > > aggregate summary stats.
> > >
> > >
> > > -bborie
> > >
> > >
> > >
> > > On Wed, Aug 21, 2013 at 12:26 PM, Hailey Eckstrand
> > > <haileyeckstrand at gmail.com> wrote:
> > >
> > >
> > > Hello,
> > > 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.
> > >
> > > My query:
> > >
> > > select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1,
> geom,
> > > true))).mean
> > > FROM aet_1,(select * from zben_allwshds where bid
> = 61)
> > > foo
> > > WHERE ST_Intersects(rast, geom)
> > > GROUP BY bid
> > > order by bid;
> > > NOTICE: The two rasters provided have no
> intersection.
> > > Returning no band raster
> > > CONTEXT: PL/pgSQL function "st_clip" line 39 at
> > > assignment
> > > NOTICE: Could not find raster band of index 1 when
> > setting
> > > pixel value. Nodata value not set. Returning original raster
> > > CONTEXT: PL/pgSQL function "st_clip" line 41 at
> > > assignment
> > > NOTICE: Invalid band index (must use 1-based).
> Returning
> > > NULL
> > > CONTEXT: PL/pgSQL function
> > "raster_summarystatsstate"
> > > line 9 at assignment
> > > SQL function "raster_summarystatsstate" statement 1
> > > bid | mean
> > > -----+------
> > > 61 |
> > >
> > >
> > > Query which returns all the Summary Stats (order
> of stats
> > is
> > > count | sum | mean | stddev | min | max)
> > >
> > > select bid,(ST_SummaryStatsAgg(ST_Clip(rast,1,
> geom,
> > > true)))
> > > FROM aet_1,(select * from zben_allwshds where bid
> = 61)
> > > foo
> > > WHERE ST_Intersects(rast, geom)
> > > GROUP BY bid
> > > order by bid;
> > > NOTICE: The two rasters provided have no
> intersection.
> > > Returning no band raster
> > > CONTEXT: PL/pgSQL function "st_clip" line 39 at
> > > assignment
> > > NOTICE: Could not find raster band of index 1 when
> > setting
> > > pixel value. Nodata value not set. Returning original raster
> > > CONTEXT: PL/pgSQL function "st_clip" line 41 at
> > > assignment
> > > NOTICE: Invalid band index (must use 1-based).
> Returning
> > > NULL
> > > CONTEXT: PL/pgSQL function
> > "raster_summarystatsstate"
> > > line 9 at assignment
> > > SQL function "raster_summarystatsstate" statement 1
> > > bid | st_summarystatsagg
> > > -----+--------------------
> > > 61 | (,237,,,3,3)
> > >
> > >
> > >
> > > The raster was loaded into postgresql with the
> > raster2pgsql
> > > tool:
> > >
> > > raster2pgsql -s 3005 -I -C -M aet_1.asc -F -t
> 100x100
> > > public.aet_1 | psql -d NWBC
> > >
> > >
> > > gdalinfo of the loaded raster:
> > >
> > > gdalinfo aet_1.asc
> > > Driver: AAIGrid/Arc/Info ASCII Grid
> > > Files: aet_1.asc
> > > aet_1.asc.aux.xml
> > > Size is 2331, 2895
> > > Coordinate System is `'
> > > Origin =
> > > (173461.904000000009546,1929781.912000000011176)
> > > Pixel Size = (400.000000000000000,-
> > > 400.000000000000000)
> > > Corner Coordinates:
> > > Upper Left ( 173461.904, 1929781.912)
> > > Lower Left ( 173461.904, 771781.912)
> > > Upper Right ( 1105861.904, 1929781.912)
> > > Lower Right ( 1105861.904, 771781.912)
> > > Center ( 639661.904, 1350781.912)
> > > Band 1 Block=2331x1 Type=Float32,
> > ColorInterp=Undefined
> > > Min=0.000 Max=9.000
> > > Minimum=0.000, Maximum=9.000, Mean=1.955,
> > > StdDev=2.375
> > > NoData Value=-9999
> > > Metadata:
> > > STATISTICS_MAXIMUM=9
> > > STATISTICS_MEAN=1.9546895049572
> > > STATISTICS_MINIMUM=0
> > > STATISTICS_STDDEV=2.3747620319231
> > >
> > > Incorrect QGIS Raster Statistics:
> > >
> > >
> > > STATISTICS_MAXIMUM=1.2341209168929e+033
> > >
> > > STATISTICS_MEAN=5.9293249344197e+031
> > >
> > > STATISTICS_MINIMUM=-9999
> > >
> > > STATISTICS_STDDEV=2.4574481861697e+032
> > >
> > > NODATA VALUE= -32768
> > >
> > >
> > >
> > > PostGIS raster info:
> > >
> > > rid serial NOT NULL,
> > > rast raster,
> > > filename text,
> > > CONSTRAINT aet_1_pkey PRIMARY KEY (rid),
> > > CONSTRAINT enforce_height_rast CHECK
> (st_height(rast)
> > =
> > > 100),
> > > CONSTRAINT enforce_max_extent_rast CHECK
> > > (st_coveredby(st_convexhull(rast),
> > >
> > '0103000020BD0B000001000000050000001D5A643BAF2C0541FCA9F1D2
> > >
> > EB7D27411D5A643BAF2C0541FED478E935723D41448B6CE7954B3141FE
> > >
> > D478E935723D41448B6CE7954B3141FCA9F1D2EB7D27411D5A643BAF2C
> > > 0541FCA9F1D2EB7D2741'::geometry)),
> > > CONSTRAINT enforce_nodata_values_rast CHECK
> > > (_raster_constraint_nodata_values(rast)::numeric(16,10)[] = '{-
> > > 9999}'::numeric(16,10)[]),
> > > CONSTRAINT enforce_num_bands_rast CHECK
> > > (st_numbands(rast) = 1),
> > > CONSTRAINT enforce_out_db_rast CHECK
> > > (_raster_constraint_out_db(rast) = '{f}'::boolean[]),
> > > CONSTRAINT enforce_pixel_types_rast CHECK
> > > (_raster_constraint_pixel_types(rast) = '{32BF}'::text[]),
> > > CONSTRAINT enforce_same_alignment_rast CHECK
> > > (st_samealignment(rast,
> > >
> > '0100000000000000000000794000000000000079C01D5A643BAF2C0541
> > >
> > FED478E935723D4100000000000000000000000000000000BD0B000001
> > > 000100'::raster)),
> > > CONSTRAINT enforce_scalex_rast CHECK
> > > (st_scalex(rast)::numeric(16,10) = 400::numeric(16,10)),
> > > CONSTRAINT enforce_scaley_rast CHECK
> > > (st_scaley(rast)::numeric(16,10) = (-400)::numeric(16,10)),
> > > CONSTRAINT enforce_srid_rast CHECK
> (st_srid(rast) =
> > > 3005),
> > > CONSTRAINT enforce_width_rast CHECK
> (st_width(rast)
> > =
> > > 100)
> > >
> > >
> > >
> > > Install info:
> > >
> > > PostgreSQL Version 9.1
> > >
> > >
> > > postgis_full_version:
> > > 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
> > > (1 row)
> > >
> > >
> > > QGIS version:
> > > 1.8.0-Lisboa
> > > Compiled against GDAL/OGR
> > > 1.9.2
> > > Running against GDAL/OGR
> > > 1.9.2
> > > GEOS Version
> > > 3.3.8
> > > PostgreSQL Client Version
> > > 8.3.10
> > >
> > >
> > > Sorry for all the extra information, I just don't
> know what
> > > people need to know to help. Thanks in advance!
> > >
> > >
> > > Hailey
> > >
> > >
> > >
> > > _______________________________________________
> > > postgis-users mailing list
> > > postgis-users at lists.osgeo.org
> > >
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-
> > > users
> > >
> > >
> > >
> > >
> > >
> > > _______________________________________________
> > > postgis-users mailing list
> > > postgis-users at lists.osgeo.org
> > >
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> > >
> > >
> > >
> >
> > _______________________________________________
> > postgis-users mailing list
> > postgis-users at lists.osgeo.org
> > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130821/e4011885/attachment.html>
More information about the postgis-users
mailing list