To preface, I'm not sure if it matters, but I'm using PostgreSQL 11 with PostGIS 2.5 (through Docker on Windows) and accessing it with Python through PsycoPG2.
My raster has the recommended index on image USING GIST (ST_ConvexHull(image)).
So far, it seems like I should have stuck with the non-PostGIS solution that using OpenCV for speed. At least until I inserted a ST_Rescale(image, 1.0) into the query (it would be replaced with a logical value for individual queries, I just had it as a placeholder).
Now I have:
for line in mos.conn.fetchall("""EXPLAIN ANALYZE
WITH foo as (SELECT
ST_MakeEmptyRaster(%(width)s, %(height)s,
%(xoffset)s, %(yoffset)s,
%(xscale)s, %(yscale)s,
%(xskew)s, %(yskew)s,
ST_SRID(image)
) AS rast
FROM image INNER JOIN mosaic_part ON image.image_pkey=mosaic_part.image_id
WHERE mosaic_id = %(mos_id)s AND ST_SRID(image) != 0
LIMIT 1)
SELECT ST_Clip(ST_Rescale(image, 1.0), ST_Envelope(rast))
FROM image INNER JOIN mosaic_part ON image.image_pkey=mosaic_part.image_id
CROSS JOIN foo
WHERE mosaic_id = %(mos_id)s AND ST_SRID(image) = ST_SRID(rast) AND ST_Intersects(image, rast)
""", {"mos_id": mos.id, "width": size[1], "height": size[0],
"xoffset": matrix[0,2], "yoffset": matrix[1,2],
"xscale": matrix[0,0], "yscale": matrix[1,1],
"xskew": matrix[0,1], "yskew": matrix[1,0]}):
print(line[0])
yielding:
Nested Loop (cost=2.02..12.43 rows=1 width=32) (actual time=40.618..549.907 rows=21 loops=1)
CTE foo
-> Limit (cost=0.56..0.84 rows=1 width=32) (actual time=3.616..3.621 rows=1 loops=1)
-> Merge Join (cost=0.56..203.59 rows=703 width=32) (actual time=3.614..3.615 rows=1 loops=1)
Merge Cond: (image_1.image_pkey = mosaic_part_1.image_id)
-> Index Scan using image_pkey on image image_1 (cost=0.28..98.92 rows=1176 width=26) (actual time=1.877..1.878 rows=1 loops=1)
Filter: (st_srid(image) <> 0)
-> Index Only Scan using mosaic_part_pkey on mosaic_part mosaic_part_1 (cost=0.28..89.42 rows=707 width=4) (actual time=0.026..0.026 rows=1 loops=1)
Index Cond: (mosaic_id = '[...]'::text)
Heap Fetches: 1
-> Nested Loop (cost=0.90..10.70 rows=1 width=58) (actual time=22.573..378.841 rows=66 loops=1)
-> CTE Scan on foo (cost=0.00..0.02 rows=1 width=32) (actual time=3.621..3.627 rows=1 loops=1)
-> Index Scan using idx_image_bounds on image (cost=0.90..10.67 rows=1 width=26) (actual time=18.932..375.092 rows=66 loops=1)
Index Cond: ((image)::geometry && (foo.rast)::geometry)
Filter: ((st_srid(foo.rast) = st_srid(image)) AND _st_intersects(st_convexhull(image), st_convexhull(foo.rast)))
Rows Removed by Filter: 29
-> Index Only Scan using mosaic_part_pkey on mosaic_part (cost=0.28..0.38 rows=1 width=4) (actual time=0.006..0.006 rows=0 loops=66)
Index Cond: ((mosaic_id = '[...]'::text) AND (image_id = image.image_pkey))
Heap Fetches: 21
Planning time: 0.538 ms
Execution time: 550.020 ms
Where it was previously:
for line in mos.conn.fetchall("""EXPLAIN ANALYZE
WITH foo as (SELECT
ST_MakeEmptyRaster(%(width)s, %(height)s,
%(xoffset)s, %(yoffset)s,
%(xscale)s, %(yscale)s,
%(xskew)s, %(yskew)s,
ST_SRID(image)
) AS rast
FROM image INNER JOIN mosaic_part ON image.image_pkey=mosaic_part.image_id
WHERE mosaic_id = %(mos_id)s AND ST_SRID(image) != 0
LIMIT 1)
SELECT ST_Clip(image, ST_Envelope(rast))
FROM image INNER JOIN mosaic_part ON image.image_pkey=mosaic_part.image_id
CROSS JOIN foo
WHERE mosaic_id = %(mos_id)s AND ST_SRID(image) = ST_SRID(rast) AND ST_Intersects(image, rast)
""", {"mos_id": mos.id, "width": size[1], "height": size[0],
"xoffset": matrix[0,2], "yoffset": matrix[1,2],
"xscale": matrix[0,0], "yscale": matrix[1,1],
"xskew": matrix[0,1], "yskew": matrix[1,0]}):
print(line[0])
yielding:
Nested Loop (cost=2.02..12.18 rows=1 width=32) (actual time=32.914..4972.957 rows=21 loops=1)
CTE foo
-> Limit (cost=0.56..0.84 rows=1 width=32) (actual time=3.441..3.446 rows=1 loops=1)
-> Merge Join (cost=0.56..203.59 rows=703 width=32) (actual time=3.439..3.439 rows=1 loops=1)
Merge Cond: (image_1.image_pkey = mosaic_part_1.image_id)
-> Index Scan using image_pkey on image image_1 (cost=0.28..98.92 rows=1176 width=26) (actual time=1.747..1.748 rows=1 loops=1)
Filter: (st_srid(image) <> 0)
-> Index Only Scan using mosaic_part_pkey on mosaic_part mosaic_part_1 (cost=0.28..89.42 rows=707 width=4) (actual time=0.020..0.021 rows=1 loops=1)
Index Cond: (mosaic_id = '[...]'::text)
Heap Fetches: 1
-> Nested Loop (cost=0.90..10.70 rows=1 width=58) (actual time=23.453..363.824 rows=66 loops=1)
-> CTE Scan on foo (cost=0.00..0.02 rows=1 width=32) (actual time=3.446..3.458 rows=1 loops=1)
-> Index Scan using idx_image_bounds on image (cost=0.90..10.67 rows=1 width=26) (actual time=19.988..360.250 rows=66 loops=1)
Index Cond: ((image)::geometry && (foo.rast)::geometry)
Filter: ((st_srid(foo.rast) = st_srid(image)) AND _st_intersects(st_convexhull(image), st_convexhull(foo.rast)))
Rows Removed by Filter: 29
-> Index Only Scan using mosaic_part_pkey on mosaic_part (cost=0.28..0.38 rows=1 width=4) (actual time=0.006..0.006 rows=0 loops=66)
Index Cond: ((mosaic_id = '[...]'::text) AND (image_id = image.image_pkey))
Heap Fetches: 21
Planning time: 0.493 ms
Execution time: 4973.083 ms
The eventual plan is to resample and union the results to the size of the empty raster foo.rast. The resample task seems to be extremely slow in PostGIS compared to getting the images and doing it myself, which will be a later question. If there's a better way to do this in general, I'm all ears. I'm new to PostGIS.
