Andreas,

Would it be possible to share code and/or sample data?
With me St_SummaryStats seems to work correct, even when I set all my raster values to 1.

Chrs,
Tom

On 24-11-2011 13:42, Andreas Forø Tollefsen wrote:
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 <[email protected] <mailto:[email protected]>>

    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 <[email protected]
    <mailto:[email protected]>>

        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 <[email protected]
        <mailto:[email protected]>>

            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 <[email protected]
            <mailto:[email protected]>>

                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 list
                [email protected]  
<mailto:[email protected]>
                http://postgis.refractions.net/mailman/listinfo/postgis-users


                _______________________________________________
                postgis-users mailing list
                [email protected]
                <mailto:[email protected]>
                http://postgis.refractions.net/mailman/listinfo/postgis-users




            _______________________________________________
            postgis-users mailing list
            [email protected]  
<mailto:[email protected]>
            http://postgis.refractions.net/mailman/listinfo/postgis-users


            _______________________________________________
            postgis-users mailing list
            [email protected]
            <mailto:[email protected]>
            http://postgis.refractions.net/mailman/listinfo/postgis-users






_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users

_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users

Reply via email to