The Challenge

I have a PostGIS Database which contains all of SRTM v4.1. I need to calculate slope for the entire globe as part of an error modeling process. The DEM is tiled at 100x100.

Whenever I try to calculate slope on the tiled data, I get edge effects. I've written about it and how to avoid it for small areas, but I can't just ST_Union the globe and expect it all to be fine. Even for just the state of Colorado I get a memory allocation error.

The Plan

  1. For each tile in the table, calculate the bounding box.
  2. Expand that bounding box by 1 pixel on each side.
  3. Using ST_Intersects, select all tiles that intersect the bounding box.
  4. Union that 9-tile set, and calculate the slope.
  5. Clip the resulting much larger raster with the original 100x100 bounding box
  6. Insert my 100x100 slope tile into a new table.

Attempted Implementations

First Try

Based on tilt's answer, I've tried this query:

INSERT INTO tmp_slope(rast)
SELECT ST_Clip(ST_Slope(ST_Union(b.rast), 1, '32BF', 'DEGREES', 111120), ST_Envelope(a.rast))
FROM dem_foco_100 a
INNER JOIN dem_foco_100 b
ON (ST_Intersects(a.rast, b.rast))
GROUP BY a.rast;

Unfortunately, as you can see from the links in my comments on his answer, this is still producing edge effects for me. It seems that this is the case because ST_Intersects doesn't return all neighboring tiles as expected.

Expanding the Bounding Box

Instead, I've created a query based on expanding the bounding box of the current tile with the JOIN route proposed by tilt to come up with something that works for my subset.

