Download topography without GMT

[boriskaus]

Why

import_topo lived in the GMT extension, so getting topography meant taking on GMT — a heavy dependency for what is really one call. Two problems follow from it:

• It downgrades other packages. GMT needs GDAL_jll, whose HDF5_jll bounds ([3.0.3 – 300.500.101] or [305.1300.100 – 305.1300.301]) exclude the 2.2.x that PETSc_jll requires. Installing GMT next to LaMEM therefore walks LaMEM back a minor version — 0.5.3 → 0.4.16 — without saying why.
• GMT's share directory is not reliably relocatable. libgmt has the BinaryBuilder sandbox path /workspace/destdir/share compiled in, and accepts it on nothing more than access(). Anyone who has a real /workspace/destdir (Yggdrasil contributors, for instance) gets File /workspace/destdir/share/google/turbo.cpt was not found and a failed precompile.

What

The grids are served over plain HTTP, so the only thing GMT was really needed for was decoding them. OpenJpeg_jll does that, and its tree — JpegTurbo, Libtiff, LittleCMS, libpng, zlib, zstd, XZ, LERC — shares nothing with GDAL or HDF5. So import_topo moves out of the GMT extension and into the package, with OpenJpeg_jll and Scratch as hard dependencies. import_GeoTIFF stays in the extension, so GMT remains an optional dependency for it.

julia> Topo = import_topo([-18.7, -17.1, 28.0, 29.2], res="03s")     # no `using GMT`
julia> Topo = import_topo(lon=[4,20], lat=[37,49])                   # unchanged
julia> Topo = import_topo(lims, file="@earth_relief_15s")            # still accepted
julia> Topo = import_topo(lims, dataset="earth_gebco", res="15s")

Tiles are cached in a scratch space, so the same region is not downloaded twice.

What the server does, and what had to be read rather than assumed

Everything that varies between sets comes from gmt_data_server.txt, not from hardcoded tables:

• tile size varies: 30° for 01m, 60° for 02m, 90° for 03m, 180° for 04m/05m, 15° for 30s, 10° for 15s, 1° for 03s;
• the integers are scaled — earth_relief_15s stores half metres, so reading them as metres doubles every elevation;
• registration is not free: SRTM is gridline-only, 15s is pixel-only, and gridline tiles repeat the edge they share with the next tile, so stacking them naively duplicates a row;
• latitude tiles are anchored at the south pole, not at multiples of the tile size — that is what broke 02m–05m;
• SRTM covers land only: it writes a flat zero for the sea and omits all-ocean tiles entirely (14297 tiles, not 64800). Both are filled from the coarser set the server names for the purpose, which is why GMT's 03s has bathymetry at all.

Verification against GMT

Over La Palma ([-18.7, -17.1, 28.0, 29.2]), comparing with gmtread directly:

Resolution	earth_relief	earth_gebco	earth_synbath
01m	0.0 m	0.0 m	—
02m, 03m, 04m, 05m	0.0 m	—	—
30s	0.0 m	—	—
15s	0.0 m	0.0 m	0.0 m
03s	exact on land; 0.68 m mean in the sea	—	—

Also 0.0 m for the Alps at 15s and Iceland at 30s, and the existing suite's expected value sum(Topo.depth.val) = 1076.7045 is reproduced exactly, so this is drop-in.

The one remaining difference is 03s over water, where GMT's interpolation of the 15s filler differs from the bilinear one here: max 251 m, mean 0.68 m, with 0.6% of cells more than 10 m apart, all at the coastline. On land it is exact — the summit lands at −17.885, 28.754, which is Roque de los Muchachos.

Tests

test/test_import_topo.jl: tile naming including the SRTM special case, pixel-vs-gridline sizes, the sets that exist in one registration only, the file= spelling, georeferencing (the summit position catches a transposed or mis-assembled grid), that the sea is not flat, and that the keyword and vector forms agree. 27/27 locally.

Compatibility with existing callers

Every call form already in the repo, the docs and the tutorials is accepted unchanged:
import_topo([8,9,50,51]), the lat=/lon= keywords, file="@earth_relief_01m", the
.grd/.nc suffix the tutorials append, and the maxattempts that topo_prefetch.jl
passes (now honoured as the number of times each mirror is tried). test_GMT.jl,
test_WaterFlow.jl and topo_prefetch.jl all pass against this implementation.

Note

Coarse resolutions (06m and up) are single NetCDF files rather than tiles, and reading those still needs NCDatasets — which stays a weak dependency. The tiled resolutions, which are the interesting ones, need nothing beyond the new hard dependencies.

