[postgis-users] raster, stats conditioned to a set of values

[Pierre Racine]

Yes. You have to:

1) make sure there is a spatial index on your raster table:

CREATE INDEX myrasters_rast_st_convexhull_idx ON myrasters USING gist( ST_ConvexHull(rast) );

2) add ST_Intersects(rast, geom) in the where clause

3) ST_Union() the intersected tiles by polygon id before doing stats

SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast, 1, 'BAND1')) rast
FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2
WHERE ST_Intersects(rast, rast);

Pierre

> -----Original Message-----
> From: juli g. pausas [mailto:juli.g.pausas at uv.es]
> Sent: Friday, July 03, 2015 12:53 PM
> To: Pierre Racine
> Cc: PostGIS Users Discussion
> Subject: Re: [postgis-users] raster, stats conditioned to a set of values
> 
> Well I think I narrowed the problem; the problem is that my map (in this
> example, rastertmp.prova) has many tiles. The original map is very big and I
> cannot loaded to progress without tiles. Apparently the intersection I'm
> doing does not work with a maps with tails. I do not get the error if I use a
> small map with no tiles ... Is it possible to do that intersection with a map
> with tiles?
> 
> Thanks
> 
> 
> Juli
> --
> CIDE, CSIC  |  www.uv.es/jgpausas  |
> 
> 
> 
> On Fri, Jul 3, 2015 at 5:09 PM, juli g. pausas <juli.g.pausas at uv.es> wrote:
> 
> 
> 	Hi
> 
> 	Thanks for this help, it make a lot of sense, but I do not get it to
> work.
> 	I'm new in postgis so to understand this query properly I'm trying
> to run it by steps:
> 
> 
> 	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):
> 
> 
> 	CREATE TABLE rastertmp.provapositiu (rid serial, rast raster);
> 	INSERT INTO rastertmp.provapositiu
> 	    SELECT rid, ST_RECLASS (ST_Band(rast, 1), 1, '[1-10000]:1', '2BUI',
> -9)
> 	    FROM rastertmp.prova;
> 
> 
> 	That works well. To select the positive values, I intersect the original
> raster (prova) with this one (provapositiu):
> 
> 
> 	CREATE TABLE rastertmp.provanet (rid serial, rast raster);
> 	INSERT INTO rastertmp.provanet
> 	    SELECT t1.rid, ST_Intersection(t1.rast, 1, t2.rast, 1, 'BAND1')
> 	    FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2;
> 
> 
> 	And this does not work. Error message:
> 
> 
> 	NOTICE:  The two rasters provided have no intersection.  Returning
> no band raster
> 	CONTEXT:  PL/pgSQL function
> st_intersection(raster,integer,raster,integer,text,double precision[]) line 20
> at assignment
> 
> 
> 	They should intersect, the second comes from the first ...
> 
> 	I thought I could use MapAlgebra; my attempt:
> 
> 	CREATE TABLE rastertmp.provanet (rid serial, rast raster);
> 	INSERT INTO rastertmp.provanet
> 	    SELECT t1.rid, ST_MapAlgebra(t1.rast, 1, t2.rast, 1, '[rast1.val]',
> '16BSI', 'SECOND', '-9', '-9', '-9')
> 	    FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2;
> 
> 
> 	This runs but it  gives a much bigger map than the original ....
> 
> 	Any help?
> 
> 	Thanks
> 
> 
> 
> 
> 	Juli
> 	--
> 	CIDE, CSIC  |  www.uv.es/jgpausas  |
> 
> 
> 
> 	On Thu, Jul 2, 2015 at 4:03 PM, Pierre Racine
> <Pierre.Racine at sbf.ulaval.ca> wrote:
> 
> 
> 		> But my main problem is that I would like to do this (e.g.,
> the query above),
> 		> but only for the pixels in which Band2 = 0.  Any idea? any
> clue?
> 
> 		1) isolate the pixels with 0. For this you can use
> ST_MapAlgebra() or ST_Reclass(). I use ST_MapAlgebra() in the following
> query
> 		2) compute the intersection with band 1
> 		3) compute the stats as you did
> 
> 		SELECT region_cod, (res).*
> 		FROM
> 		  (SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1,
> p.geom, true)) AS res
> 		    FROM gis_wd.wd_regiones AS p, rastertmp.ndvitmp AS r
> 		    WHERE ST_Intersects(r.rast, p.geom)
> 		    AND p.region_cod = 'PA1214'
> 		   ) AS foo WHERE (res).value  > 0;
> 
> 
> 		WITH bands AS ( -- reclass the second band to 0 and 1
> 		SELECT ST_MapAlgebra(ST_Band(rast, 2), '16BSI'::text, 'CASE
> WHEN [rast] < 0 or [rast] > 0 THEN NULL ELSE 1 END') band2,
> 		             ST_Band(rast, 1) band1
> 		) rastintersect AS ( -- compute the intersection of band 1
> and band 2
> 		  SELECT ST_Intersection(band1, band2, 'BAND1') rast FROM
> bands
> 		)
> 		SELECT region_cod, (res).*
> 		FROM
> 		  (SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1,
> p.geom, true)) AS res
> 		    FROM gis_wd.wd_regiones AS p, rastintersect AS r
> 		    WHERE ST_Intersects(r.rast, p.geom)
> 		    AND p.region_cod = 'PA1214'
> 		   ) AS foo;
> 
> 
> 
>