[postgis-users] raster, stats conditioned to a set of values
juli g. pausas
juli.g.pausas at uv.es
Mon Jul 6 02:35:57 PDT 2015
Thanks! adding ST_Union helped me for other queries that were producing the
results by tile, but doesn't seems to work for intersection:
INSERT INTO rastertmp.provanet
SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast, 1,
'BAND1'))
FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2
WHERE ST_Intersects(t1.rast, t2.rast);
ERROR: column "t1.rid" must appear in the GROUP BY clause or be used in an
aggregate function
LINE 2: SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast...
Adding a GROUP BY t1.rid, doesn't solve the problem, still get an error.
Both rasters t1 and t2 have tiles.
Juli
--
*CIDE, CSIC* | www.uv.es/jgpausas |
On Fri, Jul 3, 2015 at 8:32 PM, Pierre Racine <Pierre.Racine at sbf.ulaval.ca>
wrote:
> 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;
> >
> >
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20150706/cc1c8259/attachment.html>
More information about the postgis-users
mailing list