[postgis-users] raster intersect with polygon resulting in false NULL value
Pierre Racine
Pierre.Racine at sbf.ulaval.ca
Wed Aug 21 13:45:51 PDT 2013
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
>
>
>
More information about the postgis-users
mailing list