2

Basic Question

I am currently trying to calculate hypsometric areas of the Earth using a GeoTIFF of the ETOPO1 Global Relief map. Specifically, I require a polygon tracing all areas with a bathymetry from 0 to -X00 meters. (I require a polygon (or multipolygon) of this area for later stages in the analysis, therefore solutions for only calculating the area are inadequate).

I have attempted to do this using both QGIS and PostGIS, but I am running into a variety of problems.

My PostGIS workflow

# Import the GeoTIFF into PostGIS using raster2pgsql
raster2pgsql -s 4326 -t 50x50 -f rast -I -Y my.tif etopo_2008_raw | psql MYDB

# Generate contours of the desired area using gdal since there is no 
# ST_contours type function in PostGIS that I am aware of.
gdal_contour -a ELEV -i METERS -f 'ESRI Shapefile' 'PG: dbname=MYDB host=LOCALHOST user=postgres password=PASSWORD port=PORT mode=2 schema=SCHEMA column=rast table=etopo_2008_raw' ~/PATH/etopo_contours

# Load the output shapefile back into postgres
shp2pgsql -I -s 4326 etopo_contours etopo_contours | psql MYDB

Next I use the following SQL query to 1) subset out the relvant contours, do lines to polygons then do the symmetrical difference. I apologize for how hairy this sql is.

CREATE TABLE hypsometric_area_X00 AS 
WITH A AS (SELECT st_makevalid((st_dump(st_polygonize(geom))).geom) AS geom FROM etopo_contours WHERE elev='0'),
B AS (SELECT st_makevalid((st_dump(st_polygonize(geom))).geom) AS geom FROM etopo_contours WHERE elev='-X00')
SELECT ST_SymDifference(C.geom,D.geom) AS geom FROM (SELECT *, row_number() OVER (ORDER BY 1) AS rn FROM A) AS C JOIN (SELECT *,row_number() OVER (ORDER BY 1) AS RN FROM B) AS D ON C.rn=D.rn;

This query does work, but does not give a correct result. This is presumably because the st_polygonize is creating bizarre geometries from the contours. Here is a picture of the output for the 0 to -300 case. Australia is definitely not -300 meters below sea-level.

Incorrect PostGIS Output

Incorrect PostGIS Output

Alternative QGIS Workflow

The QGIS workflow is conceptually identical to the PostGIS workflow.

  • Import the GeoTIFF using add raster layer
  • Use Raster->Conversion->Contour to extract the contours. (Note this uses the exact same gdal query as above.)
  • Subset out the contours that equal = 0 and -X00 as separate line layers.
  • Use Vector->Geometry Tools->Lines to Polygons to create polygons

The resulting polygons look much more sensible than what we got from PostGIS. The 0 Contour As Purple and -300 Contour As Blue

The last step should be relatively simple, use either Vector->Geoprocessing Tools->Difference or Vector->Geoprocessing Tools->Symmetrical Difference to find the 0 to -X00 polygon. However, this attempt consistently crashes or freezes, even if the "ignore incorrect geometries selection is chosen". I've tried it on both a Mac and Windows machine - both running QGIS 2.18.

Summary of Final Question

I am looking for any alternative workflow or modification to the current workflow to make this work in QGIS or PostGIS. (I prefer a PostGIS answer). To clarify the desired final product, I want just the blue part of picture 2 as a separate layer.

P.S. My intern spent a week trying to get this to work in ArcGIS to no avail. I prefer a PostGIS solution because the rest of the project is coded that way, so for consistency and reproducibility it would be better.

4
  • I accidentally deleted my former comment, but it was rubbish anyways. If you suspect your polygonized contours could be the problem, check if your contours are correctly closed with ST_IsClosed. If not, you could e.g. add a new point at the end of each contour line that matches its starting point or use something like ST_ConvexHull to get your polygons from each contour line (I think you can use ST_ConvexHull as a non-aggregate but not sure).
    – geozelot
    Aug 2 '17 at 10:42
  • Great idea to use ST_ConvexHull I will investigate and report back.
    – Andy
    Aug 2 '17 at 14:48
  • I know your requirement is PostGIS solution, but another open source GIS software "SAGA" has module to calculate hypsometric curves. saga-gis.org/saga_tool_doc/2.2.0/ta_morphometry_5.html
    – PKG
    Aug 2 '17 at 16:31
  • Well, I checked and it seems, ConvexHull isn´t the tool you are looking for. But the data, based on your precedure, has indeed, and as expected, several invalid geometries for polygonizing, i.e. countours that are not closed. Update your table and add the starting point to each linestring (via ST_AddPoint) that was identified as not closed (via ST_IsClosed) before polygonizing. Revalidate again that they are corerctly noded. Collect all resulting polygons and subtract (ST_Difference or ST_SymDifference). Worked for me, even my simpler query took ages though.
    – geozelot
    Aug 3 '17 at 9:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.