mthielma · Sep 29, 2026 · Sep 29, 2026 · Sep 29, 2026 · Sep 29, 2026
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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