Avoiding collisions between contour lines and waterbody polygons

Following up on my earlier posts on contour lines, I have devised a method and a script which takes care of another issue that comes with generating contours from a DEM. The script can be found at the bottom of this page.

The problem can be seen above. Contours tend to get close to or even intersect with waterbodies. This is not unexpected as the lakes may be generalized in their shape or not conform exactly to the DEM for other reasons.

These contours are not very effective in their task and contribute to an overall disorganized appearance of the map. Cleaning up the contours by hand is time-consuming. An automated way of fixing this would be desirable.

Result

Above are the resulting contours. There is a consistent minimum distance to all waterbodies and no collisions.

Theory of operation

This is how the script works under the hood:

1. For each vector polygon, sample the DEM and save the median value to a column.
2. Round the median values using a step function so that the resulting values will lay perfectly between the values where contour lines will be drawn. When using a contour interval of 10 meters, the lake polygons will have elevation values of 5, 15, 25 and so on. This will ensure a consistent distance from lake to contour for all lakes.
3. Buffer the polygons, distance depending on the target map scale. This will be the allowed minimum distance in meters between lakes and contours.
4. Rasterize the buffered polygons burning in the rounded values.
5. Finally, merge the original DEM with the newly rasterized file, laying the lakes on top.

The output file can then be used to generate contours.

NB: There will be some jagged edges around lakes, I recommend adding a 1 px Gaussian blur to the output DEM to smooth these out.

Shell script

lakeFix.sh

#!/bin/sh

if [ -z "$1" ]; then
    echo "Missing arguments. Syntax:"
    echo "   lakeFix.sh dem.tif lakes.gpkg scale[50/25] desired_contour_interval"
    exit 2
fi

set -o errexit
set -o nounset
set -o pipefail

dem="$1"
lakes="$2"
scale="$3"
desired_interval="$4"
halfInterval=$(echo $4/2 | bc -l)

#Get meters per pixel of DEM
pixelsize=$(gdalinfo $dem | grep Pixel | sed 's/Pixel Size = //g' | sed 's/(//g' | sed 's/)//g' | awk -F ',' '{ print $1 }')

#Get extent of DEM
gdaltindex extent.geojson $dem
EXTENT=$(ogrinfo -al -so extent.geojson | grep Extent | sed 's/Extent: //g' | sed 's/(//g' | sed 's/)//g' | sed 's/ - /, /g')
rm extent.geojson

# Hack to make qgis_process happy
export PROJ_LIB=/Applications/QGIS.app/Contents/Resources/proj/

# For each lake polygon, get the median value from the DEM, save to "_median"-column
/Applications/QGIS.app/Contents/MacOS/qgis_process.app/Contents/MacOS/qgis_process run  \
native:zonalstatisticsfb --distance_units=meters --area_units=m2 --ellipsoid=EPSG:7019  \
--INPUT=${lakes} --INPUT_RASTER=${dem} --RASTER_BAND=1 --COLUMN_PREFIX=_ --STATISTICS=3 \
--OUTPUT=lakes_zonal.gpkg

# Set buffer size depending on scale
if [ $scale = 25 ]
then
    bufferDistance=8
elif [ $scale = 50 ]
then
    bufferDistance=16
fi

# Buffer lake polygons
/Applications/QGIS.app/Contents/MacOS/qgis_process.app/Contents/MacOS/qgis_process run \
native:buffer --distance_units=meters --area_units=m2 --ellipsoid=EPSG:7019 \
--INPUT='lakes_zonal.gpkg|layername=lakes_zonal' \
--DISTANCE=$bufferDistance --SEGMENTS=5 --END_CAP_STYLE=0 --JOIN_STYLE=0 --MITER_LIMIT=2 --DISSOLVE=false --OUTPUT=lakes_zonal_buf.gpkg

# Restrict lake elevation values to be as far as possible from contour interval values, save value to column roundedElevation
/Applications/QGIS.app/Contents/MacOS/qgis_process.app/Contents/MacOS/qgis_process run \
native:fieldcalculator --distance_units=meters --area_units=m2 --ellipsoid=EPSG:7019 \
--INPUT='lakes_zonal_buf.gpkg|layername=lakes_zonal_buf' --FIELD_NAME=roundedElevation --FIELD_TYPE=0 --FIELD_LENGTH=0 \
--FIELD_PRECISION=0 --FORMULA="(round(((\"_median\"-${halfInterval})/${desired_interval}))*${desired_interval})+${halfInterval}" \
--OUTPUT=lakes_zonal_rounded.gpkg

unset PROJ_LIB

# Burn a raster with our buffered polygons with modified elevation values using same extent and pixel size as the input DEM
gdal_rasterize -te $EXTENT -a "roundedElevation" -a_nodata 0 -tr $pixelsize $pixelsize lakes_zonal_rounded.gpkg lakes_zonal.tif

# Merge the new raster with the original
gdal_merge.py $dem lakes_zonal.tif -o dem_waterfix.tif
listgeo -tfw dem_waterfix.tif

rm lakes_zonal*

Happy mapping! Till next time.