🤖 Generated with Claude Code

[github-actions]

Your PR requires formatting changes to meet the project's style guidelines.
Please consider running Runic (git runic main) to apply these changes.

Click here to view the suggested changes.

diff --git a/src/import_topo.jl b/src/import_topo.jl
index b953464..eef9b29 100644
--- a/src/import_topo.jl
+++ b/src/import_topo.jl
@@ -57,11 +57,13 @@ function server_entry(dataset::AbstractString, res::AbstractString, reg::Abstrac
         length(f) < 11 && continue
         name = rstrip(f[2], ['/'])
         replace(name, ".grd" => "") == want || continue
-        return (name = name,
-                scale = something(tryparse(Float64, f[5]), 1.0),
-                offset = something(tryparse(Float64, f[6]), 0.0),
-                tile_deg = something(tryparse(Int, f[8]), 0),
-                filler = strip(f[11]))
+        return (
+            name = name,
+            scale = something(tryparse(Float64, f[5]), 1.0),
+            offset = something(tryparse(Float64, f[6]), 0.0),
+            tile_deg = something(tryparse(Int, f[8]), 0),
+            filler = strip(f[11]),
+        )
     end
     return nothing
 end
@@ -110,8 +112,10 @@ The name of the tile whose south-west corner is at (`lat`, `lon`).
 The 3 arcsecond set is the odd one out: it is SRTM repackaged, and is named after that rather
 than after the dataset it is served under.
 """
-function tile_filename(lat::Integer, lon::Integer, dataset::AbstractString,
-                       res::AbstractString, reg::AbstractString)
+function tile_filename(
+        lat::Integer, lon::Integer, dataset::AbstractString,
+        res::AbstractString, reg::AbstractString
+    )
     ns = lat < 0 ? "S" : "N"
     ew = lon < 0 ? "W" : "E"
     corner = string(ns, lpad(abs(lat), 2, '0'), ew, lpad(abs(lon), 3, '0'))
@@ -222,13 +226,15 @@ Read a global `.grd` grid, returning `(lon, lat, z)`. Filled in by the `NCDatase
 extension; the tiled resolutions do not go through here.
 """
 function read_topo_netcdf(file::AbstractString)
-    return error("""
+    return error(
+        """
         reading the global grid $(basename(file)) needs NCDatasets. Either
 
             using NCDatasets
 
         or ask for one of the tiled resolutions (05m and finer), which need nothing extra.
-        """)
+        """
+    )
 end
 
 """
@@ -238,8 +244,10 @@ Fetch one file, trying each mirror in turn. Returns `true` if it was downloaded,
 the file is not on the server -- which is not an error: the 3 arcsecond set covers land only,
 so a tile that is entirely ocean simply does not exist.
 """
-function download_tile(relative_url::AbstractString, dest::AbstractString;
-                      maxattempts::Integer = 5, timeout::Real = TILE_TIMEOUT)
+function download_tile(
+        relative_url::AbstractString, dest::AbstractString;
+        maxattempts::Integer = 5, timeout::Real = TILE_TIMEOUT
+    )
     isfile(dest) && filesize(dest) > 0 && return true       # cached
 
     tmp = dest * ".part"                                    # never leave a half file behind
@@ -290,8 +298,10 @@ end
 Assemble the tiles that cover `limits` and cut out the requested region. Tiles that are not
 on the server are taken to be at sea level, which is what the land-only 3 arcsecond set needs.
 """
-function grid_from_tiles(limits, dataset::AbstractString, res::AbstractString, reg::AbstractString;
-                         maxattempts::Integer = 5)
+function grid_from_tiles(
+        limits, dataset::AbstractString, res::AbstractString, reg::AbstractString;
+        maxattempts::Integer = 5
+    )
     lonmin, lonmax, latmin, latmax = limits
     # Not every set is published in both registrations: the SRTM-derived ones are gridline
     # only, 15s is pixel only. GMT serves whichever the set actually has rather than
@@ -317,8 +327,10 @@ function grid_from_tiles(limits, dataset::AbstractString, res::AbstractString, r
     # tile contributes `step` new rows/columns and only the last one of a row contributes
     # its final edge as well.
     step = tile_reg == "g" ? n - 1 : n
-    Z = zeros(Int16, length(lats) * step + (tile_reg == "g"),
-                     length(lons) * step + (tile_reg == "g"))
+    Z = zeros(
+        Int16, length(lats) * step + (tile_reg == "g"),
+        length(lons) * step + (tile_reg == "g")
+    )
     dir = topo_cache_dir()
     sub = tile_directory(dataset, res, tile_reg)
 
@@ -331,8 +343,10 @@ function grid_from_tiles(limits, dataset::AbstractString, res::AbstractString, r
             continue
         end
         A = decode_jp2(file, n)
-        Z[((j - 1) * step + 1):((j - 1) * step + n),
-          ((i - 1) * step + 1):((i - 1) * step + n)] = A[end:-1:1, :]
+        Z[
+            ((j - 1) * step + 1):((j - 1) * step + n),
+            ((i - 1) * step + 1):((i - 1) * step + n),
+        ] = A[end:-1:1, :]
     end
 
     # the coordinates of each sample: cell centres for pixel registration, cell corners
@@ -374,8 +388,10 @@ set that the server names for the purpose.
 function fill_gaps!(Z, lon, lat, filler)
     any(iszero, Z) || return Z
     fdataset, fres, freg = filler
-    flon, flat, FZ = grid_from_tiles([minimum(lon), maximum(lon), minimum(lat), maximum(lat)],
-                                     fdataset, fres, freg)
+    flon, flat, FZ = grid_from_tiles(
+        [minimum(lon), maximum(lon), minimum(lat), maximum(lat)],
+        fdataset, fres, freg
+    )
     for j in axes(Z, 1), i in axes(Z, 2)
         iszero(Z[j, i]) || continue
         Z[j, i] = bilinear(flon, flat, FZ, lon[i], lat[j])
@@ -402,7 +418,7 @@ function bilinear(xs::AbstractVector, ys::AbstractVector, A::AbstractMatrix, x::
     ty = clamp(ty, 0.0, 1.0)
 
     return (1 - tx) * (1 - ty) * A[j, i] + tx * (1 - ty) * A[j, i + 1] +
-           (1 - tx) * ty * A[j + 1, i] + tx * ty * A[j + 1, i + 1]
+        (1 - tx) * ty * A[j + 1, i] + tx * ty * A[j + 1, i + 1]
 end
 
 """
@@ -412,8 +428,10 @@ The coarse resolutions are published as one NetCDF file for the whole globe, whi
 enough to download and cut down here. Reading it needs `NCDatasets`, which is a weak
 dependency -- the tiled resolutions, which are the interesting ones, need nothing extra.
 """
-function grid_from_single_file(limits, dataset::AbstractString, res::AbstractString, reg::AbstractString;
-                               maxattempts::Integer = 5)
+function grid_from_single_file(
+        limits, dataset::AbstractString, res::AbstractString, reg::AbstractString;
+        maxattempts::Integer = 5
+    )
     lonmin, lonmax, latmin, latmax = limits
     name = "$(dataset)_$(res)_$(reg).grd"
     file = joinpath(topo_cache_dir(), name)
@@ -459,10 +477,12 @@ julia> Topo = import_topo(lon=[4,20], lat=[37,49])
 
 Tiles are cached, so importing the same region again does not download them a second time.
 """
-function import_topo(limits; dataset::AbstractString = "earth_relief",
-                     res::AbstractString = "01m", reg::AbstractString = "g",
-                     file::Union{Nothing, AbstractString} = nothing,
-                     maxattempts::Integer = 5)
+function import_topo(
+        limits; dataset::AbstractString = "earth_relief",
+        res::AbstractString = "01m", reg::AbstractString = "g",
+        file::Union{Nothing, AbstractString} = nothing,
+        maxattempts::Integer = 5
+    )
 
     # accept the `file="@earth_relief_01m"` spelling that this function used to take
     if !isnothing(file)
@@ -505,6 +525,6 @@ function parse_topo_file(file::AbstractString)
     m = match(r"^(.*)_(\d+[dms])(?:_([gp]))?$", s)
     isnothing(m) && error("cannot make sense of the topography file name \"$file\"")
     # a name without a trailing _g/_p leaves the registration open, and GMT then serves
-            # whichever the set has, preferring gridline -- the same default as `reg`
-            return String(m[1]), String(m[2]), isnothing(m[3]) ? "g" : String(m[3])
+    # whichever the set has, preferring gridline -- the same default as `reg`
+    return String(m[1]), String(m[2]), isnothing(m[3]) ? "g" : String(m[3])
 end

[mthielma]

Looks good and should resolve the recurring GMT issues.