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/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..b95346461 --- /dev/null +++ b/src/import_topo.jl @@ -0,0 +1,510 @@ +# 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 + +""" +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. +""" +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; timeout = TILE_TIMEOUT) + 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; + 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` 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 + +""" + 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; + 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 + # 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; maxattempts = maxattempts) + 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; + maxattempts::Integer = 5) + lonmin, lonmax, latmin, latmax = limits + name = "$(dataset)_$(res)_$(reg).grd" + file = joinpath(topo_cache_dir(), name) + download_tile("$dataset/$name", file; maxattempts = maxattempts) || + 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, + maxattempts::Integer = 5) + + # 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; 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 + 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, '@') + # 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 + # 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..5e7775d5a --- /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 wrongly 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 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