Skip to content

Commit 61ca54a

Browse files
fix GMT import_topo
1 parent e0d90c4 commit 61ca54a

File tree

1 file changed

+32
-33
lines changed

1 file changed

+32
-33
lines changed

ext/GMT_utils.jl

Lines changed: 32 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ println("Loading GMT routines within GMG")
1818

1919

2020
"""
21-
Topo = import_topo(limits; file::String="@earth_relief_01m.grd", maxattempts=5)
21+
Topo = import_topo(limits; file::String="@earth_relief_01m", maxattempts=5)
2222
2323
Uses `GMT` to download the topography of a certain region, specified with limits=[lon_min, lon_max, lat_min, lat_max].
2424
Sometimes download fails because of the internet connection. We do `maxattempts` to download it.
@@ -65,45 +65,44 @@ julia> write_paraview(Topo,"Topo_Alps")
6565
"Topo_Alps.vts"
6666
```
6767
"""
68-
function import_topo(limits; file::String="@earth_relief_01m.grd", maxattempts=5)
68+
function import_topo(limits; file::String="@earth_relief_01m", maxattempts=5)
6969

70-
# Correct if negative values are given (longitude coordinates that are west)
71-
ind = findall(limits[1:2] .< 0);
70+
# Correct if negative values are given (longitude coordinates that are west)
71+
ind = limits[1:2] .< 0
7272

73-
if (limits[1] < 0) && (limits[2] < 0)
74-
limits[ind] .= 360 .+ limits[ind];
75-
limits[1:2] = sort(limits[1:2])
76-
end
73+
if (limits[1] < 0) && (limits[2] < 0)
74+
limits[ind] .= 360 .+ limits[ind]
75+
limits[1:2] = sort(limits[1:2])
76+
end
7777

78-
# Download topo file - add a few attempts to do so
79-
G = [];
80-
attempt = 0
81-
while attempt<maxattempts
82-
try
83-
G = gmtread(file, limits=limits, grid=true);
84-
break
85-
catch
86-
@warn "Failed downloading GMT topography on attempt $attempt/$maxattempts"
87-
sleep(5) # wait a few sec
88-
end
89-
attempt += 1
90-
end
91-
if isempty(G)
92-
error("Could not download GMT topography data")
78+
# Download topo file - add a few attempts to do so
79+
local G
80+
attempt = 0
81+
while attempt < maxattempts
82+
try
83+
G = gmtread(file, limits=limits, grid=true);
84+
break
85+
catch
86+
@warn "Failed downloading GMT topography on attempt $attempt/$maxattempts"
87+
sleep(5) # wait a few sec
9388
end
89+
attempt += 1
90+
end
91+
if isempty(G)
92+
error("Could not download GMT topography data")
93+
end
9494

95-
# Transfer to GeoData
96-
nx,ny = size(G.z,2), size(G.z,1)
97-
Lon,Lat,Depth = lonlatdepth_grid(G.x[1:nx],G.y[1:ny],0);
98-
Depth[:,:,1] = 1e-3*G.z';
99-
Topo = GeoData(Lon, Lat, Depth, (Topography=Depth*km,))
95+
# Transfer to GeoData
96+
nx, ny = size(G.z,2), size(G.z,1)
97+
Lon,Lat,Depth = lonlatdepth_grid(G.x[1:nx],G.y[1:ny],0);
98+
@views Depth[:,:,1] = 1e-3*G.z';
99+
Topo = GeoData(Lon, Lat, Depth, (Topography=Depth*km,))
100100

101-
return Topo
102-
101+
return Topo
103102
end
104103

105104
"""
106-
import_topo(; lat::Vector{2}, lon::Vector{2}, file::String="@earth_relief_01m.grd", maxattempts=5)
105+
import_topo(; lat::Vector{2}, lon::Vector{2}, file::String="@earth_relief_01m", maxattempts=5)
107106
108107
Imports topography (using GMT), by specifying keywords for latitude and longitude ranges
109108
@@ -114,11 +113,11 @@ julia> Topo = import_topo(lat=[30,40], lon=[30, 50] )
114113
```
115114
The values can also be defined as tuples:
116115
```julia
117-
julia> Topo = import_topo(lon=(-50, -40), lat=(-10,-5), file="@earth_relief_30s.grd")
116+
julia> Topo = import_topo(lon=(-50, -40), lat=(-10,-5), file="@earth_relief_30s")
118117
```
119118
120119
"""
121-
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)
120+
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)
122121

123122

124123
"""

0 commit comments

Comments
 (0)