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

[juli g. pausas]

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