[postgis-users] Raster pixel value

Andreas Forø Tollefsen andreasft at gmail.com
Thu Nov 24 04:42:23 PST 2011


I double checked the values that I get from this query and they do not seem
to be correct.
For polygon cells that have only pixels with the value 1 inside only have a
mean of 0.85 (see image).
The correct value should be a mean of 1. After some visual controls I am
not sure if these numbers are correct.
There are many polygon cells with only pixel value 1, but no polygon cell
with a higher mean than .87.

Any idea why this could be? Both the raster and polygons are in srid 4326.

Andreas



2011/11/24 Andreas Forø Tollefsen <andreasft at gmail.com>

> I know. RTFM before asking dumb questions. Set the ST_Summary(rast, false)
> so it does not exclude_nodata_value? Right?
>
> Andreas
>
>
> 2011/11/24 Andreas Forø Tollefsen <andreasft at gmail.com>
>
>> Great. I updated to latest rev and added the st_union.sql to my function
>> list. Now it works.
>> However, I need to ask one more question.
>> When running the script, it calculates the average pixel value inside the
>> polygon. It does not take into account the NULL DATA values.
>> Is there any way i can make the script calculate the average pixel value
>> but count NULL DATA values as 0?
>>
>> For instance, one pixel in a cell belonging to the Comoros have a value
>> (0.95). The rest is NULL. Hence, it seems a bit wrong that the whole
>> polygon get 0.95 when almost none of its area have any elevation?
>>
>> Any ideas?
>>
>> Thank you so much for this info Tom!
>>
>> Andreas
>>
>> 2011/11/24 Tom van Tilburg <tom.van.tilburg at gmail.com>
>>
>>>  Yes, I think so.
>>> Current windows build is *r8221
>>> *
>>> Also, I think ST_Union(raster) (which I used in the example) is not
>>> included in this version yet.
>>> You would have to download a prototype from.
>>>
>>> http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_union.sql
>>>
>>> Chrs,
>>>  Tom
>>>
>>>
>>> On 24-11-2011 12:07, Andreas Forø Tollefsen wrote:
>>>
>>> Great, Thank you so much for this.
>>> However, I do not seem to have the ST_MapAlgebraExpr(rast, rast.....)
>>> function. Was this implemented in rev 8001?
>>>
>>>  Andreas
>>>
>>> 2011/11/24 Tom van Tilburg <tom.van.tilburg at gmail.com>
>>>
>>>>  Andreas,
>>>>
>>>> We did approx. the same thing for non-quadrate polygons. Perhaps it
>>>> might be useful:
>>>>
>>>> Step 1: Make a raster from every polygon, based on the grid
>>>> specifications of the elevation raster. Here is also the solution: the
>>>> raster cells will only be created for cells that have their midpoint
>>>> *inside* your geometry.
>>>>     ST_AsRaster(a.geom, b.rast, '<pixtype>')
>>>>
>>>> Step 2: Overlay the elevation raster with the raster you just created
>>>> and keep only the values of the elevation raster
>>>> ST_MapAlgebraExpr(
>>>>             <raster from geometry>
>>>>             ,b.rast
>>>>             ,'rast2' -- <-- keep only raster 2 value
>>>>             , '<pixtype>','INTERSECTION','0','0',0
>>>>         )
>>>>
>>>> Step 3: get the mean from the statistics on the resulting raster
>>>> (ST_SummaryStats(
>>>>     (ST_Union( -- < --- we did a UNION because we occasionaly had
>>>> vectors crossing tiled rasters
>>>>         <overlay raster>
>>>>     )).rast
>>>> )).mean As avg_height
>>>>
>>>> That did the trick. Complete script is below.
>>>>
>>>> I suspect your method of doing a ST_Intersection for every pix. makes
>>>> it slower because it creates a geometry first that you do not really need.
>>>>
>>>> Cheers,
>>>>  Tom
>>>>
>>>> ----------
>>>> FULL SCRIPT
>>>>
>>>>
>>>> SELECT
>>>> a.gid As id,
>>>> (ST_SummaryStats(
>>>>     (ST_Union(
>>>>         ST_MapAlgebraExpr(
>>>>             ST_AsRaster(a.geom, b.rast, '32BF')
>>>>             ,b.rast
>>>>             ,'rast2', '32BF','INTERSECTION','0','0',0
>>>>         )
>>>>     )).rast
>>>> )).mean As avg_height
>>>>
>>>> FROM
>>>> polygons.grid a LEFT JOIN
>>>> rasters.elev b
>>>>     ON ST_Intersects(a.geom, b.rast)
>>>> GROUP BY a.gid
>>>>
>>>>
>>>>
>>>> On 24-11-2011 11:14, Andreas Forø Tollefsen wrote:
>>>>
>>>>  Hi,
>>>>
>>>>  I am trying to calculate the average pixel value in a elevation
>>>> raster inside quadrate polygons.
>>>> However, I am not getting the correct values from my query:
>>>>
>>>>  SELECT gid, AVG(((foo.geomval).val)) as avgmnt
>>>> FROM (SELECT p.gid, ST_Intersection(p.cell, r.rast) AS geomval FROM
>>>> mountain r, priogrid_land p WHERE ST_Intersects(p.cell, r.rast, ) AND p.gid
>>>> =186124) AS foo
>>>> GROUP BY gid ORDER BY gid;
>>>>
>>>>  The problem here is that the ST_Intersects(geom, rast) takes into
>>>> consideration the pixels that is outside, but touches the border of the
>>>> quadrate polygons. Then, the average values for each quadrate polygon is
>>>> affected by pixels inside other polygons. This will potentially lead to a
>>>> flawed result.
>>>> So what I want is to be able to calculate the average value for the
>>>> pixels INSIDE the polygon excluding those outside.
>>>>
>>>>  How can i restrict the AVG pixel value to be calculated only for
>>>> pixels that is inside the polygon, and not the pixels that touch the
>>>> outside of the border?
>>>>
>>>>  Thanks!
>>>>
>>>>  Best,
>>>> Andreas
>>>>
>>>>
>>>>
>>>>
>>>>  _______________________________________________
>>>> 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/20111124/6ff5ace4/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: raster_vector_mnt.jpg
Type: image/jpeg
Size: 168534 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20111124/6ff5ace4/attachment.jpg>


More information about the postgis-users mailing list