Add terrain (hillshade and contour) generation scripts

terrain/Dockerfile

@@ -0,0 +1,26 @@
+FROM debian:bookworm
+
+RUN apt update
+RUN apt install -y python3 python3-pip python3-venv libsqlite3-dev gdal-bin git jq
+RUN ln -s /usr/bin/python3 /usr/bin/python
+
+RUN mkdir /root/src
+WORKDIR /root/src
+
+RUN git clone https://github.com/felt/tippecanoe.git
+WORKDIR /root/src/tippecanoe
+RUN make -j && make install
+
+WORKDIR /root/src
+RUN python -m venv .
+RUN bin/pip install mercantile shapely
+
+COPY generate.sh .
+COPY make_tile.py .
+COPY mapzen-dem.wms.xml .
+COPY filter_tiles.py .
+COPY land.geojson .
+COPY entrypoint.sh .
+RUN chmod +x entrypoint.sh
+
+CMD /root/src/entrypoint.sh

terrain/README.md

@@ -0,0 +1,25 @@
+The scripts in this directory can be used to generate the hillshade and contour line tilesets used by OpenTrailMap.
+
+The main entrypoint is `generate.sh` which builds all tiles in both tilesets. Edit the variables at the top of this script to change the min and max zoom and the bounding box for which tiles are built.
+
+Internally, `generate.sh` calls `make_tile.py` repeatedly. `make_tile.py` takes a tile ID, fetches the required DEM data for that region, and builds a hillshade tile and contour tile. There are variables at the top of this file that can be changed to configure how tiles are built (e.g. contour intervals and simplification threshold).
+
+DEM data is fetched from Mapzen's global DEM dataset, hosted on AWS. The file `mapzen-dem.wms.xml` describes how to read this dataset in a format that GDAL tools such as `gdalwarp` understand. `make_tile.py` uses this to fetch the DEM data for a single-tile region (plus required buffer) on the fly, and discards the data once it builds the output tiles. This reduces the disk space required to build tiles.
+
+The scripts have the following dependencies:
+- gdal
+- tippecanoe
+- mercantile
+- jq
+- python
+
+You can install these dependencies yourself, or use the provided Dockerfile to build a container that has them all. When building in Docker, be sure to mount a directory from the host system to `/root/src/build` in the container, as this is where the output tiles will be written.
+
+Running `generate.sh` will produce two tilesets in the working directory called `contours` and `hillshade`, which are directory trees. If desired, these directory trees can be converted to MBTiles or PMTiles archives using a variety of tools, but this is not handled automatically by the scripts in this repository.
+
+Tip: building the tiles creates a lot of intermediate artifacts on disk. These are created in temporary directories and automatically cleaned up afterwards. But you can speed up the process significantly by ensuring that they are written to an in-memory 'ramfs' rather than disk. By default tempdirs are created in `/tmp`, but if that directory isn't a ramfs on your machine (e.g. macOS), you can set the `TMPDIR` environment variable to somewhere else before running `generate.sh`, and tempdirs will be created there instead.
+
+The current tilesets (s3://osmus-tile/hillshade.pmtiles and s3://osmus-tile/contours-feet.pmtiles) were built on 2024-09-07, on a c7g.8xlarge EC2 instance in AWS. The tiled area covers the bbox `[-180, -60, 180, 85]`, which includes all land except Antarctica. The `filter_tiles.py` script was used to additionally exclude tiles in the oceans (i.e. not intersecting `land.geojson`). Hillshade tiles were built for zooms 1-12 and contour lines were built for zooms 8-12.
+
+The build took around 36 hours and cost about $50. It required about 128GB of disk space, though I recommend allocating 256GB or even 512GB since at then end you need to package the tiles into PMTiles, and I found the easiest way to do that was to create MBTiles first, so there were large intermediate artifacts to deal with. Converting the tiles from a directory tree to a PMTiles archive is I/O bound and does not require multiple cores, so you can downsize the EC2 instance to something like a c7g.medium to save some money during this phase.
+

terrain/entrypoint.sh

@@ -0,0 +1,3 @@
+#!/bin/bash
+source bin/activate
+./generate.sh

terrain/filter_tiles.py

@@ -0,0 +1,29 @@
+"""
+Read tile IDs from STDIN, filter them to just those that intersect a GeoJSON geometry,
+and write just the intersecting ones to STDOUT.
+
+Useful in conjunction with `mercantile tiles` command, to filter tiles to an area of
+interest and avoid needlessly tiling unwanted areas e.g. oceans.
+
+Usage: python filter_tiles.py GEOJSON_FILE 
+"""
+
+import json
+import sys
+
+import mercantile
+import shapely
+
+geojson_file = sys.argv[1]
+with open(geojson_file) as f:
+    geojson = shapely.from_geojson(f.read())
+    shapely.prepare(geojson)
+
+for line in sys.stdin:
+    x, y, z = json.loads(line)
+    bounds = mercantile.bounds((x, y, z))
+    if shapely.intersects(geojson, shapely.box(*bounds)):
+        print(json.dumps([x, y, z]), file=sys.stdout)
+
+
+

terrain/generate.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# This script generates hillshade and contour line tiles from Mapzen's DEM dataset
+
+### Configuration parameters
+
+SOURCE_DATASET=mapzen-dem.wms.xml
+MINZOOM=10
+MAXZOOM=12
+# BOUNDS="[-124.73460, 45.56631, -116.87882, 48.99251]" # US-WA
+# BOUNDS="[-180, -85.05, 180, 85.05]" # let's GOOOO!
+BOUNDS="[-180, -60, 180, 85.05]" # exclude antarctica
+
+# Z7 + Z8 took about 9h on my laptop. But osmx extract for a planet file was running at the same time
+
+
+### GDAL settings (these make fetching and reprojecting data faster)
+
+export GDAL_CACHEMAX=4096 # in MB
+export GDAL_MAX_CONNECTIONS=8 # make GDAL use more parallelism when fetching (default is 2)
+
+for zoom in $(seq $MINZOOM $MAXZOOM); do
+  echo $BOUNDS \
+  | mercantile tiles $zoom \
+  | python filter_tiles.py land.geojson \
+  | jq -r '[.[2], .[0], .[1]] | join(" ")' \
+  | xargs --verbose -d '\n' -n 1 -P 16 sh -c "python make_tile.py $SOURCE_DATASET \$0 >/dev/null 2>&1"
+done
+

terrain/make_tile.py

@@ -0,0 +1,173 @@
+"""
+Handles building a contour and hillshade tile for a single Z/X/Y tile ID.
+
+Usage: python make_tile.py SOURCE_DEM Z X Y
+"""
+
+import json
+import os
+import sys
+import tempfile
+import subprocess
+import shutil
+
+# non-stdlib dependencies
+import mercantile
+
+OUTPUT_DIR = "build"
+
+# buffer controls how much margin around the tile region will be added when fetching
+# DEM data. a one-pixel margin is required for hillshades to avoid edge artifacts.
+# this margin also affects how far the contour line vectors extend beyond the edges
+# of the MVT tile. having a bit of margin here helps avoid both rendering artifacts
+# and issues where line simplification causes contour lines to be disconnected at
+# tile boundaries
+BUFFER = 0.01
+
+# controls the width/height of the hillshade tiles. also affects the size of DEM
+# data that is fetched to build the tile. these could be decoupled but in practice
+# it makes sense to have them be the same (there's no reason to fetch more pixels of
+# DEM data, it doesn't improve the hillshade or contour line quality).
+RESOLUTION = 512
+
+# controls compression of hillshade JPEG images; 0-100, higher means bigger filesize
+# but fewer artifacts in the images.
+JPEG_QUALITY = 75
+
+# controls the vertical exaggeration used for the hillshading. exaggeration is
+# necessary to make mountains visible at very low zooms, because when considered
+# on this large scale, Earth's surface is very flat.
+#
+# there's no magic reason for the formula given below, it was tweaked until the
+# hillshading seemed to have relatively uniform contrast at different zoom levels.
+VERTICAL_EXAGGERATION = { z: 2 ** max(0.7 * (8 - z), 1.0) for z in range(1, 13) }
+#     1: 21.0,
+#     2: 13.0,
+#     3: 8.0,
+#     4: 5.0,
+#     5: 4.0,
+#     6: 3.0,
+#     7: 2.0,
+#     8: 1.75,
+#     9: 1.5,
+#     10: 1.25,
+#     11: 1.0,
+#     12: 1.0,
+# }
+
+# sets the contour interval for each zoom level.
+# only zoom levels in this dict will have contour tiles built! if make_tile.py is
+# called with a zoom level that's not in this dict, it will only build a hillshade
+# tile (useful for low zooms where contours aren't needed).
+CONTOUR_INTERVALS = {
+    # zoom level : (major, minor) contour interval
+    8: (4000, 800),
+    9: (2000, 400),
+    10: (2000, 400),
+    11: (800, 160),
+    12: (800, 160),
+}
+
+SIMPLIFICATION_TOLERANCE = 3
+
+# Read CLI args
+source_dem = os.path.abspath(sys.argv[1])
+z, x, y = [int(v) for v in sys.argv[2:]]
+
+# Compute both an exact bbox and buffered bbox for the given tile (in EPSG:3857)
+bounds = mercantile.xy_bounds(x, y, z)
+size = abs(bounds[2] - bounds[0])
+shift = size * BUFFER
+buffered_bounds = mercantile.Bbox(bounds.left - shift, bounds.bottom - shift, bounds.right + shift, bounds.top + shift)
+buffered_resolution = round(RESOLUTION * (1.0 + BUFFER))
+
+def run(argv, *args, **kwargs):
+    """Wrapper around subprocess.run() which prints its arguments first"""
+    print(" ".join(argv), file=sys.stderr)
+    return subprocess.run(argv, *args, **kwargs)
+
+# Make a temporary directory in which all intermediate files are placed.
+# It'll be deleted when the `with` block closes, so be sure to move final
+# artifacts out first.
+with tempfile.TemporaryDirectory() as tmpdir:
+    old_cwd = os.getcwd()
+    os.chdir(tmpdir)
+    print(tmpdir)
+
+    ### Build the hillshade raster tile ###
+
+    # fetch the DEM data for the tile
+    run([
+        "gdalwarp",
+        "-te_srs", "EPSG:3857", "-te", *[str(v) for v in buffered_bounds],
+        "-ts", str(buffered_resolution), str(buffered_resolution),
+        "-r", "cubic",
+        "-overwrite", source_dem, "dem.tif"
+     ])
+
+    # apply vertical exaggeration (for hillshade only)
+    exaggeration = VERTICAL_EXAGGERATION[z]
+    run(["gdal_translate", "-scale", *[str(v) for v in [0, 1, 0, exaggeration]], "dem.tif", "dem.exaggerated.tif"])
+
+    # create a hillshade from the fetched data
+    run(["gdaldem", "hillshade", "-igor", "dem.exaggerated.tif", "hillshade.tif"])
+
+    # trim the margin from the hillshade tile and convert to JPEG
+    run([
+        "gdalwarp",
+        "-te_srs", "EPSG:3857", "-te", *[str(v) for v in bounds],
+        "-ts", str(RESOLUTION), str(RESOLUTION),
+        "-r", "cubic",
+        "-co", f"JPEG_QUALITY={JPEG_QUALITY}",
+        "hillshade.tif", "hillshade.jpg"
+    ])
+
+    # put the hillshade tile in the output directory
+    os.makedirs(os.path.join(old_cwd, OUTPUT_DIR, "hillshade", str(z), str(x)), exist_ok=True)
+    shutil.move("hillshade.jpg", os.path.join(old_cwd, OUTPUT_DIR, "hillshade", str(z), str(x), str(y) + ".jpg"))
+
+    ### Build the contour vector tile ###
+
+    if z in CONTOUR_INTERVALS:
+        major_interval, minor_interval = CONTOUR_INTERVALS[z]
+
+        # create a version of the DEM in WGS 84 (EPGS:4326)
+        # note: this needs to be real, not a VRT, because gdal_contour doesn't
+        # understand scale factors in VRT datasets.
+        run(["gdalwarp", "-t_srs", "EPSG:4326", "-r", "cubic", "dem.tif", "dem.4326.vrt"])
+
+        # input DEM is in meters but we need feet, so apply a scale factor
+        run(["gdal_translate", "-scale", *[str(v) for v in [0, 0.3048, 0, 1]], "dem.4326.vrt", "dem.4326.ft.tif"])
+
+        # compute the contours
+        run(["gdal_contour", "-a", "ele", "-i", str(minor_interval), "dem.4326.ft.tif", "contours-raw.geojson"])
+        # add 'idx' bool property and filter out contours at or below sea level
+        res = run([
+            "jq", "-c",
+            f"(.features[].properties |= (.idx = (fmod(.ele; {major_interval}) == 0))) | .features[] |= select(.properties.ele > 0)",
+            "contours-raw.geojson",
+        ], capture_output=True)
+        # jq can only write to stdout, so manually write its output to a file
+        with open("contours.geojson", "wb") as f:
+            f.write(res.stdout)
+
+        simplification_tolerance = 1 if z == max(CONTOUR_INTERVALS.keys()) else SIMPLIFICATION_TOLERANCE
+        run([
+            "tippecanoe",
+            "-e", "contours", # output directory name
+            "-l", "contours", # layer name in encoded MVT
+            "--minimum-zoom", str(z), "--maximum-zoom", str(z),
+            "--buffer", "1", # units are "screen pixels" i.e. 1/256 of tile width
+            "--simplification", str(simplification_tolerance),
+            "-y", "ele", "-y", "idx", # include only these two attributes
+            "--no-tile-compression",
+            "contours.geojson"
+        ])
+
+        # put the contour tile in the output directory
+        os.makedirs(os.path.join(old_cwd, OUTPUT_DIR, "contours", str(z), str(x)), exist_ok=True)
+        shutil.move(
+            os.path.join("contours", str(z), str(x), str(y) + ".pbf"),
+            os.path.join(old_cwd, OUTPUT_DIR, "contours", str(z), str(x), str(y) + ".mvt"),
+        )
+

terrain/mapzen-dem.wms.xml

@@ -0,0 +1,25 @@
+<GDAL_WMS>
+  <Service name="TMS">
+    <ServerUrl>https://s3.amazonaws.com/elevation-tiles-prod/geotiff/${z}/${x}/${y}.tif</ServerUrl>
+  </Service>
+  <DataWindow>
+    <UpperLeftX>-20037508.34</UpperLeftX>
+    <UpperLeftY>20037508.34</UpperLeftY>
+    <LowerRightX>20037508.34</LowerRightX>
+    <LowerRightY>-20037508.34</LowerRightY>
+    <TileLevel>13</TileLevel>
+    <TileCountX>1</TileCountX>
+    <TileCountY>1</TileCountY>
+    <YOrigin>top</YOrigin>
+  </DataWindow>
+  <Projection>EPSG:3857</Projection>
+  <BlockSizeX>512</BlockSizeX>
+  <BlockSizeY>512</BlockSizeY>
+  <BandsCount>1</BandsCount>
+  <DataType>Float32</DataType>
+  <ZeroBlockHttpCodes>403,404</ZeroBlockHttpCodes>
+  <DataValues>
+    <NoDataValue>-3.4028234663852886e+38</NoDataValue>
+  </DataValues>
+  <Cache></Cache>
+</GDAL_WMS>