[postgis-users] Joining multiple rasters in a single query

[David Haynes]

Hello,

I have a question regarding the use of joins when clipping multiple raster
datasets in a single query. I have two raster datasets with the same
spatial resolution in separate tables. I perform the clips on them
separately in their own CTEs and then in the last step I would like to do
to use the map algebra to multiply the rasters together. What I am finding
is that the LEFT join is giving the correct answer and the INNER join is
not. I find this strange as both sets should be equal.

Are these assumptions correct?
The clip function should return the same number of tiles, with the same
spatial extents if the same boundary and raster datasets (spatial extent
and resolution) are used? To perform a join on raster tiles one needs to
determine if they have the same spatial extent. Currently I am using the id
fields and the ST_Within(a,b) and ST_Within(b,a) functions. (
http://postgis.net/docs/manual-1.4/ST_Equals.html)

Consider the toy queries below.
With r1 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster1 r on ST_Intersects(r.rast, p.geom)
),
r2 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster2 r on ST_Intersects(r.rast, p.geom)
),
Select r1.id, r1.name, (ST_SummaryStatsAgg(ST_MapAlgebra(r1.rast, 1,
r2.rast, 1, '[rast1]*[rast2]', '32BF'),1, True)).*
FROM r1 left join r2 on r1.id = r2.id and ST_Within(r1.rast, r2.rast)
and ST_Within(r2.rast, r1.rast)


With r1 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster1 r on ST_Intersects(r.rast, p.geom)
),
r2 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster2 r on ST_Intersects(r.rast, p.geom)
),
Select r1.id, r1.name, (ST_SummaryStatsAgg(ST_MapAlgebra(r1.rast, 1,
r2.rast, 1, '[rast1]*[rast2]', '32BF'),1, True)).*
FROM r1 inner join r2 on r1.id = r2.id and ST_Within(r1.rast, r2.rast)
and ST_Within(r2.rast, r1.rast)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20170112/d97c089b/attachment.html>