Joining multiple rasters in a single query

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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)

_______________________________________________
postgis-users mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/postgis-users
Loading...