<html>
  <head>
    <meta content="text/html; charset=ISO-8859-1"
      http-equiv="Content-Type">
  </head>
  <body bgcolor="#FFFFFF" text="#000000">
    Andreas,<br>
    <br>
    Would it be possible to share code and/or sample data?<br>
    With me St_SummaryStats seems to work correct, even when I set all
    my raster values to 1.<br>
    <br>
    Chrs,<br>
    Tom<br>
    <br>
    On 24-11-2011 13:42, Andreas Forø Tollefsen wrote:
    <blockquote
cite="mid:CAGMz7DnCx+++dxC2HK-FGObrHrNH0LbgXex7ZLatoDqOTd1xdQ@mail.gmail.com"
      type="cite">I double checked the values that I get from this query
      and they do not seem to be correct.
      <div>For polygon cells that have only pixels with the value 1
        inside only have a mean of 0.85 (see image).</div>
      <div>The correct value should be a mean of 1. After some visual
        controls I am not sure if these numbers are correct.</div>
      <div>There are many polygon cells with only pixel value 1, but no
        polygon cell with a higher mean than .87.</div>
      <div><br>
      </div>
      <div>Any idea why this could be? Both the raster and polygons are
        in srid 4326.</div>
      <div><br>
      </div>
      <div>Andreas</div>
      <div><br>
      </div>
      <div><br>
        <br>
        <div class="gmail_quote">2011/11/24 Andreas Forø Tollefsen <span
            dir="ltr"><<a moz-do-not-send="true"
              href="mailto:andreasft@gmail.com">andreasft@gmail.com</a>></span><br>
          <blockquote class="gmail_quote" style="margin:0 0 0
            .8ex;border-left:1px #ccc solid;padding-left:1ex;">
            I know. RTFM before asking dumb questions. Set the
            ST_Summary(rast, false) so it does not exclude_nodata_value?
            Right?
            <div><br>
            </div>
            <div><font color="#888888">Andreas</font>
              <div>
                <div class="h5"><br>
                  <br>
                  <div class="gmail_quote">
                    2011/11/24 Andreas Forø Tollefsen <span dir="ltr"><<a
                        moz-do-not-send="true"
                        href="mailto:andreasft@gmail.com"
                        target="_blank">andreasft@gmail.com</a>></span><br>
                    <blockquote class="gmail_quote" style="margin:0 0 0
                      .8ex;border-left:1px #ccc solid;padding-left:1ex">Great.
                      I updated to latest rev and added the st_union.sql
                      to my function list. Now it works.
                      <div>However, I need to ask one more question. </div>
                      <div>When running the script, it calculates the
                        average pixel value inside the polygon. It does
                        not take into account the NULL DATA values.</div>
                      <div>Is there any way i can make the script
                        calculate the average pixel value but count NULL
                        DATA values as 0?</div>
                      <div><br>
                      </div>
                      <div>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?</div>
                      <div><br>
                      </div>
                      <div>Any ideas?</div>
                      <div><br>
                      </div>
                      <div>Thank you so much for this info Tom!</div>
                      <div>
                        <div>
                          <div><br>
                          </div>
                          <div>Andreas<br>
                            <br>
                            <div class="gmail_quote">2011/11/24 Tom van
                              Tilburg <span dir="ltr"><<a
                                  moz-do-not-send="true"
                                  href="mailto:tom.van.tilburg@gmail.com"
                                  target="_blank">tom.van.tilburg@gmail.com</a>></span><br>
                              <blockquote class="gmail_quote"
                                style="margin:0 0 0 .8ex;border-left:1px
                                #ccc solid;padding-left:1ex">
                                <div bgcolor="#FFFFFF" text="#000000">
                                  Yes, I think so.<br>
                                  Current windows build is <b>r8221<br>
                                  </b><br>
                                  Also, I think ST_Union(raster) (which
                                  I used in the example) is not included
                                  in this version yet.<br>
                                  You would have to download a prototype
                                  from. <br>
                                  <a moz-do-not-send="true"
