[postgis-users] Raster pixel value

[Andreas Forø Tollefsen]

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/624eb25b/attachment.html>