<div dir="ltr"><div><div>Hi<br></div>Thanks for this help, it make a lot of sense, but I do not get it to work.<br>I'm new in postgis so to understand this query properly I'm trying to run it by steps:<br><br></div>I have a raster, rastertmp.prova, that have positive and negative values, and I'd like to use the positive values only. As suggested, I generate a raster that is 1 (for positive values) or Null (otherwise):<br><div><br>CREATE TABLE rastertmp.provapositiu (rid serial, rast raster);<br>INSERT INTO rastertmp.provapositiu<br> SELECT rid, ST_RECLASS (ST_Band(rast, 1), 1, '[1-10000]:1', '2BUI', -9) <br> FROM rastertmp.prova;<br><br></div><div>That works well. To select the positive values, I intersect the original raster (prova) with this one (provapositiu):<br></div><div><br>CREATE TABLE rastertmp.provanet (rid serial, rast raster);<br>INSERT INTO rastertmp.provanet<br> SELECT t1.rid, ST_Intersection(t1.rast, 1, t2.rast, 1, 'BAND1')<br> FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2;<br><br></div><div>And this does not work. Error message: <br></div><div><br>NOTICE: The two rasters provided have no intersection. Returning no band raster<br>CONTEXT: PL/pgSQL function st_intersection(raster,integer,raster,integer,text,double precision[]) line 20 at assignment<br><br></div><div>They should intersect, the second comes from the first ...<br></div><div>I thought I could use MapAlgebra; my attempt:<br><br>CREATE TABLE rastertmp.provanet (rid serial, rast raster);<br>INSERT INTO rastertmp.provanet<br> SELECT t1.rid, ST_MapAlgebra(t1.rast, 1, t2.rast, 1, '[rast1.val]', '16BSI', 'SECOND', '-9', '-9', '-9') <br> FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2;<br><br></div><div>This runs but it gives a much bigger map than the original ....<br></div><div>Any help?<br></div><div>Thanks<br><br></div></div><div class="gmail_extra"><br clear="all"><div><div class="gmail_signature"><div dir="ltr"><div><div dir="ltr">Juli<div><span style="color:rgb(51,51,153)">--</span><br><font color="#333399"><b>CIDE, CSIC</b> | <a href="http://www.uv.es/jgpausas" target="_blank">www.uv.es/jgpausas</a> | <br><br></font></div></div></div></div></div></div>
<br><div class="gmail_quote">On Thu, Jul 2, 2015 at 4:03 PM, Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca" target="_blank">Pierre.Racine@sbf.ulaval.ca</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">> But my main problem is that I would like to do this (e.g., the query above),<br>
> but only for the pixels in which Band2 = 0. Any idea? any clue?<br>
<br>
</span>1) isolate the pixels with 0. For this you can use ST_MapAlgebra() or ST_Reclass(). I use ST_MapAlgebra() in the following query<br>
2) compute the intersection with band 1<br>
3) compute the stats as you did<br>
<span class=""><br>
SELECT region_cod, (res).*<br>
FROM<br>
(SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1, p.geom, true)) AS res<br>
FROM gis_wd.wd_regiones AS p, rastertmp.ndvitmp AS r<br>
WHERE ST_Intersects(r.rast, p.geom)<br>
AND p.region_cod = 'PA1214'<br>
) AS foo WHERE (res).value > 0;<br>
<br>
<br>
</span>WITH bands AS ( -- reclass the second band to 0 and 1<br>
SELECT ST_MapAlgebra(ST_Band(rast, 2), '16BSI'::text, 'CASE WHEN [rast] < 0 or [rast] > 0 THEN NULL ELSE 1 END') band2,<br>
ST_Band(rast, 1) band1<br>
) rastintersect AS ( -- compute the intersection of band 1 and band 2<br>
SELECT ST_Intersection(band1, band2, 'BAND1') rast FROM bands<br>
<span class="">)<br>
SELECT region_cod, (res).*<br>
FROM<br>
(SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1, p.geom, true)) AS res<br>
</span> FROM gis_wd.wd_regiones AS p, rastintersect AS r<br>
<span class=""> WHERE ST_Intersects(r.rast, p.geom)<br>
AND p.region_cod = 'PA1214'<br>
</span> ) AS foo;<br>
<br>
</blockquote></div><br></div>