href="http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_union.sql"
                                    target="_blank">http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_union.sql</a><br>
                                  <br>
                                  Chrs,<br>
                                  <font color="#888888">  Tom</font>
                                  <div>
                                    <div><br>
                                      <br>
                                      On 24-11-2011 12:07, Andreas Forø
                                      Tollefsen wrote:
                                      <blockquote type="cite">Great,
                                        Thank you so much for this.
                                        <div>However, I do not seem to
                                          have the
                                          ST_MapAlgebraExpr(rast,
                                          rast.....) function. Was this
                                          implemented in rev 8001?</div>
                                        <div><br>
                                        </div>
                                        <div>Andreas<br>
                                          <br>
                                          <div class="gmail_quote">2011/11/24
                                            Tom van Tilburg <span
                                              dir="ltr"><<a
                                                moz-do-not-send="true"
                                                href="mailto:tom.van.tilburg@gmail.com"
                                                target="_blank">tom.van.tilburg@gmail.com</a>></span><br>
                                            <blockquote
                                              class="gmail_quote"
                                              style="margin:0 0 0
                                              .8ex;border-left:1px #ccc
                                              solid;padding-left:1ex">
                                              <div bgcolor="#FFFFFF"
                                                text="#000000"> Andreas,<br>
                                                <br>
                                                We did approx. the same
                                                thing for non-quadrate
                                                polygons. Perhaps it
                                                might be useful:<br>
                                                <br>
                                                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.<br>
                                                    ST_AsRaster(a.geom,
                                                b.rast,
                                                '<pixtype>')<br>
                                                <br>
                                                Step 2: Overlay the
                                                elevation raster with
                                                the raster you just
                                                created and keep only
                                                the values of the
                                                elevation raster<br>
                                                ST_MapAlgebraExpr(<br>
                                                            <raster
                                                from geometry><br>
                                                            ,b.rast<br>
                                                            ,'rast2' --
                                                <-- keep only raster
                                                2 value<br>
                                                            ,
                                                '<pixtype>','INTERSECTION','0','0',0<br>
                                                        )<br>
                                                <br>
                                                Step 3: get the mean
                                                from the statistics on
                                                the resulting raster<br>
                                                (ST_SummaryStats(<br>
                                                    (ST_Union( -- <
                                                --- we did a UNION
                                                because we occasionaly
                                                had vectors crossing
                                                tiled rasters<br>
                                                        <overlay
                                                raster><br>
                                                    )).rast<br>
                                                )).mean As avg_height  <br>
                                                <br>
                                                That did the trick.
                                                Complete script is
                                                below.<br>
                                                <br>
                                                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.<br>
                                                <br>
                                                Cheers,<br>
                                                 Tom<br>
                                                <br>
                                                ----------<br>
                                                FULL SCRIPT<br>
                                                <br>
                                                <br>
                                                SELECT <br>
                                                a.gid As id,<br>
                                                (ST_SummaryStats(<br>
                                                    (ST_Union(<br>
                                                       
                                                ST_MapAlgebraExpr(<br>
                                                           
                                                ST_AsRaster(a.geom,
                                                b.rast, '32BF')<br>
                                                            ,b.rast<br>
                                                            ,'rast2',
                                                '32BF','INTERSECTION','0','0',0<br>
                                                        )<br>
                                                    )).rast<br>
                                                )).mean As avg_height<br>
                                                <br>
                                                FROM <br>
                                                polygons.grid a LEFT
                                                JOIN <br>
                                                rasters.elev b<br>
                                                    ON
                                                ST_Intersects(a.geom,
                                                b.rast)<br>
                                                GROUP BY a.gid
                                                <div>
                                                  <div><br>
                                                    <br>
                                                    <br>
                                                    On 24-11-2011 11:14,
                                                    Andreas Forø
                                                    Tollefsen wrote: </div>
                                                </div>
                                                <blockquote type="cite">
                                                  <div>
                                                    <div>Hi,
                                                      <div><br>
                                                      </div>
                                                      <div>I am trying
                                                        to calculate the
                                                        average pixel
                                                        value in a
                                                        elevation raster
                                                        inside quadrate
                                                        polygons.</div>
                                                      <div>However, I am
                                                        not getting the
                                                        correct values
                                                        from my query:</div>
                                                      <div><br>
                                                      </div>
                                                      <div>
                                                        <div>SELECT gid,
                                                          AVG(((foo.geomval).val))
                                                          as avgmnt </div>
                                                        <div>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 </div>
                                                        <div>GROUP BY
                                                          gid ORDER BY
                                                          gid;</div>
                                                      </div>
                                                      <div><br>
                                                      </div>
                                                      <div>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.</div>
                                                      <div>So what I
                                                        want is to be
                                                        able to
                                                        calculate the
                                                        average value
                                                        for the pixels
                                                        INSIDE the
                                                        polygon
                                                        excluding those
                                                        outside.</div>
                                                      <div><br>
                                                      </div>
                                                      <div>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?</div>
                                                      <div><br>
                                                      </div>
                                                      <div>Thanks!</div>
                                                      <div><br>
                                                      </div>
                                                      <div>Best,</div>
                                                      <div>Andreas</div>
                                                      <div><br>
                                                      </div>
                                                      <div><br>
                                                      </div>
                                                      <br>
                                                      <fieldset></fieldset>
                                                      <br>
                                                    </div>
                                                  </div>
                                                  <pre>_______________________________________________
postgis-users mailing list
<a moz-do-not-send="true" href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a>
<a moz-do-not-send="true" href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a>
</pre>
                                                </blockquote>
                                                <br>
                                              </div>
                                              <br>
_______________________________________________<br>
                                              postgis-users mailing list<br>
                                              <a moz-do-not-send="true"
href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a><br>
                                              <a moz-do-not-send="true"
href="http://postgis.refractions.net/mailman/listinfo/postgis-users"
                                                target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
                                              <br>
                                            </blockquote>
                                          </div>
                                          <br>
                                        </div>
                                        <br>
                                        <fieldset></fieldset>
                                        <br>
                                        <pre>_______________________________________________
postgis-users mailing list
<a moz-do-not-send="true" href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a>
<a moz-do-not-send="true" href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a>
</pre>
                                      </blockquote>
                                      <br>
                                    </div>
                                  </div>
                                </div>
                                <br>
_______________________________________________<br>
                                postgis-users mailing list<br>
                                <a moz-do-not-send="true"
                                  href="mailto:postgis-users@postgis.refractions.net"
                                  target="_blank">postgis-users@postgis.refractions.net</a><br>
                                <a moz-do-not-send="true"
                                  href="http://postgis.refractions.net/mailman/listinfo/postgis-users"
                                  target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
                                <br>
                              </blockquote>
                            </div>
                            <br>
                          </div>
                        </div>
                      </div>
                    </blockquote>
                  </div>
                  <br>
                </div>
              </div>
            </div>
          </blockquote>
        </div>
        <br>
      </div>
      <br>
      <fieldset class="mimeAttachmentHeader"></fieldset>
      <br>
      <pre wrap="">_______________________________________________
postgis-users mailing list
<a class="moz-txt-link-abbreviated" href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a>
<a class="moz-txt-link-freetext" href="http://postgis.refractions.net/mailman/listinfo/postgis-users">http://postgis.refractions.net/mailman/listinfo/postgis-users</a>
</pre>
    </blockquote>
    <br>
  </body>
</html>