postgis_raster - 2

Things change quickly in PostGIS land. Bborie Park has been really hard at work beefing up the PostGIS raster functionality. Just today he committed a change that makes multi-band union raster much easier to write. Recall in last example, Faster Raster Union. I showed this lengthy syntax to get back a mulit-band raster after unioning a multiband raster. This is because in prior versions, the unadorned: ST_Union(rast) was really short-hand for ST_Union(rast,1) or ST_Union(rast,1,'LAST'). This has always driven me nuts. Well today Bborie changed the behavior to be more intuitive. It's a breaking change if you have more than one band in your raster, but honestly I think a behavior people naively assumed anyway. So ST_Union(rast), now means, if you have more than one band in your raster, all your bands will be in the union operation.

So in short, to get: , we can simply write this query:

SELECT ST_Union(rast)
FROM (SELECT  ST_Clip(rast , geom) As rast
	FROM test_100_100 
           CROSS JOIN 
     ST_Buffer(ST_GeomFromText('POINT(231234.436173996 887025.024778253)',26986)
    ,50, 'quad_segs=2') As geom
WHERE ST_Intersects(rast,  geom ) ) As foo

The LAST, FIRST, etc. band number stuff is still useful for many use cases. I'll try to explain that in subsequent articles. They really come into play when we start talking about overlapping raster tiles and then you might also want to combine these with the ORDERED AGG feature supported in PostgreSQL 9.0+ to better control the order of raster tile processing.

Again this feature is available in latest windows experimental buildbot builds. To use just copy over same named folders into your EnterpriseDb Windows install and CREATE EXTENSION postgis VERSION "2.1.0SVN";. Copying over all the files will make 2.1.0SVN the default so CREATE EXTENSION postgis; will do the same. I advice, if you want to work with both 2.0.1 and 2.1.0SVN, you should explicitly specify the version in your CREATE EXTENSION statements so you know what you are getting. Regardless though your geos wil be upgrade to 3.4.0dev (which seems to work fine for me in production) for both 2.0.1 and 2.1.0SVN so you may just want to use on a testing instance.

One of the new enhancements in PostGIS 2.1 and the one I'm most excited about so far is the new faster raster union functionality. Raster Unioning I think is one of the most important functions for raster especially if you have your tiles chunked since for certain cases, you are going to want to glue them back together again.

In PostGIS 2.0.1, the raster ST_Union function piggy-backed a little on plpgsql, in PostGIS 2.1.0, it has been completely reimplemented and reimplemented in C. It's still going thru some optimizations, but as of PostGIS 2.1.0 r10361, you can expect union speed to be about 4-8 times faster. Here is an example. This example still uses ST_Clip which hasn't been optimized yet so you want see the complete speed improvement.

Testing

To test the speed difference, I used a MassGIS 2008 aerial tile I had converted from MrSid to JPEG and loaded with this command:

You can download the test data from 23128870.zip. If you are on Windows and running 9.2 (32-bit/64-bit) or 9.0/9.1 postgresql 64-bit, Winnie the PostGIS windows build bot has fresh binaries you can experiment with you can download from http://postgis.net/windows_downloads

This is I most favorite artsy query to run which produces an output:

SELECT ST_AddBand(NULL,ARRAY[ST_Union(rast,1), ST_Union(rast,2), ST_Union(rast,3) ]) 
FROM (SELECT  ST_Clip(rast,geom) As rast 
	FROM test_100_100 
           CROSS JOIN ST_Buffer(ST_GeomFromText('POINT(231237.436173996 887025.024778253)',26986),50, 'quad_segs=2') As geom
WHERE ST_Intersects(rast,  geom ) ) As foo;

and takes about 4.5 seconds on PostGIS 2.0.1, PostgreSQL 9.2.1 windows 7 64-bit with default factory settings and has to run thru 16 tiles each tile being 100x100.

If you run this query in PostGIS 2.1.0, it goes down to about 1,560ms.

SELECT ST_Union(rast,ARRAY[ROW(1, 'LAST'), ROW(2, 'LAST'), ROW(3, 'LAST')]::unionarg[])
FROM (SELECT  ST_Clip(rast , geom) As rast
	FROM test_100_100 
           CROSS JOIN ST_Buffer(ST_GeomFromText('POINT(231234.436173996 887025.024778253)',26986),50, 'quad_segs=2') As geom
WHERE ST_Intersects(rast,  geom ) ) As foo

Although the 2.0.1 syntax is still supported, the second query (preferred 2.1 syntax) runs about 300ms faster for this particular query.

The speed difference is more startling if we compare unclipped speeds. For unclipped speeds, the 2.1.0 ST_Union still finishes at about 1,560ms, while the PostGIS 2.0.1 goes up to 6,712ms.