[postgis-users] Adding bands to a PostGIS raster

Pierre Racine Pierre.Racine at sbf.ulaval.ca
Thu Oct 17 13:10:50 PDT 2013 

If your raster are not too big you can load them not tiled, resample and tile them afterward with ST_Tile().

You can also load everything tiled and add the proper WHERE clause to make sure only matching tiles are added to the base raster table.

The important thing is that the table should be tiled before you intent any further spatial query for performance reason.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-
> bounces at lists.osgeo.org] On Behalf Of Sebastian Schutte
> Sent: Thursday, October 17, 2013 5:41 AM
> To: PostGIS Users Discussion
> Subject: Re: [postgis-users] Adding bands to a PostGIS raster
> 
> Dear Pierre,
> 
> Thanks you so much for pointing me to the problem. If I understand
> correctly, a feasible solution would look like this: I could take all
> the raster data that I have (with different resolutions but the same
> coverage and origin) and load them into the DB without any tiling. After
> that, I could decide on some base resolution (say 1 decimal degree just
> for giggles). I would then resample the data to match that resolution
> and proceed with adding raster bands to my template raster. After that,
> I would break up this template raster into tiles for performance gains.
> 
> Does that sound reasonable?
> 
> Thanks again,
> Sebastian
> 
> On 10/16/2013 04:23 PM, Pierre Racine wrote:
> > Sebastian,
> >
> > Your problem is that you did not manage the tiling properly.
> >
> > 1) You loaded dbe.elevation50km as 3045 tiles (with the -t raster2pgsql
> option) and dbe.elevation20km apparently as 198 tiles (this is odd because
> there should be more tiles in the lower resolution dataset. Are they
> covering the same extent?)
> >
> > Both coverages should have the same alignment (tiles corners should be
> aligned and have the same cell size), have the same number of tiles and
> covering the same extent. You should make sure to get that BEFORe trying
> to merge the bands together.
> >
> > 2) You then want to add bands to the corresponding tiles without
> specifying any condition for this to happen:
> >
> > SELECT     ST_AddBand(toadd.rast, elev.rast, 1) INTO  dbe.raster
> > FROM dbe.test_elev_resample2 AS toadd, dbe.test_elev_resample AS
> elev
> >
> > What this query did was to add all the tiles of the first coverage to all the
> tiles of the second! Ending up with 3045 * 198 = 602910 tiles!
> >
> > When you want to match overlapping (corresponding) tiles you have to
> reduce these combination with a WHERE clause ensuring to match them.
> One trick is to use the upper left x and y. So you would add:
> >
> > WHERE ST_UpperLeftX(elev.rast) = ST_UpperLeftX(toadd.rast) AND
> ST_UpperLeftY(elev.rast) = ST_UpperLeftY(toadd.rast)
> >
> > So I think what you missed was the understand of the tile structures
> inside your tables.
> >
> > Pierre
> >
> >> -----Original Message-----
> >> From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-
> >> bounces at lists.osgeo.org] On Behalf Of Sebastian Schutte
> >> Sent: Wednesday, October 16, 2013 6:36 AM
> >> To: postgis-users at lists.osgeo.org
> >> Subject: [postgis-users] Adding bands to a PostGIS raster
> >>
> >> 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"
> >>
> >> _______________________________________________
> >> postgis-users mailing list
> >> postgis-users at lists.osgeo.org
> >> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> > _______________________________________________
> > postgis-users mailing list
> > postgis-users at lists.osgeo.org
> > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> 
> --
> Sebastian Schutte
> International Conflict Research
> ETH Zurich
> http://www.icr.ethz.ch/people/schutte
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

More information about the postgis-users mailing list