[postgis-users] raster intersect with polygon resulting in false NULL value
Tim Keitt
tkeitt at utexas.edu
Wed Aug 21 14:01:35 PDT 2013
Sorry. Should have said QGIS, not postgis. The QGIS stats clearly indicate
the wrong no data value.
THK
On Wed, Aug 21, 2013 at 2:54 PM, Bborie Park <dustymugs at gmail.com> wrote:
> Her table constraints indicate that the NODATA value is properly set to
> -9999.
>
> -bborie
>
>
> On Wed, Aug 21, 2013 at 12:51 PM, Tim Keitt <tkeitt at utexas.edu> wrote:
>
>> 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.
>>
>> THK
>>
>>
>> On Wed, Aug 21, 2013 at 2: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),
>>> '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
>>>
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at lists.osgeo.org
>>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>>
>>>
>>
>>
>> --
>> http://www.keittlab.org/
>>
>> _______________________________________________
>> 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
>
>
--
http://www.keittlab.org/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130821/1aec46f1/attachment.html>
More information about the postgis-users
mailing list