[raster] possible performance issue with ST_MapAlgebra

[robe]

In the docs I have this example: ​http://postgis.net/documentation/manual-2.1SVN/RT_ST_MapAlgebraExpr2.html

which is much faster now with r10798 upgrade --on my 64-bit windows 2008 it used to take about 3.6 seconds -- now it's down to 490ms (SUPER improvement)

(
WITH pr AS
-- Note the order of operation: we clip all the rasters to dimensions of our region
(SELECT ST_Clip(rast,ST_Expand(geom,50) ) As rast, g.geom
	FROM aerials.o_2_boston AS r INNER JOIN
-- union our parcels of interest so they form a single geometry we can later intersect with
		(SELECT ST_Union(ST_Transform(the_geom,26986)) AS geom 
		  FROM landparcels WHERE pid IN('0303890000', '0303900000')) As g
		ON ST_Intersects(rast::geometry, ST_Expand(g.geom,50))
),
-- we then union the raster shards together
-- ST_Union on raster is kinda of slow but much faster the smaller you can get the rasters
-- therefore we want to clip first and then union
prunion AS
(SELECT ST_AddBand(NULL, ARRAY[ST_Union(rast,1),ST_Union(rast,2),ST_Union(rast,3)] ) As clipped,geom
FROM pr
GROUP BY geom)
-- return our final raster which is the unioned shard with 
-- with the overlay of our parcel boundaries
-- add first 2 bands, then mapalgebra of 3rd band + geometry
SELECT ST_AddBand(ST_Band(clipped,ARRAY[1,2])
	, ST_MapAlgebraExpr(ST_Band(clipped,3), ST_AsRaster(ST_Buffer(ST_Boundary(geom),2),clipped, '8BUI',250),
	 '[rast2.val]', '8BUI', 'FIRST', '[rast2.val]', '[rast1.val]') ) As rast
FROM prunion;
)

However I thought I could just replace my ST_MapAlgebraExpr with ST_MapAlgebra so I'm not using a deprecated function, and while the output looks the same as the original, my time for this query is: 2320 ms which granted is better than the older timing of ST_MapAlgebraExpr, but not nearly as good as the new timing of ST_MapAlgebraExpr.

(
WITH pr AS
-- Note the order of operation: we clip all the rasters to dimensions of our region
(SELECT ST_Clip(rast,ST_Expand(geom,50) ) As rast, g.geom
	FROM aerials.o_2_boston AS r INNER JOIN
-- union our parcels of interest so they form a single geometry we can later intersect with
		(SELECT ST_Union(ST_Transform(the_geom,26986)) AS geom 
		  FROM landparcels WHERE pid IN('0303890000', '0303900000')) As g
		ON ST_Intersects(rast::geometry, ST_Expand(g.geom,50))
),
-- we then union the raster shards together
-- ST_Union on raster is kinda of slow but much faster the smaller you can get the rasters
-- therefore we want to clip first and then union
prunion AS
(SELECT ST_AddBand(NULL, ARRAY[ST_Union(rast,1),ST_Union(rast,2),ST_Union(rast,3)] ) As clipped,geom
FROM pr
GROUP BY geom)
-- return our final raster which is the unioned shard with 
-- with the overlay of our parcel boundaries
-- add first 2 bands, then mapalgebra of 3rd band + geometry
SELECT ST_AddBand(ST_Band(clipped,ARRAY[1,2])
	, ST_MapAlgebra(ST_Band(clipped,3), ST_AsRaster(ST_Buffer(ST_Boundary(geom),2),clipped, '8BUI',250),
	 '[rast2.val]', '8BUI', 'FIRST', '[rast2.val]', '[rast1.val]') ) As rast
FROM prunion;
)

Is there more to it than just swapping the function names?

[Bborie Park]

There shouldn't be anything more than swapping the function names. What makes it even stranger is that none of the deprecated Map Algebra functions have been worked on in -trunk since 2.0 go branched. Is the source raster used in that example available anywhere? I'd like to replicate that test.

[robe]

I suspect the speed improvement might be the performance improvment changes you made to ST_Clip and ST_Union. I think last I tested ST_Clip had not been converted to pure C.

Not sure which one it is since I have it all in one table chunked. Basically I loaded all of Boston tiles from MassGIS.

with a download transform sid to jpeg like this:

wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22528915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22528930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22678900.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22678915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22678930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22678945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22678975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22678990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22679005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22828900.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22828915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22828930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22828945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22828975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22828990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22829005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978870.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978885.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978900.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22978990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22979005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/22979020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128870.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128885.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128900.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23128975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23129005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23129020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278885.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278900.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23278990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23279005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428900.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23428990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23429005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23429035.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23429050.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23578915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23578930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23578945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23578960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23578975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23578990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23579005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23579020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23579035.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23579050.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23728915.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23728930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23728945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23728960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23728975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23728990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23729005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23729020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23729035.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23878930.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23878960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23878975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23878990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23879005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23879020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/23879035.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24028945.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24028960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24028975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24028990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24029005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24029020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24029035.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24029050.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24178960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24178975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24179005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24179020.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24179035.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24179050.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24328960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24328975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24329005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24478960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24478975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24478990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24479005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24628960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24628975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24628990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24778960.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24778975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24778990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24928975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24928990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/24929005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/25078975.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/25078990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/25079005.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/25228990.zip
wget http://gis.amherstma.gov/massgis/sid_tiles_30cm/25229020.zip

-- unzip command here which uses 7zip on windows --- 
for /f %%a IN ('dir /b bostonaerials2008\*.sid') do "gdal_translate" -of JPEG bostonaerials2008\%%a bostonaerials2008\%%a.jpg


raster2pgsql -I -e -F -C -Y -s 26986 -t 500x500 -l 2 *.jpg aerials.boston | psql

I'll see if I can get a smaller set that exhibits the same behavior. i'm sure you don't want to load all that stuff and the table is kinda big.

parcel layer is on MASSGIS somewhere too though I use one from Boston assessing.

[Bborie Park]

Wow. That'd be a lot of data which is ideal for performance testing. I do have CA aerial imagery and parcel datasets so I'll see about trying it myself.

I'd really like to see what the cold cache performance for that table is in 2.0 and -trunk using single user mode in PostgreSQL. Granted, there is no equivalent to "echo 3 > /proc/sys/vm/drop_caches" in Windows.

Oh. Also the URLs for wget are no longer valid ;-) I see Paul's recent blog entry coming true...

​http://blog.cleverelephant.ca/2012/12/whats-so-hard-about-that-download.html

[robe]

No US is not going the way of Canada. Especially not after I've dissed Canadians for their ability to create complicated systems that go nowhere.

In this case it's just because MassGIS was upgrading servers. Here is their interface for downloading aerials -- even gives you the wget commands you'd need to be efficient.

​http://www.mass.gov/anf/research-and-tech/it-serv-and-support/application-serv/office-of-geographic-information-massgis/datalayers/colororthos2008.html

Once you get to the map interface there is a download spreadsheet that has all the wget commands you'll need for the location you chose.

[Bborie Park]

Thanks! I'm going to test with some Los Angeles parcels on mine using that same query and see what the difference in performance is like. If that doesn't work, I'll use MassGIS.

[robe]

dustymugs - good job. I guess Duncan noticed the improvements too :)

​http://duncanjg.wordpress.com/2012/12/05/ten-fold-speed-up-using-postgis2-1-svn/

[Bborie Park]

I'm suspecting there are some some expression-based ST_MapAlgebra performance issues. I hope to work on this at the Boston Code Sprint...

[Bborie Park]

Fixed in r11222

[robe]

I think this is still an issue.

[robe]

when writing up book playing around with climate data from: ​http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/tmean_26.zip (tiled to 256x256)

IF I do an operation like:

-- takes 7580 ms
INSERT INTO tmean_cel(rid, rast)
SELECT rid,
    ST_MapAlgebra(rast, 1, '32BF'::text, 
     '[rast.val]/10.0'::text, -999)
  FROM tmean AS t LIMIT 10;

-- takes 6480 ms
INSERT INTO tmean_cel(rid, rast)
SELECT rid,
    ST_MapAlgebraExpr(rast, 1, '32BF'::text, 
     '[rast.val]/10.0'::text, -999)
  FROM tmean AS t LIMIT 10;

I'll try to prepare a self-contained example to demonstrate the point

[Bborie Park]

Which minor version of PostGIS are you using?

[robe]

2.1.1

[Bborie Park]

Did you pull all the rasters from the zip file into the tmean table? Can I get the raster2pgsql command?

[robe]

I misspoke I guess I'm running 2.2.0 dev of some variant here.

POSTGIS="2.2.0dev r11899" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24 GDAL_DATA not found" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY (topology procs from "2.2.0dev r11808" need upgrade) RASTER

I pulled all from tmean_26.zip into the table.

raster2pgsql -s 4326 -I -C -M tmean/*.bil -F -t 256x256 tmean | psql

same set of data as #2555

[robe]

yah don't think we'll get to this before 2.1.2 release

[robe]

Milestone renamed