INSERT INTO tmp_slope(rast)
FROM (SELECT * FROM dem_foco_100) a
INNER JOIN dem_foco_100 b
ON ST_Intersects(
  ST_Envelope(b.rast),  -- All raster tiles
  ST_MakeEnvelope(      -- Expanded bouding box
    -- WEST
    (ST_UpperLeftX(a.rast) - ST_PixelWidth(a.rast)),
    -- SOUTH
    (ST_UpperLeftY(a.rast) - ((ST_Height(a.rast) + 1) * ST_PixelHeight(a.rast))),
    -- EAST
    (ST_UpperLeftX(a.rast) + ((ST_Width(a.rast) + 1) * ST_PixelWidth(a.rast))),
    -- NORTH
    (ST_UpperLeftY(a.rast) + ST_PixelHeight(a.rast)),
    -- Use the OG SRID
GROUP BY a.rast;

The Cursor Approach

I've also tried a cursor using the same bounding box logic:

    data CURSOR FOR select * from gdem_100;
    data_row record;
    OPEN data;

    FETCH data INTO data_row;

        INSERT INTO cur_slope(rast)
        SELECT ST_Clip(ST_Slope(ST_Union(a.rast), 1, '32BF', 'DEGREES', 111120), ST_Envelope(data_row.rast))
        FROM gdem_100 a
        WHERE ST_Intersects(
          ST_Envelope(a.rast),  -- All raster tiles
          ST_MakeEnvelope(      -- Expanded bouding box
            -- WEST
            (ST_UpperLeftX(data_row.rast) - ST_PixelWidth(data_row.rast)),
            -- SOUTH
            (ST_UpperLeftY(data_row.rast) - ((ST_Height(data_row.rast) + 1) * ST_PixelHeight(data_row.rast))),
            -- EAST
            (ST_UpperLeftX(data_row.rast) + ((ST_Width(data_row.rast) + 1) * ST_PixelWidth(data_row.rast))),
            -- NORTH
            (ST_UpperLeftY(data_row.rast) + ST_PixelHeight(data_row.rast)),
            -- Use the OG SRID

    CLOSE data;

Both of these queries do what I want for small areas (< 600 tiles). Unfortunately, both fail on my full test area (~256000 tiles). Neither actually errored out, but I cancelled them after about 40 hours.

I need to cover the full area in an hour or less of run time for this to be feasible to run on the full GDEM table, as the land area of my full test dataset is ~1/90th the area of North America alone.

How can I write a performant query or function to calculate edge-less slope for the entire world?


There is not really a need to put this into plpgsql. I guess it even causes the problem of loading everything into memory first. Here's the query that works for me (on a considerably smaller dataset). As a bonus I've added funcionality that exports your data to a tiff file, though I can imagine this may also

--Putting in next to lines because I always forget about them myself. Newer versions of postgres will complain without it.
SET postgis.enable_outdb_rasters = True;
SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';

--Creating a brand new table but you might as well make a insert or update of it
CREATE TABLE tmp.tmp_out AS 
WITH sloped AS (
     ST_Clip( --only keep the tile that we started with by cutting it out based on ST_Envelope of the original raster (a.rast)
       ST_Slope( --calculate slope on set of tiles
         ST_Union(b.rast) --combine all neighbour tiles and our own tile
    ) AS slope
    FROM gdem_100 a 
    --The trick here is to do a join on the same table, and then see what tiles in your table are neigbours
    INNER JOIN gdem_100  b ON ST_Intersects(a.rast, b.rast)
    GROUP BY a.rast -- Don't forget to group by, or psql will complain
--This lo_from_bytea is only for exporting, I've commented out now because it may cause confusion and is not part of the answer to your question.
/*SELECT lo_from_bytea(0, ST_AsTiff(ST_Union(slope))) as loid*/
SELECT slope AS rast
FROM sloped;
--You still have to find a way to get your raster out of the database...    

--same here, only for exporting
SELECT lo_export(loid, '/tmp/slope.tiff')
FROM tmp.tmp_out;

SELECT lo_unlink(loid)
FROM tmp.tmp_out;

You will notice that the query is relatively slow (much slower than what you would get with using GDAL) but stable.

I assume by the way that you've build an index on your raster table.

edit: The advantage of using a database for a large amount of data is that it will take care of planning the calculations (e.g. doing a union) and looking up data that is (spatially) related (St_Intersects). Therefore you don't have to load all your data into memory at once, the query planner will make sure it happens one row at a time. The disadvantage is that the process will likely be slower than calculating all of this in memory.

I have also added some comments within the query, hope it helps.

  • Yeah, table is indexed using GIST on the convex hull of the tiles. Can you explain a little bit about what this query is doing and why it's effective when mine isn't? I'm pretty new to SQL and a solid conceptualization of joins still evades me. – nronnei May 5 '17 at 14:56
  • My understanding of joins is that they create an "intermediate table" (just don't know what to call it, a view?) that includes all of the rows that exist in the other table, in our case, if there's a tile in each location. Since it's the same table, why do we do that? Why do we duplicate the data? – nronnei May 5 '17 at 15:06
  • @nronnei: JOINS are the most important part of relational databases. It is important that you understand them so please check some basic tutorials on internet. Spatial joins like this can be difficult to start with (I still make mistakes with them). As I said: the planner will make sure for you that data isn't unnecessarily duplicated. There are things like intermediate tables but that doesn't happen with joins like this. – tilt May 5 '17 at 15:21
  • I've been unable to confirm that this works because my connection breaks before the process completes for my test area. The DB is running on 1 core and is localhost, not a remote server. SSL SYSCALL error: EOF detected The connection to the server was lost. Attempting reset: Failed. According to what I've seen, this could be caused by a number of things including memory and disk issues. – nronnei May 6 '17 at 21:15
  • I've tried your query on much smaller subset, this doesn't work. I still get edge effects. Here are the results of my query, and the one you've offered. It appears that the edge effects have just scaled out. – nronnei May 8 '17 at 0:34

I faced the same issue with a big coverage, I solved it using plpsql. I found out necessary to use a buffer (maybe not the best approach, but it worked) in order to make sure all the surrounding tiles were included.

--DROP FUNCTION elevation.create_slope();
CREATE OR REPLACE FUNCTION elevation.create_slope() 
rast   raster) AS
    r elevation.srtm%rowtype;
    FOR r IN SELECT * FROM elevation.srtm as srtm
        RETURN query select ST_Slope(ST_Union(rr.rast), 1, r.rast,'32BF','PERCENT',111120) as rast
            from elevation.srtm as rr
            -- I used a small distance within otherwise some tiles may not be included
            where ST_DWithin(r.rast,rr.rast,0.001)
            GROUP BY r.rid;
LANGUAGE 'plpgsql' ;

create table elevation.slope as
        SELECT * FROM elevation.create_slope();

This took around 1h:30m to compute a coverage with around 1 500 000 tiles of 15x15 px.

In the end don't forget to add serial primary key; GIST; constrains etc to your new table!

Your Answer

