Building contour elevation lines with GDAL and PostGIS

I’m trying to build a complete contour elevation map for the Philippines. I did a couple of them before using various tools. This time, I’m trying to explore how efficient it is to use only CLI tools and cover a much larger area (the whole country). My choice of tools for this process are GDAL and PostGIS.

1. I downloaded DEM tiles from the SRTM CGIAR website (http://srtm.csi.cgiar.org/Index.asp).

2. In order to generate the contour for all the DEM files, I created a VRT to “stitch” all the DEM tiles. This way, I get a continuous contour line along the edges of each DEM tile.
gdalbuildvrt srtm_cgiar_v4.vrt \
1.tif 2.tif 3.tif 4.tif 5.tif \
6.tif 7.tif 8.tif 9.tif 10.tif

3. Using gdal_contour, I created a single contour shapefile at 10 meters interval where the elevation value is in the “elev” column.
gdal_contour -a elev -i 10 srtm_cgiar_v4.vrt srtm_cgiar_v4_10m.shp

4. To add the new shapefile to my existing PosGIS DB, I used a two step process where your first convert the shapefile into a sql file using shp2psql and then with psql, add the the sql file as a new table.
shp2pgsql -s 4326 -I -W "latin1" \
srtm_cgiar_v4_10m.shp > srtm_cgiar_v4_10m.sql
psql -h localhost -d philippines -p 5438 -f srtm_cgiar_v4_10m.sql

5. Finally we run some maintenance process since we inserted a lot of new data in the db.
psql -h localhost -d philippines -p 5438
VACUUM VERBOSE srtm_cgiar_v4_10m;
REINDEX TABLE srtm_cgiar_v4_10m;

Now, I have a complete 10m contour elevation layer I can use in QGIS.
sample_contour

One big problem is that data is so big (1M records) making it very slow to load in QGIS (plus the fact that my test db server is really low in capacity).

In most cases I don’t need the full data, only a subset of what I’m currently working on. To solve this, I’m trying several PostGIS sqls to test the fastest way of loading the layer in QGIS.

Like, clipping a small area using a polygon layer:
SELECT srtm_cgiar_v4_10m.*
INTO srtm_cgiar_v4_10m_ph
FROM srtm_cgiar_v4_10m
INNER JOIN test_clip on
ST_INTERSECTS (adm0.the_geom, srtm_cgiar_v4_10m.the_geom);

Or, using VIEWS:
CREATE OR REPLACE VIEW test_view AS
SELECT srtm_cgiar_v4_10m.*
FROM srtm_cgiar_v4_10m
INNER JOIN test_clip on
ST_INTERSECTS (test_clip.the_geom, srtm_cgiar_v4_10m.the_geom);

Or, in QGIS DB Manager plugin, select only the data within a bounding box:
SELECT * FROM srtm_cgiar_v4_10m
WHERE ST_INTERSECTS
(ST_SetSRID( ST_MakeBox2D
(ST_Point(120.56060,14.27307),
ST_Point(121.70101,15.14615)),
4326 ), the_geom);

So far I liked, the QGIS DB Manager approach since I don’t need to create temporary tables or views. Surely, there are better approaches. Any advice?
I’m fairly new to PostGIS and the sql-way of managing geodata so be patient if my query syntax are not elegant :-)

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s