From c9e71ccf74f79cff197c9667cc0f90d06727b0f4 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 29 Sep 2026 09:42:28 +0200 Subject: [PATCH 1/3] Download topography without GMT `import_topo` needed `GMT`, which is a heavy dependency for one call: it pulls in GDAL_jll, whose HDF5 bounds clash with the PETSc stack (installing GMT alongside LaMEM silently downgrades LaMEM by a minor version), and GMT's own share directory is not reliably relocatable. The grids it fetches are served over plain HTTP, so the only thing GMT was really needed for was decoding them. Do that with `OpenJpeg_jll`, whose tree is small and shares nothing with GDAL. `import_topo` therefore moves out of the GMT extension and into the package itself; `import_GeoTIFF` stays behind, so GMT remains an optional dependency for that. Everything that varies between the sets is read from the server's own index rather than assumed: - tiles are 30, 60, 90, 180, 15, 10 or 1 degrees wide depending on the set; - the integers in them are scaled (earth_relief_15s stores half metres); - some sets exist in only one registration -- SRTM is gridline, 15s is pixel -- and gridline tiles repeat the edge they share with the next tile; - latitude tiles are anchored at the south pole, not at multiples of their size; - SRTM covers land only, writing a flat zero for the sea and omitting all-ocean tiles, so both are filled from the coarser set the server names for it. Checked against GMT over La Palma: every tiled resolution of earth_relief (01m-05m, 30s, 15s) and of earth_gebco and earth_synbath agrees to 0.0 m, as do the Alps and Iceland. 03s agrees exactly on land; in the sea it differs by 0.68 m on average, where GMT's interpolation of the filler differs from ours. The existing test's expected value, sum(Topo.depth.val) = 1076.7045, is reproduced. Coarse resolutions (06m and up) are single NetCDF files and still need `NCDatasets`, which stays a weak dependency. Co-Authored-By: Claude Opus 5 (1M context) --- Project.toml | 4 + ext/GMT_utils.jl | 151 +--------- ext/NCDatasets_utils.jl | 21 +- src/GeophysicalModelGenerator.jl | 10 +- src/import_topo.jl | 482 +++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/test_import_topo.jl | 88 ++++++ 7 files changed, 599 insertions(+), 158 deletions(-) create mode 100644 src/import_topo.jl create mode 100644 test/test_import_topo.jl diff --git a/Project.toml b/Project.toml index 1a1ab5c9e..9bba72bf4 100644 --- a/Project.toml +++ b/Project.toml @@ -20,8 +20,10 @@ LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +OpenJpeg_jll = "643b3616-a352-519d-856d-80112ee9badc" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Scratch = "6c6a2e73-6563-6170-7368-637461726353" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -64,8 +66,10 @@ LightXML = "0.8, 0.9" MeshIO = "0.5" NCDatasets = "0.14" NearestNeighbors = "0.2 - 0.4" +OpenJpeg_jll = "2.5" ParallelTestRunner = "2" Parameters = "0.9 - 0.12" +Scratch = "1" SpecialFunctions = "1.0, 2" StaticArrays = "1" Statistics = "1" diff --git a/ext/GMT_utils.jl b/ext/GMT_utils.jl index f88f69beb..c9c000304 100644 --- a/ext/GMT_utils.jl +++ b/ext/GMT_utils.jl @@ -1,7 +1,7 @@ # NOTE: these are useful routines that are only made available when the GMT package is already loaded in the REPL module GMT_utils -import GeophysicalModelGenerator: import_topo, import_GeoTIFF +import GeophysicalModelGenerator: import_GeoTIFF # We do not check `isdefined(Base, :get_extension)` as recommended since # Julia v1.9.0 does not load package extensions when their dependency is @@ -16,155 +16,6 @@ using GeophysicalModelGenerator: lonlatdepth_grid, GeoData, UTMData, km, remove_ println("Loading GMT routines within GMG") -# Official GMT data server mirrors, see https://www.generic-mapping-tools.org/mirrors. -# `import_topo` cycles through these if a download fails, so that a single -# mirror being offline does not break the download (or the tests that check it). -# Ordered by measured response; "singapore" is omitted as it was unreachable. -const GMT_DATA_SERVERS = ["oceania", "australia", "sdsc-opentopography", "noaa", "brasil", "china"] - - -""" - Topo = import_topo(limits; file::String="@earth_relief_01m", maxattempts=5) - -Uses `GMT` to download the topography of a certain region, specified with limits=[lon_min, lon_max, lat_min, lat_max]. -Sometimes download fails because of the internet connection. We do `maxattempts` to download it. - -Note: -==== -- latitude values in the southern hemisphere should have a minus sign (e.g., -2.8) -- longitude values that are "west" should *either* come with a minus sign *or* are defined by values >180 - -| Dataset | Resolution | Description | -|:---------------- | ------------ | ----------------------------------------------------------- | -| "@earth\\_relief\\_01s" | 1 arc sec | SRTM tiles (14297 tiles, land only, 60S-60N) [NASA/USGS] | -| "@earth\\_relief\\_03s" | 3 arc sec | SRTM tiles (14297 tiles, land only, 60S-60N) [NASA/USGS] | -| "@earth\\_relief\\_15s" | 15 arc sec | SRTM15+ [David Sandwell, SIO/UCSD] | -| "@earth\\_relief\\_30s" | 30 arc sec | SRTM30+ [Becker et al., 2009, SIO/UCSD] | -| "@earth\\_relief\\_01m" | 1 arc min | ETOPO1 Ice surface [NEIC/NOAA] | -| "@earth\\_relief\\_02m" | 2 arc min | ETOPO2v2 Ice surface [NEIC/NOAA] | -| "@earth\\_relief\\_03m" | 3 arc min | ETOPO1 after Gaussian spherical filtering (5.6 km fullwidth)| -| "@earth\\_relief\\_04m" | 4 arc min | ETOPO1 after Gaussian spherical filtering (7.5 km fullwidth)| -| "@earth\\_relief\\_05m" | 5 arc min | ETOPO1 after Gaussian spherical filtering (9 km fullwidth) | -| "@earth\\_relief\\_06m" | 6 arc min | ETOPO1 after Gaussia30n spherical filtering (10 km fullwidth) | -| "@earth\\_relief\\_10m" | 10 arc min | ETOPO1 after Gaussian spherical filtering (18 km fullwidth) | -| "@earth\\_relief\\_15m" | 20 arc min | ETOPO1 after Gaussian spherical filtering (28 km fullwidth) | -| "@earth\\_relief\\_20m" | 20 arc min | ETOPO1 after Gaussian spherical filtering (37 km fullwidth) | -| "@earth\\_relief\\_30m" | 30 arc min | ETOPO1 after Gaussian spherical filtering (55 km fullwidth) | -| "@earth\\_relief\\_60m" | 60 arc min | ETOPO1 after Gaussian spherical filtering (111 km fullwidth)| - -*Note*: this routine is only available once the GMT.jl package is loaded in the REPL - -# Example -```julia-repl -julia> Topo = import_topo([4,20,37,49]); -GeoData - size : (960, 720, 1) - lon ϵ [ 4.0 : 19.983333333333334] - lat ϵ [ 37.0 : 48.983333333333334] - depth ϵ [ -3.8725 km : 4.2495 km] - fields: (:Topography,) -``` -And you can save this to Paraview with -```julia -julia> write_paraview(Topo,"Topo_Alps") -1-element Vector{String}: - "Topo_Alps.vts" -``` -""" -function import_topo(limits; file::String = "@earth_relief_01m", maxattempts = 5) - - # Correct if negative values are given (longitude coordinates that are west) - ind = findall(limits[1:2] .< 0) - - if (limits[1] < 0) && (limits[2] < 0) - limits[ind] .= 360 .+ limits[ind] - limits[1:2] = sort(limits[1:2]) - end - - # Download topo file. If the GMT data server is unreachable, GMT stops - # auto-downloading and reports "Remote download is currently deactivated". - # - # First try the server the user/GMT is already configured with, so that the - # normal case does not touch any global state. Only if that fails do we fall - # back to the official mirrors (https://www.generic-mapping-tools.org/mirrors), - # so a single server being offline does not fail the download. The download - # itself is still exercised -- only *which* mirror serves it changes. - # - # Note: `gmtset` changes GMT's *global* configuration (it writes a gmt.conf in - # the working directory), so switching mirrors is visible to anything else - # using GMT concurrently. We therefore restore the original setting afterwards - # and only switch when a download has actually failed. - local G - local last_err = nothing - - # `gmtset` writes/updates a gmt.conf in the working directory. Remember whether - # one already existed so we can leave the directory as we found it. - conf_file = joinpath(pwd(), "gmt.conf") - had_conf = isfile(conf_file) - conf_backup = had_conf ? read(conf_file) : nothing - - try - for attempt in 1:maxattempts - if attempt > 1 - # Previous attempt failed -> try the next mirror in the list - server = GMT_DATA_SERVERS[mod1(attempt - 1, length(GMT_DATA_SERVERS))] - gmtset(GMT_DATA_SERVER = server) - @info "Retrying GMT topography download from mirror \"$server\"" - end - try - G = gmtread(file, limits = limits, grid = true) - break - catch e - last_err = e - @warn "Failed downloading GMT topography ($file) on attempt $attempt/$maxattempts" - if attempt < maxattempts - sleep(5) # brief pause before trying the next mirror - end - end - end - finally - # Restore the original configuration, so our mirror choice does not leak - # into other code (or other tests) that use GMT afterwards. - try - if had_conf - conf_backup !== nothing && write(conf_file, conf_backup) - elseif isfile(conf_file) - rm(conf_file, force = true) - end - catch - end - end - if !(@isdefined G) - error("Could not download GMT topography \"$file\" after $maxattempts attempts (mirrors tried: $(GMT_DATA_SERVERS)). Last error: $last_err") - end - - # Transfer to GeoData - nx, ny = size(G.z, 2), size(G.z, 1) - Lon, Lat, Depth = lonlatdepth_grid(G.x[1:nx], G.y[1:ny], 0) - @views Depth[:, :, 1] = 1.0e-3 * G.z' - Topo = GeoData(Lon, Lat, Depth, (Topography = Depth * km,)) - - return Topo -end - -""" - import_topo(; lat::Vector{2}, lon::Vector{2}, file::String="@earth_relief_01m", maxattempts=5) - -Imports topography (using GMT), by specifying keywords for latitude and longitude ranges - -# Example -========= -```julia -julia> Topo = import_topo(lat=[30,40], lon=[30, 50] ) -``` -The values can also be defined as tuples: -```julia -julia> Topo = import_topo(lon=(-50, -40), lat=(-10,-5), file="@earth_relief_30s") -``` - -""" -import_topo(; lat = [37, 49], lon = [4, 20], file::String = "@earth_relief_01m", maxattempts = 5) = import_topo([lon[1], lon[2], lat[1], lat[2]], file = file, maxattempts = maxattempts) - """ data_GMT = import_GeoTIFF(fname::String; fieldname=:layer1, negative=false, iskm=true, NorthernHemisphere=true, constantDepth=false, removeNaN_z=false, removeNaN_field=false) diff --git a/ext/NCDatasets_utils.jl b/ext/NCDatasets_utils.jl index 77717f8c5..2b150ee7d 100644 --- a/ext/NCDatasets_utils.jl +++ b/ext/NCDatasets_utils.jl @@ -8,7 +8,7 @@ end using NCDatasets: nc_create, NC_NETCDF4, NC_CLOBBER, NC_NOWRITE, nc_def_dim, nc_def_compound, nc_insert_compound, nc_def_var, nc_put_var, nc_close, NC_INT, nc_unsafe_put_var, libnetcdf, check, ncType, nc_open, nc_inq_vartype, nc_inq_compound_nfields, nc_inq_compound_size, nc_inq_compound_name, nc_inq_compound_fieldoffset, nc_inq_compound_fieldndims, nc_inq_compound_fielddim_sizes, nc_inq_compound_fieldname, nc_inq_compound_fieldindex, nc_inq_compound_fieldtype, nc_inq_compound, nc_inq_varid, nc_get_var!, nc_insert_array_compound, nc_def_vlen -import GeophysicalModelGenerator: write_ASAGI, read_ASAGI, tomo_2_GeoData +import GeophysicalModelGenerator: write_ASAGI, read_ASAGI, tomo_2_GeoData, read_topo_netcdf using GeophysicalModelGenerator: CartData, GeoData, lonlatdepth_grid, xyz_grid """ @@ -221,4 +221,23 @@ function tomo_2_GeoData(filename::String; vel_type::String = "vs") end + +""" + read_topo_netcdf(file) + +Read one of the global `.grd` grids that the GMT data server publishes for its coarse +resolutions, and return `(lon, lat, z)`. The tiled resolutions do not come through here -- +they are JPEG2000 and are decoded without NCDatasets. +""" +function read_topo_netcdf(file::AbstractString) + return NCDataset(file) do ds + lon = Array(ds["lon"][:]) + lat = Array(ds["lat"][:]) + z = Array(ds["z"][:, :]) + # the file stores z as (lon, lat); the rest of the code wants (lat, lon) + (lon, lat, permutedims(coalesce.(z, 0.0f0))) + end +end + + end diff --git a/src/GeophysicalModelGenerator.jl b/src/GeophysicalModelGenerator.jl index f5eee56c4..5f8e5715b 100644 --- a/src/GeophysicalModelGenerator.jl +++ b/src/GeophysicalModelGenerator.jl @@ -49,6 +49,8 @@ include("surface_functions.jl") include("movies_from_pics.jl") include("sea_lvl.jl") include("WaterFlow.jl") +include("import_topo.jl") +export import_topo # Add optional routines (only activated when the packages are loaded) @@ -78,18 +80,12 @@ Reads a seismic tomography NetCDF file and returns a GeoData structure. Requires function tomo_2_GeoData end export tomo_2_GeoData -""" - import_topo -Optional routine that imports topography. It requires you to load `GMT` -""" -function import_topo end - """ import_GeoTIFF Optional routine that imports GeoTIFF images. It requires you to load `GMT` """ function import_GeoTIFF end -export import_topo, import_GeoTIFF +export import_GeoTIFF # GLMakie routines diff --git a/src/import_topo.jl b/src/import_topo.jl new file mode 100644 index 000000000..01856adf8 --- /dev/null +++ b/src/import_topo.jl @@ -0,0 +1,482 @@ +# Downloads topography from the GMT data server without going through GMT itself. +# +# The server publishes its grids in two forms: the coarse resolutions (06m and coarser) as +# single NetCDF files, and the finer ones as tiles of JPEG2000. Both are plain HTTP +# downloads, so the only thing GMT was really needed for was decoding -- which `OpenJpeg_jll` +# does. That matters because GMT pulls in GDAL_jll, whose HDF5 bounds conflict with the +# PETSc-based stack that LaMEM uses, and because GMT's own share directory is not reliably +# relocatable. + +using Downloads, OpenJpeg_jll, Scratch + +""" +The GMT data server mirrors. The first that answers is used; the others are tried in turn if +a download fails, so one server being offline does not fail the import. +""" +const GMT_SERVERS = [ + "https://oceania.generic-mapping-tools.org/server/earth", + "https://data.generic-mapping-tools.org/server/earth", + "https://au.generic-mapping-tools.org/server/earth", + "https://sa.generic-mapping-tools.org/server/earth", +] + +""" + tile_size_deg(res) + +The width and height, in degrees, of one tile at resolution `res`, or `nothing` if that +resolution is published as a single global file rather than as tiles. +""" +function tile_size_deg(dataset::AbstractString, res::AbstractString, reg::AbstractString) + e = server_entry(dataset, res, reg) + isnothing(e) && return nothing # not tiled: one file for the globe + endswith(e.name, ".grd") && return nothing + return e.tile_deg +end + +""" + server_entry(dataset, res, reg) + +The line that `gmt_data_server.txt` holds for this set, as a named tuple. Everything that +varies from set to set -- how large a tile is, how the integers in it scale to metres, which +coarser set fills the gaps -- is published there, so it is read rather than assumed. +""" +function server_entry(dataset::AbstractString, res::AbstractString, reg::AbstractString) + index = server_index_file() + isnothing(index) && return nothing + want = "$(dataset)_$(res)_$(reg)" + for line in eachline(index) + startswith(line, "#") && continue + f = split(line, '\t') + 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])) + end + return nothing +end + +""" + server_index_file() + +The server's own index of its datasets, downloaded once and cached. +""" +function server_index_file() + index = joinpath(topo_cache_dir(), "gmt_data_server.txt") + (isfile(index) && filesize(index) > 0) && return index + for server in GMT_SERVERS + root = replace(server, "/server/earth" => "") + try + Downloads.download("$root/gmt_data_server.txt", index) + return index + catch + end + end + return nothing +end + +""" + tile_size_px(res, deg, reg) + +How many pixels a tile has along each side: the tile width in degrees divided by the spacing +that `res` names (`m` = arcminutes, `s` = arcseconds), plus one for gridline registration. +""" +function tile_size_px(res::AbstractString, deg::Integer, reg::AbstractString) + unit = last(res) + n = parse(Int, res[1:(end - 1)]) + spacing_deg = unit == 'm' ? n / 60 : n / 3600 + npix = round(Int, deg / spacing_deg) + # Pixel registration ("p") samples cell centres, so a tile of `deg` degrees holds exactly + # `npix` of them. Gridline registration ("g") samples the cell corners instead and so + # repeats the shared edge: one more row and column. Getting this wrong shifts everything. + return reg == "g" ? npix + 1 : npix +end + +""" + tile_filename(lat, lon, dataset, res, reg) + +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) + ns = lat < 0 ? "S" : "N" + ew = lon < 0 ? "W" : "E" + corner = string(ns, lpad(abs(lat), 2, '0'), ew, lpad(abs(lon), 3, '0')) + res in ("03s", "01s") && return "$corner.SRTMGL$(res == "03s" ? 3 : 1).jp2" + return "$corner.$(dataset)_$(res)_$(reg).jp2" +end + +""" + tile_directory(dataset, res, reg) + +The directory on the server that holds the tiles of this set. +""" +function tile_directory(dataset::AbstractString, res::AbstractString, reg::AbstractString) + res in ("03s", "01s") && return "earth_relief_$(res)_g" + return "$(dataset)_$(res)_$(reg)" +end + +""" + dataset_scale(dataset, res, reg) + +The scale factor and offset that turn the integers in a tile into metres. The server +publishes these in `gmt_data_server.txt`, one line per set: `earth_relief_15s_p` is stored in +half metres, for instance, while the SRTM sets are stored directly in metres. Reading them +rather than assuming keeps this correct if a set is ever re-encoded. + +Falls back to `(1, 0)` if the index cannot be fetched, which is right for every set that is +not scaled. +""" +function dataset_scale(dataset::AbstractString, res::AbstractString, reg::AbstractString) + e = server_entry(dataset, res, reg) + return isnothing(e) ? (1.0, 0.0) : (e.scale, e.offset) +end + +""" + filler_dataset(dataset, res, reg) + +Some sets cover only part of the globe and are published with a coarser set to fill the rest: +the 3 arcsecond relief is SRTM, which is land only, and the server names +`earth_relief_15s_p` as what belongs in the sea. GMT resamples that filler into the gaps, and +so do we -- otherwise the ocean comes out flat at zero. + +Returns `nothing` when a set needs no filler. +""" +function filler_dataset(dataset::AbstractString, res::AbstractString, reg::AbstractString) + e = server_entry(dataset, res, reg) + isnothing(e) && return nothing + (isempty(e.filler) || e.filler == "-") && return nothing + m = match(r"^(.*)_(\d+[dms])_([gp])$", e.filler) + return isnothing(m) ? nothing : (String(m[1]), String(m[2]), String(m[3])) +end + +""" + native_registration(dataset, res, reg) + +Which registration this set is actually published in. Most come in both, and then the +caller's `reg` is used; the SRTM-derived sets are gridline only and 15s is pixel only, so +asking for the other one would 404 on every tile. +""" +function native_registration(dataset::AbstractString, res::AbstractString, reg::AbstractString) + res in ("03s", "01s") && return "g" + index = joinpath(topo_cache_dir(), "gmt_data_server.txt") + isfile(index) || return reg + have = String[] + for line in eachline(index) + startswith(line, "#") && continue + f = split(line, '\t') + length(f) < 4 && continue + name = replace(rstrip(f[2], ['/']), ".grd" => "") + startswith(name, "$(dataset)_$(res)_") && push!(have, String(strip(f[4]))) + end + isempty(have) && return reg + return reg in have ? reg : first(have) +end + +""" + subset_indices(coords, lo, hi, px) + +The samples that make up `lo..hi`, following what GMT returns. With gridline registration the +samples are the cell corners and GMT widens the request out to the nodes that enclose it; +with pixel registration they are the cell centres and a cell belongs to the region when its +centre does. +""" +function subset_indices(coords::AbstractVector, lo::Real, hi::Real, px::Real, reg::AbstractString) + if reg == "g" + # nodes sit on the cell corners, and GMT widens the request out to the nodes that + # enclose it -- but no further, so the tolerance stops just short of a full step + tol = px - 1.0e-9 + return findall(c -> (lo - tol) < c < (hi + tol), coords) + else + # nodes are cell centres, and a cell belongs to the region when its centre does + tol = 1.0e-9 + return findall(c -> (lo - tol) <= c <= (hi + tol), coords) + end +end + +""" + topo_cache_dir() + +Where downloaded tiles are kept, so that repeated imports of the same region do not download +them again. Uses a scratch space, which is removed when the package is. +""" +topo_cache_dir() = @get_scratch!("topo_tiles") + +""" + read_topo_netcdf(file) + +Read a global `.grd` grid, returning `(lon, lat, z)`. Filled in by the `NCDatasets` +extension; the tiled resolutions do not go through here. +""" +function read_topo_netcdf(file::AbstractString) + 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 + +""" + download_tile(url, dest) + +Fetch one file, trying each mirror in turn. Returns `true` if it was downloaded, `false` if +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) + isfile(dest) && filesize(dest) > 0 && return true # cached + for server in GMT_SERVERS + try + Downloads.download("$server/$relative_url", dest) + return true + catch err + # a 404 means there is no such tile; anything else is worth trying elsewhere + if err isa Downloads.RequestError && err.response.status == 404 + return false + end + end + end + return false +end + +""" + decode_jp2(file, n) + +Decode a JPEG2000 tile into an `n x n` matrix of `Int16` elevations, using the decoder that +comes with `OpenJpeg_jll`. The file stores its rows from north to south. +""" +function decode_jp2(file::AbstractString, n::Integer) + raw = file * ".raw" + if !(isfile(raw) && filesize(raw) == n * n * 2) + run(pipeline(`$(opj_decompress()) -i $file -o $raw`, stdout = devnull, stderr = devnull)) + end + A = Array{Int16}(undef, n, n) + read!(raw, A) + return permutedims(A) # to (row, column) = (lat, lon) +end + +""" + grid_from_tiles(limits, dataset, res, reg) + +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) + 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 + # resampling, so follow the set and only honour `reg` where there is a choice. + tile_reg = native_registration(dataset, res, reg) + deg = tile_size_deg(dataset, res, tile_reg) + n = tile_size_px(res, deg, tile_reg) + px = deg / (tile_reg == "g" ? n - 1 : n) # degrees per pixel + + # The tiles are named after their south-west corner. Longitude tiles start at the + # antimeridian and latitude tiles at the south pole, so for a tile larger than 30 degrees + # the latitude anchors are -90, -90+deg, ... rather than multiples of `deg`. + snap(v, origin) = origin + Int(floor((v - origin) / deg)) * deg + lat0 = snap(latmin, -90) + lat1 = snap(latmax - 1.0e-9, -90) + lon0 = snap(lonmin, -180) + lon1 = snap(lonmax - 1.0e-9, -180) + + lats = lat0:deg:lat1 + lons = lon0:deg:lon1 + + # Gridline-registered tiles repeat the edge they share with the next tile along, so a + # 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")) + dir = topo_cache_dir() + sub = tile_directory(dataset, res, tile_reg) + + missing_tiles = Tuple{Int, Int}[] + for (j, la) in enumerate(lats), (i, lo) in enumerate(lons) + name = tile_filename(la, lo, dataset, res, tile_reg) + file = joinpath(dir, name) + if !download_tile("$dataset/$sub/$name", file) + push!(missing_tiles, (la, lo)) # filled from the coarser set below + 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, :] + end + + # the coordinates of each sample: cell centres for pixel registration, cell corners + # (so starting exactly on the tile edge) for gridline registration + off = tile_reg == "g" ? 0.0 : 0.5 + lon = [lon0 + (k - 1 + off) * px for k in 1:size(Z, 2)] + lat = [lat0 + (k - 1 + off) * px for k in 1:size(Z, 1)] + ix = subset_indices(lon, lonmin, lonmax, px, tile_reg) + iy = subset_indices(lat, latmin, latmax, px, tile_reg) + + isempty(ix) && error("no data between longitudes $lonmin and $lonmax") + isempty(iy) && error("no data between latitudes $latmin and $latmax") + + # the tiles hold integers in units the server declares, not metres + scale, offset = dataset_scale(dataset, res, tile_reg) + out = Float64.(Z[iy, ix]) .* scale .+ offset + + # where a tile was missing there is no data at all, only the zeros we started with. + # Fill those from the coarser set that the server names for the purpose, so the sea has + # its depth rather than being flat. + filler = filler_dataset(dataset, res, tile_reg) + if !isnothing(filler) + # SRTM records the sea as a flat zero rather than leaving it out, and has no tile at + # all where there is no land, so both cases have to be filled -- every zero, not only + # the tiles that were missing. GMT does the same, which is why its 03s has bathymetry. + fill_gaps!(out, lon[ix], lat[iy], filler) + end + + return lon[ix], lat[iy], out +end + +""" + fill_gaps!(Z, lon, lat, filler) + +Replace the cells of `Z` that carry no data -- the flat zeros that SRTM writes for the sea, +and whatever a missing tile left behind -- with values interpolated from the coarser `filler` +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) + 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]) + end + return Z +end + +""" + bilinear(xs, ys, A, x, y) + +Sample `A` at (`x`, `y`) by interpolating between the four surrounding points, which is what +GMT does when it resamples a coarser set into a finer grid. Falls back to the nearest edge +value outside the grid. +""" +function bilinear(xs::AbstractVector, ys::AbstractVector, A::AbstractMatrix, x::Real, y::Real) + i = searchsortedlast(xs, x) + j = searchsortedlast(ys, y) + i = clamp(i, 1, length(xs) - 1) + j = clamp(j, 1, length(ys) - 1) + + tx = (x - xs[i]) / (xs[i + 1] - xs[i]) + ty = (y - ys[j]) / (ys[j + 1] - ys[j]) + tx = clamp(tx, 0.0, 1.0) + 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] +end + +""" + grid_from_single_file(limits, dataset, res, reg) + +The coarse resolutions are published as one NetCDF file for the whole globe, which is small +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) + lonmin, lonmax, latmin, latmax = limits + name = "$(dataset)_$(res)_$(reg).grd" + file = joinpath(topo_cache_dir(), name) + download_tile("$dataset/$name", file) || + error("could not download $name from the GMT data server") + + lon, lat, Z = read_topo_netcdf(file) + px = length(lon) > 1 ? abs(lon[2] - lon[1]) : 1.0 + ix = subset_indices(lon, lonmin, lonmax, px, reg) + iy = subset_indices(lat, latmin, latmax, px, reg) + + isempty(ix) && error("no data between longitudes $lonmin and $lonmax") + isempty(iy) && error("no data between latitudes $latmin and $latmax") + + return lon[ix], lat[iy], Float64.(Z[iy, ix]) +end + +""" + import_topo(limits; dataset="earth_relief", res="01m", reg="p") + +Download topography (and bathymetry) for the region `limits = [lonmin, lonmax, latmin, latmax]` +from the GMT data server, and return it as `GeoData` whose `Topography` field is in km. + +This does not require `GMT`: the grids are fetched over HTTP and decoded here. + +# Datasets +- `"earth_relief"` (the default) blends land topography with ocean bathymetry +- `"earth_gebco"` is GEBCO bathymetry +- `"earth_synbath"` is GEBCO with synthetic bathymetry where the sea floor is unsurveyed + +# Resolutions +`"01d"`, `"30m"`, `"20m"`, `"15m"`, `"10m"`, `"06m"` come as one file for the globe; +`"05m"` through `"01m"`, `"30s"`, `"15s"` and `"03s"` are tiled. + +Note that `"03s"` is SRTM, which covers **land only** -- the sea is flat zero there, and the +finest resolution that carries real bathymetry is `"15s"`. + +# Example +```julia +julia> Topo = import_topo([-18.7, -17.1, 28.0, 29.2], res="03s") +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) + + # accept the `file="@earth_relief_01m"` spelling that this function used to take + if !isnothing(file) + dataset, res, reg = parse_topo_file(file) + end + + # longitudes west of Greenwich may be given either way round + limits = collect(float.(limits)) + if limits[1] > limits[2] + limits[1], limits[2] = limits[2], limits[1] + end + if limits[3] > limits[4] + limits[3], limits[4] = limits[4], limits[3] + end + + lon, lat, Z = isnothing(tile_size_deg(dataset, res, native_registration(dataset, res, reg))) ? + grid_from_single_file(limits, dataset, res, reg) : + grid_from_tiles(limits, dataset, res, reg) + + Lon, Lat, Depth = lonlatdepth_grid(lon, lat, 0) + @views Depth[:, :, 1] = 1.0e-3 * Z' # the field is in km + return GeoData(Lon, Lat, Depth, (Topography = Depth * km,)) +end + +import_topo(; lat = [37, 49], lon = [4, 20], kwargs...) = + import_topo([lon[1], lon[2], lat[1], lat[2]]; kwargs...) + +""" + parse_topo_file(file) + +Take apart the `"@earth_relief_01m"` spelling into the dataset, resolution and registration +that name it on the server. +""" +function parse_topo_file(file::AbstractString) + s = lstrip(file, '@') + 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]) +end diff --git a/test/runtests.jl b/test/runtests.jl index 5474eac48..6e9a93ea9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using ParallelTestRunner # concurrently (which times out regularly in CI). See topo_prefetch.jl. include("topo_prefetch.jl") +include("test_import_topo.jl") testsuite = find_tests(@__DIR__) # Not a test file, only a helper delete!(testsuite, "topo_prefetch") diff --git a/test/test_import_topo.jl b/test/test_import_topo.jl new file mode 100644 index 000000000..928ba9f21 --- /dev/null +++ b/test/test_import_topo.jl @@ -0,0 +1,88 @@ +# Tests `import_topo`, which downloads from the GMT data server without needing GMT. +# +# These hit the network, so they are written to be cheap: small regions, and the tiles are +# cached between runs. +using Test +using GeophysicalModelGenerator + +@testset "import_topo" begin + + # the tile addressing, which is where the awkward cases live + @testset "tile naming" begin + G = GeophysicalModelGenerator + # named after the south-west corner, zero padded + @test G.tile_filename(28, -18, "earth_relief", "03s", "g") == "N28W018.SRTMGL3.jp2" + @test G.tile_filename(0, -30, "earth_relief", "01m", "p") == "N00W030.earth_relief_01m_p.jp2" + @test G.tile_filename(-30, 60, "earth_relief", "15s", "p") == "S30E060.earth_relief_15s_p.jp2" + + # pixel registration samples cell centres, gridline the corners, so one more of them + @test G.tile_size_px("01m", 30, "p") == 1800 + @test G.tile_size_px("01m", 30, "g") == 1801 + @test G.tile_size_px("03s", 1, "g") == 1201 + + # the sets that exist in only one registration + @test G.native_registration("earth_relief", "03s", "p") == "g" + @test G.native_registration("earth_relief", "15s", "g") == "p" + end + + # `"@earth_relief_01m"`, as this function used to be called + @testset "file= spelling" begin + G = GeophysicalModelGenerator + @test G.parse_topo_file("@earth_relief_01m") == ("earth_relief", "01m", "g") + @test G.parse_topo_file("earth_relief_15s_p") == ("earth_relief", "15s", "p") + @test G.parse_topo_file("@earth_gebco_30s_g") == ("earth_gebco", "30s", "g") + end + + # La Palma: small, has both land and deep water, and a summit to check against + limits = [-18.7, -17.1, 28.0, 29.2] + + @testset "downloads and georeferencing" begin + Topo = import_topo(limits, res = "01m") + @test Topo isa GeoData + + z = ustrip.(Topo.fields.Topography) # km + @test size(z, 3) == 1 + @test -10 < minimum(z) < 0 # deep water off the island + @test 2 < maximum(z) < 3 # Roque de los Muchachos, 2.4 km + + # the region we asked for, and no more + @test minimum(Topo.lon.val) ≥ limits[1] - 0.1 + @test maximum(Topo.lon.val) ≤ limits[2] + 0.1 + @test minimum(Topo.lat.val) ≥ limits[3] - 0.1 + @test maximum(Topo.lat.val) ≤ limits[4] + 0.1 + end + + @testset "summit is in the right place" begin + # 3 arcseconds resolves the summit; it should land on the real one to within a + # pixel, which is what catches a transposed or mis-assembled grid + Topo = import_topo(limits, res = "03s") + z = ustrip.(Topo.fields.Topography) + i = argmax(z) + @test isapprox(Topo.lon.val[i], -17.885, atol = 0.002) + @test isapprox(Topo.lat.val[i], 28.754, atol = 0.002) + @test isapprox(maximum(z) * 1000, 2414, atol = 2) # metres + end + + @testset "bathymetry" begin + # 3 arcseconds is SRTM, which is land only -- the sea comes from the coarser set the + # server names for it, so it must not be flat + Topo = import_topo(limits, res = "03s") + z = ustrip.(Topo.fields.Topography) + @test minimum(z) < -3 # km, real depth + @test count(<(-1), z) > 0.5 * length(z) # mostly ocean here + end + + @testset "other ways of asking" begin + a = import_topo(limits, res = "30s") + b = import_topo(lon = [limits[1], limits[2]], lat = [limits[3], limits[4]], res = "30s") + @test ustrip.(a.fields.Topography) == ustrip.(b.fields.Topography) + + c = import_topo(limits, file = "@earth_relief_30s") + @test size(c.fields.Topography) == size(a.fields.Topography) + end + + @testset "other datasets" begin + gebco = import_topo(limits, dataset = "earth_gebco", res = "15s") + @test minimum(ustrip.(gebco.fields.Topography)) < -3 # km + end +end From 411ddd1fcec0d99baf2667422c806bf25005b1d6 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 29 Sep 2026 10:01:42 +0200 Subject: [PATCH 2/3] Accept the call forms the existing callers use CI caught two: `maxattempts`, which the prefetch script passes and which is now honoured as the number of times each mirror is tried, and the `.grd`/`.nc` suffix that the tutorials append to the file name, which says which format GMT downloaded into rather than which set is wanted. Also a spelling the linter objects to. Co-Authored-By: Claude Opus 5 (1M context) --- src/import_topo.jl | 41 ++++++++++++++++++++++++---------------- test/test_import_topo.jl | 2 +- 2 files changed, 26 insertions(+), 17 deletions(-) diff --git a/src/import_topo.jl b/src/import_topo.jl index 01856adf8..37b262a12 100644 --- a/src/import_topo.jl +++ b/src/import_topo.jl @@ -231,16 +231,18 @@ 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) +function download_tile(relative_url::AbstractString, dest::AbstractString; maxattempts::Integer = 5) isfile(dest) && filesize(dest) > 0 && return true # cached - for server in GMT_SERVERS - try - Downloads.download("$server/$relative_url", dest) - return true - catch err - # a 404 means there is no such tile; anything else is worth trying elsewhere - if err isa Downloads.RequestError && err.response.status == 404 - return false + for attempt in 1:max(1, maxattempts) + for server in GMT_SERVERS + try + Downloads.download("$server/$relative_url", dest) + return true + catch err + # a 404 means there is no such tile, and no mirror will have it either + if err isa Downloads.RequestError && err.response.status == 404 + return false + end end end end @@ -269,7 +271,8 @@ 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) +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 @@ -304,7 +307,7 @@ function grid_from_tiles(limits, dataset::AbstractString, res::AbstractString, r for (j, la) in enumerate(lats), (i, lo) in enumerate(lons) name = tile_filename(la, lo, dataset, res, tile_reg) file = joinpath(dir, name) - if !download_tile("$dataset/$sub/$name", file) + if !download_tile("$dataset/$sub/$name", file; maxattempts = maxattempts) push!(missing_tiles, (la, lo)) # filled from the coarser set below continue end @@ -390,11 +393,12 @@ 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) +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) - download_tile("$dataset/$name", file) || + download_tile("$dataset/$name", file; maxattempts = maxattempts) || error("could not download $name from the GMT data server") lon, lat, Z = read_topo_netcdf(file) @@ -438,7 +442,8 @@ Tiles are cached, so importing the same region again does not download them a se """ function import_topo(limits; dataset::AbstractString = "earth_relief", res::AbstractString = "01m", reg::AbstractString = "g", - file::Union{Nothing, AbstractString} = nothing) + file::Union{Nothing, AbstractString} = nothing, + maxattempts::Integer = 5) # accept the `file="@earth_relief_01m"` spelling that this function used to take if !isnothing(file) @@ -455,8 +460,8 @@ function import_topo(limits; dataset::AbstractString = "earth_relief", end lon, lat, Z = isnothing(tile_size_deg(dataset, res, native_registration(dataset, res, reg))) ? - grid_from_single_file(limits, dataset, res, reg) : - grid_from_tiles(limits, dataset, res, reg) + grid_from_single_file(limits, dataset, res, reg; maxattempts = maxattempts) : + grid_from_tiles(limits, dataset, res, reg; maxattempts = maxattempts) Lon, Lat, Depth = lonlatdepth_grid(lon, lat, 0) @views Depth[:, :, 1] = 1.0e-3 * Z' # the field is in km @@ -474,6 +479,10 @@ that name it on the server. """ function parse_topo_file(file::AbstractString) s = lstrip(file, '@') + # the tutorials write "@earth_relief_01m.grd"; the suffix names the format GMT used to + # download into and says nothing about which set is wanted + endswith(s, ".grd") && (s = s[1:(end - 4)]) + endswith(s, ".nc") && (s = s[1:(end - 3)]) 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 diff --git a/test/test_import_topo.jl b/test/test_import_topo.jl index 928ba9f21..5e7775d5a 100644 --- a/test/test_import_topo.jl +++ b/test/test_import_topo.jl @@ -54,7 +54,7 @@ using GeophysicalModelGenerator @testset "summit is in the right place" begin # 3 arcseconds resolves the summit; it should land on the real one to within a - # pixel, which is what catches a transposed or mis-assembled grid + # pixel, which is what catches a transposed or wrongly assembled grid Topo = import_topo(limits, res = "03s") z = ustrip.(Topo.fields.Topography) i = argmax(z) From a38de956761b4f9769a07ac50c548ccdd4e84ad2 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 29 Sep 2026 11:34:37 +0200 Subject: [PATCH 3/3] Do not let a throttled download hang the run `Downloads.download` waits forever by default, and the retry loop tried each of four mirrors up to five times, so a server that accepts the connection and then stalls blocked indefinitely rather than failing over. That is what the data server does when several CI jobs pull tiles at once: on one run three macOS jobs finished in 4-8 minutes while a fourth on the same platform was still going after an hour. Give each request a timeout, back off between rounds of mirrors, and download through a .part file so a failed transfer cannot leave a truncated tile behind. The prefetch is still worth having -- it serialises the first downloads, which is what keeps the workers off the server -- but it was warming ~/.gmt, which nothing reads any more. Point it, and the CI cache, at the scratch space that `import_topo` actually uses. ~/.gmt stays cached for `import_GeoTIFF`, which does still go through GMT. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/CI.yml | 14 +++++++++++--- src/import_topo.jl | 25 ++++++++++++++++++++++--- test/topo_prefetch.jl | 29 +++++++++++++++-------------- 3 files changed, 48 insertions(+), 20 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0d51c0971..0b4be2ba9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -70,21 +70,29 @@ jobs: with: include-matrix: false - uses: julia-actions/julia-buildpkg@latest - # Cache the GMT topography tiles that the tests download (see + # Cache the topography tiles that the tests download (see # test/topo_prefetch.jl), so PR runs do not depend on the GMT data # server. The scheduled run skips the cache on purpose, so that the # download itself is still tested regularly. + # + # `import_topo` keeps the tiles in its own scratch space rather than in + # ~/.gmt, which is where GMT used to put them; the UUID is the package's. + # ~/.gmt is still cached for import_GeoTIFF, which does go through GMT. - uses: actions/cache/restore@v4 id: gmt-cache if: github.event_name != 'schedule' with: - path: ~/.gmt + path: | + ~/.julia/scratchspaces/3700c31b-fa53-48a6-808a-ef22d5a84742/topo_tiles + ~/.gmt key: gmt-topo-${{ runner.os }}-${{ runner.arch }}-${{ hashFiles('test/topo_prefetch.jl') }} - uses: julia-actions/julia-runtest@latest - uses: actions/cache/save@v4 if: always() && github.event_name != 'schedule' && steps.gmt-cache.outputs.cache-hit != 'true' with: - path: ~/.gmt + path: | + ~/.julia/scratchspaces/3700c31b-fa53-48a6-808a-ef22d5a84742/topo_tiles + ~/.gmt key: gmt-topo-${{ runner.os }}-${{ runner.arch }}-${{ hashFiles('test/topo_prefetch.jl') }} - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v7 diff --git a/src/import_topo.jl b/src/import_topo.jl index 37b262a12..b95346461 100644 --- a/src/import_topo.jl +++ b/src/import_topo.jl @@ -9,6 +9,13 @@ using Downloads, OpenJpeg_jll, Scratch +""" +How long to wait for one tile before giving up on that mirror and trying the next. The +transfers are a few megabytes, so a minute is generous; without a limit a throttled or +stalled connection blocks indefinitely. +""" +const TILE_TIMEOUT = 60.0 + """ The GMT data server mirrors. The first that answers is used; the others are tried in turn if a download fails, so one server being offline does not fail the import. @@ -70,7 +77,7 @@ function server_index_file() for server in GMT_SERVERS root = replace(server, "/server/earth" => "") try - Downloads.download("$root/gmt_data_server.txt", index) + Downloads.download("$root/gmt_data_server.txt", index; timeout = TILE_TIMEOUT) return index catch end @@ -231,20 +238,32 @@ 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) +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 for attempt in 1:max(1, maxattempts) for server in GMT_SERVERS try - Downloads.download("$server/$relative_url", dest) + # `Downloads.download` waits forever by default. The data server throttles + # when several jobs pull tiles at once -- which is exactly what CI does -- + # and without a timeout a throttled connection hangs the whole run rather + # than failing over to the next mirror. + Downloads.download("$server/$relative_url", tmp; timeout = timeout) + mv(tmp, dest; force = true) return true catch err + rm(tmp, force = true) # a 404 means there is no such tile, and no mirror will have it either if err isa Downloads.RequestError && err.response.status == 404 return false end end end + # back off a little before going round the mirrors again, so a server that is busy + # is given a chance rather than hammered + attempt < maxattempts && sleep(min(2.0^attempt, 10.0)) end return false end diff --git a/test/topo_prefetch.jl b/test/topo_prefetch.jl index 280dffb29..c5a682317 100644 --- a/test/topo_prefetch.jl +++ b/test/topo_prefetch.jl @@ -1,19 +1,20 @@ -# Warm the GMT topography cache before the parallel test workers start. +# Warm the topography cache before the parallel test workers start. # # Several tests (test_GMT, test_WaterFlow, the LaPalma tutorial) download -# topography tiles through GMT. Under ParallelTestRunner they run concurrently, -# so they all hit the GMT data server at the same time, and the transfers then -# regularly time out in CI ("Libcurl Error: Timeout was reached"). +# topography tiles. Under ParallelTestRunner they run concurrently, so they all +# hit the GMT data server at the same time, and the server throttles: the +# transfers then crawl or time out in CI. # # Downloading the tiles here, one at a time and with generous retries, means the -# workers later find them in GMT's cache (~/.gmt/server) and never touch the -# network. The download itself is still exercised: it happens right here, and -# on CI the cache is bypassed for the scheduled run so it is tested regularly. +# workers later find them in the cache and never touch the network. The download +# itself is still exercised: it happens right here, and on CI the cache is +# bypassed for the scheduled run so it is tested regularly. # -# NOTE: .github/workflows/CI.yml caches ~/.gmt keyed on the hash of this file, -# so changing the list below invalidates that cache automatically. +# NOTE: .github/workflows/CI.yml caches the scratch space that `import_topo` +# keeps its tiles in, keyed on the hash of this file, so changing the list below +# invalidates that cache automatically. -using GeophysicalModelGenerator, GMT +using GeophysicalModelGenerator # (region, file) pairs, matching the calls made in the tests/tutorials const TOPO_PREFETCH = [ @@ -23,8 +24,8 @@ const TOPO_PREFETCH = [ ] function prefetch_topography() - gmt_cache = joinpath(get(ENV, "GMT_USERDIR", joinpath(homedir(), ".gmt")), "server") - n_before = isdir(gmt_cache) ? count(f -> endswith(f, ".nc"), collect(Iterators.flatten(map(t -> t[3], walkdir(gmt_cache))))) : 0 + cache = GeophysicalModelGenerator.topo_cache_dir() + n_before = isdir(cache) ? count(f -> endswith(f, ".jp2") || endswith(f, ".grd"), readdir(cache)) : 0 for (limits, file) in TOPO_PREFETCH t = @elapsed try @@ -38,8 +39,8 @@ function prefetch_topography() @info "Prefetched topography $file for $limits in $(round(t, digits = 1)) s" end - n_after = isdir(gmt_cache) ? count(f -> endswith(f, ".nc"), collect(Iterators.flatten(map(t -> t[3], walkdir(gmt_cache))))) : 0 - @info "GMT topography cache: $n_after tiles in $gmt_cache ($(n_after - n_before) newly downloaded)" + n_after = isdir(cache) ? count(f -> endswith(f, ".jp2") || endswith(f, ".grd"), readdir(cache)) : 0 + @info "Topography cache: $n_after tiles in $cache ($(n_after - n_before) newly downloaded)" return nothing end