Skip to content

Commit

Permalink
Merge pull request #122 from albert-de-montserrat/adm/CI
Browse files Browse the repository at this point in the history
Fix the download of topographyic data with GMT
  • Loading branch information
albert-de-montserrat authored Jun 19, 2024
2 parents 4b3886a + 285d400 commit 432ac9c
Show file tree
Hide file tree
Showing 7 changed files with 39 additions and 40 deletions.
2 changes: 1 addition & 1 deletion docs/src/man/Tutorial_AlpineData.md
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ using GeophysicalModelGenerator, GMT
When loading both packages, several `GMT` routines within `GMG` will be loaded. One of these routines is the function `import_topo`, where one simply has to provide the region for which to download the topographic data and the data source.

```julia
Topo = import_topo([4,20,37,50], file="@earth_relief_01m.grd")
Topo = import_topo([4,20,37,50], file="@earth_relief_01m")
```

The data is available in different resolutions; see [here](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html) for an overview. Generally, it is advisable to not use the largest resolution if you have a large area, as the files become very large.
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/Tutorial_Jura.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ using GeophysicalModelGenerator, GMT
Download the topography with:

```julia
Topo = import_topo(lat=[45.5,47.7], lon=[5, 8.1], file="@earth_relief_03s.grd")
Topo = import_topo(lat=[45.5,47.7], lon=[5, 8.1], file="@earth_relief_03s")
```

Next, we drape the geological map on top of the geological map.
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/Tutorial_LaPalma.md
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ using GeophysicalModelGenerator, GMT, DelimitedFiles
We will use GMT to download the topography with:

```julia
Topo = import_topo(lon = [-18.2, -17.5], lat=[28.4, 29.0], file="@earth_relief_15s.grd")
Topo = import_topo(lon = [-18.2, -17.5], lat=[28.4, 29.0], file="@earth_relief_15s")
```

Next, lets load the seismicity. The earthquake data is available on [https://www.ign.es/web/ign/portal/vlc-catalogo](https://www.ign.es/web/ign/portal/vlc-catalogo).
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/tutorial_GMT_Topography.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ The nice thing about GMT is that it automatically downloads data for you for a c

```julia
julia> using GeophysicalModelGenerator, GMT
julia> Topo = import_topo([4,20,37,49], file="@earth_relief_01m.grd")
julia> Topo = import_topo([4,20,37,49], file="@earth_relief_01m")
GeoData
size : (960, 720, 1)
lon ϵ [ 4.0 : 19.983333333333334]
Expand Down
4 changes: 2 additions & 2 deletions docs/src/man/tutorial_GMT_Topography_GeologicalMap.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ In many cases, we want to add topographic data as well a information about tecto

#### 1. Download topographic data and tectonic maps of the Alpine region
The ETOPO1 data file used in this example can be downloaded here:
[https://ngdc.noaa.gov/mgg/global/global.html](https://ngdc.noaa.gov/mgg/global/global.html). For this example, we downloaded `ETOPO1_Ice_g_gmt4.grd` and stored it directly in the folder where we will be working. For the geological map, we download the data from the [SPP 4DMB repository](http://www.spp-mountainbuilding.de/data/Maps.zip) and extract the zip file (to the current folder). In this data set, a `gmt` file with the data for different tectonic units is given in `./tectonic_maps_4dmb_2020_09_17/GMT_example/alcapadi_polygons.gmt`.
[https://ngdc.noaa.gov/mgg/global/global.html](https://ngdc.noaa.gov/mgg/global/global.html). For this example, we downloaded `ETOPO1_Ice_g_gmt4` and stored it directly in the folder where we will be working. For the geological map, we download the data from the [SPP 4DMB repository](http://www.spp-mountainbuilding.de/data/Maps.zip) and extract the zip file (to the current folder). In this data set, a `gmt` file with the data for different tectonic units is given in `./tectonic_maps_4dmb_2020_09_17/GMT_example/alcapadi_polygons.gmt`.

#### 2. Create a tectonic map with orthogonal projection
To create a png with an orthogonal map projection (which we need for the png import), we do the following in julia:
Expand All @@ -28,7 +28,7 @@ julia> using GMT, NearestNeighbors, GeoParams, GeophysicalModelGenerator
```
First, define the filenames of the files you want to import:
```julia
julia> filename_topo = "./ETOPO1/ETOPO1_Ice_g_gmt4.grd"
julia> filename_topo = "./ETOPO1/ETOPO1_Ice_g_gmt4"
julia> filename_geo = "./tectonicmap_SPP.png"
```
Next, define the region that you want to visualize (note that we use the same coordinates here as we used previously for the generation of the geological map):
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/tutorial_GPS.md
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,7 @@ using GMT, Interpolations
We use the `import_topo` function to read the topography from a file:

```julia
Elevation = import_topo([3,17,42,50], file="@earth_relief_01m.grd");
Elevation = import_topo([3,17,42,50], file="@earth_relief_01m");
nothing #hide
```

Expand Down
65 changes: 32 additions & 33 deletions ext/GMT_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ println("Loading GMT routines within GMG")


"""
Topo = import_topo(limits; file::String="@earth_relief_01m.grd", maxattempts=5)
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.
Expand Down Expand Up @@ -65,45 +65,44 @@ julia> write_paraview(Topo,"Topo_Alps")
"Topo_Alps.vts"
```
"""
function import_topo(limits; file::String="@earth_relief_01m.grd", maxattempts=5)
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);
# Correct if negative values are given (longitude coordinates that are west)
ind = limits[1:2] .< 0

if (limits[1] < 0) && (limits[2] < 0)
limits[ind] .= 360 .+ limits[ind];
limits[1:2] = sort(limits[1:2])
end
if (limits[1] < 0) && (limits[2] < 0)
limits[ind] .= 360 .+ limits[ind]
limits[1:2] = sort(limits[1:2])
end

# Download topo file - add a few attempts to do so
G = [];
attempt = 0
while attempt<maxattempts
try
G = gmtread(file, limits=limits, grid=true);
break
catch
@warn "Failed downloading GMT topography on attempt $attempt/$maxattempts"
sleep(5) # wait a few sec
end
attempt += 1
end
if isempty(G)
error("Could not download GMT topography data")
# Download topo file - add a few attempts to do so
local G
attempt = 0
while attempt < maxattempts
try
G = gmtread(file, limits=limits, grid=true);
break
catch
@warn "Failed downloading GMT topography on attempt $attempt/$maxattempts"
sleep(5) # wait a few sec
end
attempt += 1
end
if isempty(G)
error("Could not download GMT topography data")
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);
Depth[:,:,1] = 1e-3*G.z';
Topo = GeoData(Lon, Lat, Depth, (Topography=Depth*km,))
# 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] = 1e-3*G.z';
Topo = GeoData(Lon, Lat, Depth, (Topography=Depth*km,))

return Topo

return Topo
end

"""
import_topo(; lat::Vector{2}, lon::Vector{2}, file::String="@earth_relief_01m.grd", maxattempts=5)
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
Expand All @@ -114,11 +113,11 @@ 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.grd")
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.grd", maxattempts=5) = import_topo([lon[1],lon[2], lat[1], lat[2]], file=file, maxattempts=maxattempts)
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)


"""
Expand Down

0 comments on commit 432ac9c

Please sign in to comment.