contours

contours/Makefile

@@ -0,0 +1,46 @@
+BOUNDS = xmin=582700 ymin=6060750 xmax=584830 ymax=6062750
+XYZ = $(patsubst %.zip,%.xyz,$(wildcard DTM_*.zip))
+SORT = sort -n -k2 -k1 -t,
+
+.PHONY: all
+all: smooth_2_5.gpkg smooth_5.gpkg
+
+smooth_%.gpkg: db/smooth_%
+	ogr2ogr $@ "PG:host=127.0.0.1 user=osm dbname=osm" $(basename $@)
+
+db/smooth_%: db/contour_% chaikin.sql
+	./managedb -- --echo-all \
+		-v ON_ERROR_STOP=1 \
+		-v src=$(notdir $(basename $<)) \
+		-v tbl=$(notdir $(basename $@)) \
+		-f chaikin.sql
+	touch $@
+
+db/contour_%: contour_%.gpkg db/.ready
+	./managedb -- -c "DROP TABLE IF EXISTS $(basename $<)"
+	ogr2ogr -f PostgreSQL "PG:host=127.0.0.1 user=osm dbname=osm" $<
+	touch $@
+
+contour_%.gpkg: all.tif
+	gdal_contour -nln $(basename $@) -i $(subst _,.,$*) -a z $^ $@
+
+all.tif: all.xyz
+	gdal_translate $< $@ \
+		-ot Float32 -a_srs EPSG:3346 \
+		-co COMPRESS=DEFLATE -co PREDICTOR=2
+
+.INTERMEDIATE: all.xyz
+all.xyz: $(XYZ)
+	echo $(XYZ)
+	$(SORT) -m $^ > $@
+
+.INTERMEDIATE: $(XYZ)
+%.xyz: %.zip
+	unzip -qq -c $< $@ | \
+		./clip.awk $(addprefix -v ,$(BOUNDS)) | \
+		$(SORT) > $@
+
+db/.ready: managedb
+	-./managedb stop; rm -fr db
+	./managedb init
+	touch $@

contours/README.md

@@ -0,0 +1,23 @@
+Countour line generator from LIDAR data
+---------------------------------------
+
+Usage:
+
+1. Download contour lines in format `6062999.75,584000.75,88.07`: coordinate
+   pair and height (in meters) to `DTM_XX_YY_N.zip` files.
+2. Adjust `BOUNDS` in the Makefile.
+3. Run `make -j$(nproc) smooth_<X>.gpkg`. This will output a geo-package with
+   contour lines every `X` meters. X can be fractional (e.g.
+   `smooth_2.5.gpkg`).
+
+Dependencies:
+
+- postgis
+- gdal-bin
+- unzip
+
+TODO:
+
+1. Replace `managedb` stack with postgis-over-docker. This will lead to a less
+   fragile postgis setup.
+2. Accept `BOUNDS` in a way that does not require to change the Makefile.

contours/chaikin.sql

@@ -0,0 +1,20 @@
+DROP TABLE IF EXISTS :tbl;
+
+CREATE TABLE :tbl (
+    fid serial NOT NULL,
+    z DOUBLE PRECISION,
+    geom geometry(MULTILINESTRING, 3346)
+);
+
+INSERT INTO :tbl (
+    SELECT
+        fid,
+        z,
+        ST_Multi(
+            ST_ChaikinSmoothing (
+                ST_SimplifyVW(geom, 50),
+                5
+            )
+        ) AS geoms
+    FROM
+        :src);

contours/clip.awk

@@ -0,0 +1,8 @@
+#!/usr/bin/awk -f
+BEGIN {
+    FS = ","
+}
+
+$1 > ymin && $1 < ymax && $2 > xmin && $2 < xmax {
+    print $2 "," $1 "," $3
+}

contours/managedb

@@ -0,0 +1,35 @@
+#!/bin/bash
+
+PATH=$PATH:/usr/lib/postgresql/12/bin
+case ${1:-} in
+    init)
+        mkdir -p db && initdb db
+        mkdir -p db/wal
+        sed -i "s/.*unix_socket_dir.*/unix_socket_directories = '.\/wal'/" \
+            db/postgresql.conf
+        pg_ctl -D db -l db/logfile start
+
+        export PGHOST=127.0.0.1
+        psql postgres -c 'CREATE ROLE osm WITH SUPERUSER LOGIN'
+        psql postgres -c 'CREATE DATABASE osm'
+        psql osm osm -c 'CREATE EXTENSION postgis'
+        ;;
+    start)
+        pg_ctl -D db -l db/logfile start
+        ;;
+    stop)
+        pg_ctl -D db -l db/logfile stop
+        ;;
+    "" | --)
+        [[ $# -gt 1 ]] && shift
+        exec env \
+            PGHOST=127.0.0.1 \
+            PGUSER=osm \
+            PGDATABASE=osm \
+            psql "$@"
+        ;;
+    *)
+        >&2 echo "Unknown command: '$*'"
+        exit 1
+        ;;
+esac