[postgis-users] Adding bands to a PostGIS raster

[Sebastian Schutte]

Dear PostGIS list,

Thanks for all the great work that many of you have invested into making
PostGIS as awesome as it is. I have a problem contstructing a raster
with multiple bands and could not find a solution in previous threads.
I need to run intersects between many large polygons and raster data
with different resolutions. In order to speed up the process, I decided
to resample the existing raster data to one standard resolution and then
add the data from multiple data sources to different bands of a single
raster. My intuition tells me that this should speed up the intersect
operations as I would simply have to run one spatial query and could
then happily select whichever data I wanted from multiple bands. So much
for the theory.

For test purposed, I uploaded a global elevation dataset into the DB
using different resolutions, one with approximately 50km² cell size near
the Equator and another one with about 20km². SRIDs are the same (4326).

First, I created an empty raster:

INSERT INTO dbe.template(rid,rast)
VALUES(1, ST_MakeEmptyRaster(360, 360, -180.0, 180.0, 2, 2, 0.0, 0.0,
4326));

I then ran ST_Resample to get the rasters into the same resolution as
the template:

SELECT    1 AS rid, ST_Resample(elev.rast, icg.rast,
'NearestNeighbour'::text, 1, false) AS rast
INTO        dbe.test_elev_resample
FROM       dbe.elevation50km AS elev,
                dbe.template AS icg;

/*Query returned successfully: 3045 rows affected, 2307 ms execution time.*/

SELECT     1 AS rid, ST_Resample(elev.rast, icg.rast,
'NearestNeighbour'::text, 1, false) AS rast
INTO         dbe.test_elev_resample2
FROM        dbe.elevation20km AS elev,
                 dbe.template AS icg

/*Query returned successfully: 198 rows affected, 2439 ms execution time.*/

Now I tried to add one band to another raster:

DROP TABLE IF EXISTS dbe.raster;
SELECT     ST_AddBand(toadd.rast, elev.rast, 1)
INTO       dbe.raster
FROM     dbe.test_elev_resample2 AS toadd,
               dbe.test_elev_resample AS elev

/*Query returned successfully: 602910 rows affected, 7449203 ms
execution time.*/

The first thing that's odd is that the operation takes so long (~2hrs)
and that so many rows are affected. When I try to access the bands, the
data is not displayed correctly in pgAdmin and I can select any band
(not just 1 and 2). Of course, adding a band from a raster to the same
raster works fine, but I can also select any band, not just 1 and 2:

DROP TABLE IF EXISTS dbe.raster_tmp;
SELECT     ST_AddBand(raster.rast,raster.rast)
INTO         dbe.raster_tmp
FROM        dbe.test_elev_resample AS raster

/*Query returned successfully: 3045 rows affected, 123 ms execution time.*/

I tried all kinds of variants of these steps and I know I am missing
something very basic here. Its probably my inadequate knowledge of how
raster data is represented in PostGIS. Do the rids have to be set in
order to add bands correctly? I'd be great if somebody could point me to
a solution.

Thanks and all the best,
Sebastian

P.S: The PostGIS version is slightly outdated and its running on Ubuntu
12.04 LTS:
PostGIS_full_version() = "POSTGIS="2.0.1 r9979" GEOS="3.3.3-CAPI-1.7.4"
PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15"
LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER"