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 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></div><div class="h5"><br><br><div class="gmail_quote">
2011/11/24 Andreas Forø Tollefsen <span dir="ltr"><<a 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><div><div><br></div><div>Andreas<br><br><div class="gmail_quote">2011/11/24 Tom van Tilburg <span dir="ltr"><<a 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 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></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 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 href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a>
<a 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 href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a><br>
            <a 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 href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a>
<a 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 href="mailto:postgis-users@postgis.refractions.net" target="_blank">postgis-users@postgis.refractions.net</a><br>
<a 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>