[postgis-users] Raster pixel count too high

[Andreas Forø Tollefsen]

Hi,

We are working on a raster summarystats script to calculate various
statistics for the pixels within fishnet polygons.

Our raster cell size is 0.0083333333333... x 0.0083333333333... degrees
while our quadrat polygons are 0.5 x 0.5 decimal degrees.
This should give us 60x60 raster pixels within each of our polygons. ArcGIS
zonal statistics returns a pixel count of 3600 in addition to other
statistics.
However, PostGIS returns 3721 pixel count.

We do not really understand why, but it seems that our query includes some
pixels that are outside of the polygon, but still touches the vertices of
the polygon and are therefore included in the calculation.
Are there any way of modifying our script to return the same result as
ArcGIS?
Thanks!

Andreas

script:

/* This query makes one raster for each PRIO-GRID cell. Clip and union is
the procedure. */
INSERT INTO nightlightsprio (gid, "year", rast)
(SELECT gid, "year", ST_Union(raster) as rast
FROM
(SELECT p.gid, n."year", ST_Clip(n.rast, p.cell) as raster
FROM nightlights n, priogridyear p
WHERE ST_Intersects(n.rast, p.cell)
AND n."year" = p."year"
)
as priorast
GROUP BY gid, "year");


/* Default BandNoDataValue is 0. Raster value 0 means no light, not no
data. Setting to NULL. This produces correct results. */
UPDATE nightlightsprio2 SET rast = ST_SetBandNoDataValue(rast, 1, NULL);


ALTER TABLE nightlightsprio2 ADD COLUMN nightlights_sum double precision,
ADD COLUMN nightlights_mean double precision,
ADD COLUMN nightlights_sd double precision,
ADD COLUMN nightlights_min double precision,
ADD COLUMN nightlights_max double precision,
ADD COLUMN nightlights_count integer;

UPDATE nightlightsprio2 SET nightlights_sum = (ST_SummaryStats(rast)).sum;
UPDATE nightlightsprio2 SET nightlights_mean = (ST_SummaryStats(rast)).mean;
UPDATE nightlightsprio2 SET nightlights_sd = (ST_SummaryStats(rast)).stddev;
UPDATE nightlightsprio2 SET nightlights_min = (ST_SummaryStats(rast)).min;
UPDATE nightlightsprio2 SET nightlights_max = (ST_SummaryStats(rast)).max;
UPDATE nightlightsprio2 SET nightlights_count =
(ST_SummaryStats(rast)).count;
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130618/0b00b405/attachment.html>