[postgis-users] Raster pixel count too high

[Andreas Forø Tollefsen]

Ok. I was worried there was something genious about your suggestion that I
could not grasp :)

Andreas

2013/6/24 Kim Bisgaard <kib at dmi.dk>

>  Hi,
>
> You are right, I mistook ST_intersects() for ST_intersection() - me
> butting out of this thread..
>
> Regards,
> Kim
>
>
> On 2013-06-23 10:46, Andreas Forø Tollefsen wrote:
>
> Thanks for your answers.
>
>
> Kim: We do not directly use ST_Intersection in our script, but ST_Clip to
> clip the raster according to our polygons, and only where they
> ST_Intersects(), not ST_Intersection();
> I still do not understand your suggestion.
>
>  Best,
> Andreas
>
> 2013/6/18 Kim Bisgaard <kib at dmi.dk>
>
>> Hi,
>>
>> Because ST_Intersection() returns neighbour pixels sharing a same value
>> as only one
>>  polygon. You thus get a bigger area and thus more points.
>>
>> Bit me once :-/
>>
>> Regards,
>> Kim
>>
>>
>>
>>
>>
>> On 2013-06-18 14:59, Andreas Forø Tollefsen wrote:
>>
>> Hi Kim,
>>
>>  Thanks for your answer. However, we want raster as an output, since we
>> want to be able to use the summarystats function.
>> Please elaborate how you think ST_PixelAsPolygons should solve out issue?
>> Thanks.
>>
>>  Andreas
>>
>> 2013/6/18 Kim Bisgaard <kib at dmi.dk>
>>
>>> Hi,
>>>
>>> Try to use 'ST_PixelAsPolygons(ST_Clip(n.rast, p.cell))'
>>> instead of 'ST_Intersects(n.rast, p.cell)'
>>>
>>> Regards,
>>> Kim
>>>
>>>
>>> On 2013-06-18 11:03, Andreas Forø Tollefsen wrote:
>>>
>>>  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;
>>>
>>>
>>>
>>>  _______________________________________________
>>> postgis-users mailing listpostgis-users at lists.osgeo.orghttp://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>>
>>>
>>> --
>>> Kim Bisgaard
>>>
>>> Application Development Division     Phone: +45 3915 7562 (direct)
>>> Danish Meteorological Institute      Fax: +45 3915 7460 (division)
>>>
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at lists.osgeo.org
>>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>>
>>>
>>
>>
>> _______________________________________________
>> postgis-users mailing listpostgis-users at lists.osgeo.orghttp://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>
>>
>> --
>> Kim Bisgaard
>>
>> Application Development Division     Phone: +45 3915 7562 (direct)
>> Danish Meteorological Institute      Fax: +45 3915 7460 (division)
>>
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org
>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>
>>
>
>
> _______________________________________________
> postgis-users mailing listpostgis-users at lists.osgeo.orghttp://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>
>
> --
> Kim Bisgaard
>
> Application Development Division     Phone: +45 3915 7562 (direct)
> Danish Meteorological Institute      Fax: +45 3915 7460 (division)
>
>
> _______________________________________________
> 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/20130624/3ef358ba/attachment.html>