[postgis-users] Raster pixel value
Andreas Forø Tollefsen
andreasft at gmail.com
Mon Dec 19 02:13:36 PST 2011
Line 59 to 77 in the st_summarystatsagg.sql.
Could it be that line 69: ($1).sum / ($1).count leads to division by zero
error if ($1).count is 0?
How could i change this so this is not the case?
-- raster_summarystatsfinal
-- Final function used by the ST_SummaryStatsAgg aggregate
CREATE OR REPLACE FUNCTION raster_summarystatsfinal(ss summarystats)
RETURNS summarystats
AS $$
DECLARE
ret summarystats;
BEGIN
ret := (($1).count,
($1).sum,
($1).sum / ($1).count,
null,
($1).min,
($1).max
)::summarystats;
RETURN ret;
END;
$$
LANGUAGE 'plpgsql';
2011/12/19 Andreas Forø Tollefsen <andreasft at gmail.com>
> Updated to latest trunk, and now it works.
> However, I ran into a new "division by zero" error. This halts the query.
> Could it be that something is divided by integer rather than decimal in
> the query, or what else causes this:
>
> ERROR: division by zero
> CONTEXT: PL/pgSQL function "raster_summarystatsfinal" line 5 at assignment
>
> ********** Error **********
>
> ERROR: division by zero
> SQL state: 22012
> Context: PL/pgSQL function "raster_summarystatsfinal" line 5 at assignment
>
> Best,
> Andreas
>
>
> 2011/12/19 Tom van Tilburg <tom.van.tilburg at gmail.com>
>
>> Andreas,
>>
>> If you got the latest ST_Clip from the repository, I found this
>> information with it:
>>
>> "Addition of C-based ST_MinPossibleValue to replace the existing
>> ST_MinPossibleVal which uses hard-coded values. Updated dependent functions
>> and scripts/plpgsql to use new function. Deleted
>> scripts/plpgsql/st_minpossibleval.sql to stop people from using it.
>> Associated ticket is #1298 <http://trac.osgeo.org/postgis/ticket/1298>."
>>
>> Probably this means you need the very newest version of rtpostgis with
>> the minpossiblevalue included in C. Not sure wheter is available for
>> windows yet.
>> At the moment I don't have time to test but I have the feeling this new
>> function in C might make things quicker and perhaps fix the error with
>> nodata values I mentioned earlier. Let me know.
>>
>> Regards,
>> Tom
>>
>>
>> On 19-12-2011 10:01, Andreas Forø Tollefsen wrote:
>>
>> Hi Tom,
>>
>> I tried both functions. The St_AreaWeightedSummaryStats() works great,
>> but it takes a lot of time to complete.
>>
>> However, the ST_Clip() i cannot manage to get working. Seems like
>> something is broken. Any idea what can cause this error?
>> I installed the latest ST_SummaryStatsAgg() and ST_Clip() from
>> http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql
>>
>> The query:
>> SELECT foo2.gid,
>> (ss).count,
>> (ss).sum,
>> (ss).mean,
>> (ss).min,
>> (ss).max
>> FROM
>> (SELECT foo.gid, ST_SummaryStatsAgg(gv) ss
>> FROM (SELECT p.gid, ST_Clip(r.rast, p.cell) gv
>> FROM access r, priogrid_land p
>> WHERE ST_Intersects(r.rast, p.cell)
>> ) foo
>> GROUP BY foo.gid
>> ) foo2
>>
>> Results in:
>>
>> ERROR: function st_minpossiblevalue(text) does not exist
>> LINE 1: ...esce(nodata, ST_BandNodataValue(rast, bandstart), ST_MinPoss...
>> ^
>> HINT: No function matches the given name and argument types. You might
>> need to add explicit type casts.
>> QUERY: SELECT coalesce(nodata, ST_BandNodataValue(rast, bandstart),
>> ST_MinPossibleValue(newpixtype))
>> CONTEXT: PL/pgSQL function "st_clip" line 31 at assignment
>>
>> ********** Error **********
>>
>> ERROR: function st_minpossiblevalue(text) does not exist
>> SQL state: 42883
>> Hint: No function matches the given name and argument types. You might
>> need to add explicit type casts.
>> Context: PL/pgSQL function "st_clip" line 31 at assignment
>>
>>
>>
>> 2011/12/6 Tom van Tilburg <tom.van.tilburg at gmail.com>
>>
>>> Andreas,
>>>
>>> I didn't have time to reproduce your problem yet. Did you have any
>>> succes by yourself on this issue?
>>> Could it have something to do with counting the non-data values as
>>> value? This is what I experience with a similar function (ST_Clip) that
>>> consequently gave me the value '0' instead of the nodata value. The result
>>> of that is that the mean was often lower than expected.
>>>
>>> Perhaps you could rewrite your previous example to something with
>>> auto-generated values in the script. That saves time in reproducing.
>>>
>>> Cheers,
>>> Tom
>>>
>>>
>>>
>>> On 25-11-2011 18:11, Andreas Forø Tollefsen wrote:
>>>
>>> Update:
>>>
>>> I think my suspicion is correct. If I do a ST_Summarystats().sum and
>>> divide this on 36 my MAX value will be 1.
>>> Hence, I think the number of values counted and the number of
>>> observations counted is not equal.
>>>
>>> New query:
>>> DROP TABLE IF EXISTS mountain_phil_cell;
>>>
>>> SELECT
>>> a.gid As gid,
>>> (ST_SummaryStats((ST_Union(ST_MapAlgebraExpr(ST_AsRaster(a.cell, b.rast,
>>> '32BF'), b.rast, 'rast2', '32BF','INTERSECTION','0','0',0))).rast,
>>> false)).sum / 36 As avgmnt
>>> INTO mountain_phil_cell
>>> FROM
>>> priogrid_land a LEFT JOIN
>>> mountain_phil b
>>> ON ST_Intersects(a.cell, b.rast)
>>> GROUP BY a.gid
>>> ORDER BY a.gid;
>>>
>>>
>>> 2011/11/25 Andreas Forø Tollefsen <andreasft at gmail.com>
>>>
>>>> A small note regarding this issue.
>>>>
>>>> My problem is that I never get a mean value of 1 even if all pixels
>>>> inside the geometry is one.
>>>>
>>>> Could this be because: 6x6 pixels goes into one polygon when visually
>>>> controlling. If each pixel has the value 1, then this will be calculated as
>>>> 36 / 36 = 1. However, if it calculates the sum to be 36 and divide by a
>>>> number higher than 36 pixels, then the result will always be below 1.
>>>> What i am thinking is that while it sums up the pixel values correctly,
>>>> it does not count only the 36 pixels, but also neighboring pixels.
>>>> Therefore: 1+1+1+1...n36 / Number of pixels higher than 36 will always lead
>>>> to a number lower than 1.
>>>>
>>>> Anyone who knows the functions well could probably answer this.
>>>>
>>>> Best regards,
>>>> Andreas
>>>>
>>>>
>>>> 2011/11/25 Andreas Forø Tollefsen <andreasft at gmail.com>
>>>>
>>>>> Could this have to do with the tiling of the raster?
>>>>> I will try to run the same query with a untiled mountain raster to see
>>>>> if that changes anything.
>>>>>
>>>>> Btw. When loading a tiled postgis raster into qgis it shows up with
>>>>> many artifacts and no data areas. The same raster untiled does not show up
>>>>> the same way.
>>>>> Qgis bug?
>>>>>
>>>>> Andreas
>>>>>
>>>>> 2011/11/25 Andreas Forø Tollefsen <andreasft at gmail.com>
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> Thanks for all of the suggestions. I will do some more testing.
>>>>>> However, as for suggestion 1 i think the pixel size should be the same as
>>>>>> the original raster or am I wrong?
>>>>>>
>>>>>> Both the mean_mnt_bin raster and the priogrid_land shapefile can be
>>>>>> downloaded as zip (2 mb) here:
>>>>>> http://gisintersect.com/mean_mnt_bin.zip
>>>>>> http://gisintersect.com/priogrid_land.zip
>>>>>>
>>>>>> Any help on getting the correct values would be very much
>>>>>> appreciated.
>>>>>>
>>>>>> My query:
>>>>>> DROP TABLE IF EXISTS mountain_cell;
>>>>>>
>>>>>> SELECT
>>>>>> a.gid As id,
>>>>>> (ST_SummaryStats((ST_Union(ST_MapAlgebraExpr(ST_AsRaster(a.cell,
>>>>>> b.rast, '32BF'), b.rast, 'rast2', '32BF','INTERSECTION','0','0',0))).rast,
>>>>>> false)).mean As avgmnt
>>>>>> INTO mountain_cell
>>>>>> FROM
>>>>>> priogrid_land a LEFT JOIN
>>>>>> mountain b
>>>>>> ON ST_Intersects(a.cell, b.rast)
>>>>>> GROUP BY a.gid
>>>>>> ORDER BY a.gid;
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>> _______________________________________________
>>> postgis-users mailing listpostgis-users at postgis.refractions.nethttp://postgis.refractions.net/mailman/listinfo/postgis-users
>>>
>>>
>>>
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at postgis.refractions.net
>>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>>
>>>
>>
>>
>> _______________________________________________
>> postgis-users mailing listpostgis-users at postgis.refractions.nethttp://postgis.refractions.net/mailman/listinfo/postgis-users
>>
>>
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20111219/6fa7a6d3/attachment.html>
More information about the postgis-users
mailing list