[postgis-users] raster intersect with polygon resulting in false NULL value

[Hailey Eckstrand]

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),
'0103000020BD0B000001000000050000001D5A643BAF2C0541FCA9F1D2EB7D27411D5A643BAF2C0541FED478E935723D41448B6CE7954B3141FED478E935723D41448B6CE7954B3141FCA9F1D2EB7D27411D5A643BAF2C0541FCA9F1D2EB7D2741'::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,
'0100000000000000000000794000000000000079C01D5A643BAF2C0541FED478E935723D4100000000000000000000000000000000BD0B000001000100'::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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130821/621a86d3/attachment.html>