<div dir="ltr"><div><div>Thanks! adding ST_Union helped me for other queries that were producing the results by tile, but doesn't seems to work for intersection:<br><br>INSERT INTO rastertmp.provanet<br> SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast, 1, 'BAND1'))<br> FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2<br> WHERE ST_Intersects(t1.rast, t2.rast);<br><br>ERROR: column "t1.rid" must appear in the GROUP BY clause or be used in an aggregate function<br>LINE 2: SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast...<br><br></div>Adding a GROUP BY t1.rid, doesn't solve the problem, still get an error.<br></div>Both rasters t1 and t2 have tiles.<br><div><div><br><div><br></div></div></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 Fri, Jul 3, 2015 at 8:32 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">Yes. You have to:<br>
<br>
1) make sure there is a spatial index on your raster table:<br>
<br>
CREATE INDEX myrasters_rast_st_convexhull_idx ON myrasters USING gist( ST_ConvexHull(rast) );<br>
<br>
2) add ST_Intersects(rast, geom) in the where clause<br>
<br>
3) ST_Union() the intersected tiles by polygon id before doing stats<br>
<br>
SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast, 1, 'BAND1')) rast<br>
<span class="">FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2<br>
</span>WHERE ST_Intersects(rast, rast);<br>
<br>
Pierre<br>
<div class="HOEnZb"><div class="h5"><br>
> -----Original Message-----<br>
> From: juli g. pausas [mailto:<a href="mailto:juli.g.pausas@uv.es">juli.g.pausas@uv.es</a>]<br>
> Sent: Friday, July 03, 2015 12:53 PM<br>
> To: Pierre Racine<br>
> Cc: PostGIS Users Discussion<br>
> Subject: Re: [postgis-users] raster, stats conditioned to a set of values<br>
><br>
> Well I think I narrowed the problem; the problem is that my map (in this<br>
> example, rastertmp.prova) has many tiles. The original map is very big and I<br>
> cannot loaded to progress without tiles. Apparently the intersection I'm<br>
> doing does not work with a maps with tails. I do not get the error if I use a<br>
> small map with no tiles ... Is it possible to do that intersection with a map<br>
> with tiles?<br>
><br>
> Thanks<br>
><br>
><br>
> Juli<br>
> --<br>
> CIDE, CSIC | <a href="http://www.uv.es/jgpausas" rel="noreferrer" target="_blank">www.uv.es/jgpausas</a> |<br>
><br>
><br>
><br>
> On Fri, Jul 3, 2015 at 5:09 PM, juli g. pausas <<a href="mailto:juli.g.pausas@uv.es">juli.g.pausas@uv.es</a>> wrote:<br>
><br>
><br>
> Hi<br>
><br>
> Thanks for this help, it make a lot of sense, but I do not get it to<br>
> work.<br>
> I'm new in postgis so to understand this query properly I'm trying<br>
> to run it by steps:<br>
><br>
><br>
> I have a raster, rastertmp.prova, that have positive and negative<br>
> values, and I'd like to use the positive values only. As suggested, I generate<br>
> a raster that is 1 (for positive values) or Null (otherwise):<br>
><br>
><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',<br>
> -9)<br>
> FROM rastertmp.prova;<br>
><br>
><br>
> That works well. To select the positive values, I intersect the original<br>
> raster (prova) with this one (provapositiu):<br>
><br>
><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>
><br>
> And this does not work. Error message:<br>
><br>
><br>
> NOTICE: The two rasters provided have no intersection. Returning<br>
> no band raster<br>
> CONTEXT: PL/pgSQL function<br>
> st_intersection(raster,integer,raster,integer,text,double precision[]) line 20<br>
> at assignment<br>
><br>
><br>
> They should intersect, the second comes from the first ...<br>
><br>
> 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]',<br>
> '16BSI', 'SECOND', '-9', '-9', '-9')<br>
> FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2;<br>
><br>
><br>
> This runs but it gives a much bigger map than the original ....<br>
><br>
> Any help?<br>
><br>
> Thanks<br>
><br>
><br>
><br>
><br>
> Juli<br>
> --<br>
> CIDE, CSIC | <a href="http://www.uv.es/jgpausas" rel="noreferrer" target="_blank">www.uv.es/jgpausas</a> |<br>
><br>
><br>
><br>
> On Thu, Jul 2, 2015 at 4:03 PM, Pierre Racine<br>
> <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>> wrote:<br>
><br>
><br>
> > But my main problem is that I would like to do this (e.g.,<br>
> the query above),<br>
> > but only for the pixels in which Band2 = 0. Any idea? any<br>
> clue?<br>
><br>
> 1) isolate the pixels with 0. For this you can use<br>
> ST_MapAlgebra() or ST_Reclass(). I use ST_MapAlgebra() in the following<br>
> query<br>
> 2) compute the intersection with band 1<br>
> 3) compute the stats as you did<br>
><br>
> SELECT region_cod, (res).*<br>
> FROM<br>
> (SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1,<br>
> 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>
> WITH bands AS ( -- reclass the second band to 0 and 1<br>
> SELECT ST_MapAlgebra(ST_Band(rast, 2), '16BSI'::text, 'CASE<br>
> 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<br>
> and band 2<br>
> SELECT ST_Intersection(band1, band2, 'BAND1') rast FROM<br>
> bands<br>
> )<br>
> SELECT region_cod, (res).*<br>
> FROM<br>
> (SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1,<br>
> p.geom, true)) AS res<br>
> FROM gis_wd.wd_regiones AS p, rastintersect AS r<br>
> WHERE ST_Intersects(r.rast, p.geom)<br>
> AND p.region_cod = 'PA1214'<br>
> ) AS foo;<br>
><br>
><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>