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: [email protected] [mailto:postgis-users- >> [email protected]] On Behalf Of Sebastian Schutte >> Sent: Wednesday, October 16, 2013 6:36 AM >> To: [email protected] >> 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 >> [email protected] >> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users > _______________________________________________ > postgis-users mailing list > [email protected] > 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 [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
