Calculating Slope for the Whole World in PostGIS

Question

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)
SELECT ST_Clip(ST_Slope(ST_Union(b.rast), 1, '32BF', 'DEGREES', 111120), ST_Envelope(a.rast)) s
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
    ST_SRID(a.rast)
  )
)
GROUP BY a.rast;

The Cursor Approach

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

CREATE OR REPLACE FUNCTION slope_cx()
RETURNS void AS
$BODY$
DECLARE
    data CURSOR FOR select * from gdem_100;
    data_row record;
BEGIN
    OPEN data;

    LOOP
    FETCH data INTO data_row;
        EXIT WHEN NOT FOUND;

        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
            ST_SRID(data_row.rast)
          )
        );

    END LOOP;
    CLOSE data;
END;
$BODY$
LANGUAGE plpgsql VOLATILE;

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?

Answer 1

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 (
    SELECT 
     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
      )
     ,ST_Envelope(a.rast)
    ) 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.

Answer 2

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() 
RETURNS TABLE (
rast   raster) AS
$BODY$
DECLARE
    r elevation.srtm%rowtype;
BEGIN
    FOR r IN SELECT * FROM elevation.srtm as srtm
    LOOP
        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;
    END LOOP;
    RETURN;
END
$BODY$
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!