[postgis-users] Adding bands to a PostGIS raster
Sebastian Schutte
sebastian.schutte at icr.gess.ethz.ch
Fri Oct 18 06:30:50 PDT 2013
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: 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
> _______________________________________________
> 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