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
Alternative QGIS Workflow
The QGIS workflow is conceptually identical to the PostGIS workflow.
- Import the GeoTIFF using add raster layer
Raster->Conversion->Contourto 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.
Vector->Geometry Tools->Lines to Polygonsto create polygons
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.