Thanks Pierre! This is all very helpful! All the best, Sebastian
On 10/17/2013 10:10 PM, Pierre Racine wrote: > 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: [email protected] [mailto:postgis-users- >> [email protected]] 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: [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 > _______________________________________________ > 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
