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

[juli g. pausas]

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;
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20150703/daa0ccb6/attachment.html>