[postgis-users] Adding bands to a PostGIS raster
Sebastian Schutte
sebastian.schutte at icr.gess.ethz.ch
Thu Oct 17 02:40:33 PDT 2013
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
More information about the postgis-users
mailing list