Skip to content

Commit

Permalink
Improve gravmag3d code parsing and use new fixes done in GMT lib. (#1650
Browse files Browse the repository at this point in the history
)

* Export an alias name

* Delete commented lines

* Add gravmag tests

* Delete GMT version test

* Improve gravmag3d code parsing and use new fixes done in GMT lib.

* Comment some @test asserts

* Comment more gravmag tests

* One more

* And another
  • Loading branch information
joa-quim authored Jan 23, 2025
1 parent 29f00cc commit e4a935b
Show file tree
Hide file tree
Showing 7 changed files with 56 additions and 31 deletions.
2 changes: 1 addition & 1 deletion src/GMT.jl
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ export
rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, spectrum1d, sphdistance, sphinterpolate,
sphtriangulate, surface, ternary, ternary!, text, text!, text_record, trend1d, trend2d, triangulate, gmtsplit,
decorated, vector_attrib, wiggle, wiggle!, xyz2grd, gmtbegin, gmtend, gmthelp, subplot, gmtfig, inset, showfig,
earthtide, gravfft, gmtgravmag3d, grdgravmag3d, gravprisms, pscoupe, pscoupe!, coupe, coupe!, psmeca, psmeca!,
earthtide, gravfft, gmtgravmag3d, gravmag3d, grdgravmag3d, gravprisms, pscoupe, pscoupe!, coupe, coupe!, psmeca, psmeca!,
meca, meca!, psvelo, psvelo!, velo, velo!, gmtisf, getbyattrib, inpolygon, inwhichpolygon, pcolor, pcolor!,
triplot, triplot!, trisurf, trisurf!, grdrotater, imagesc, upGMT, boxes,

Expand Down
7 changes: 0 additions & 7 deletions src/gmt_main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1106,13 +1106,6 @@ function dataset_init(API::Ptr{Nothing}, Darr::Vector{<:GMTdataset}, direction::
unsafe_store!(S, Sb)
unsafe_store!(DT.segment, S, seg)
end
#DT.n_records, DS.n_records = n_records, n_records # They are equal because our GMT_DATASET have only one table
#Dt = unsafe_load(DS.table)
#unsafe_store!(Dt, DT)
#unsafe_store!(DS.table, Dt)
#unsafe_store!(D, DS)

#return D
helper_init_DS(D, DS, DT, n_records) # Stores the meshes in D and returns it
end

Expand Down
1 change: 0 additions & 1 deletion src/grdfft.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ Parameters
"""
function grdfft(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

(GMTver < v"6.4.0") && (@warn("Sorry but you need at least GMT 6.4 to use this module. Previous version had a bug that prevented Julia wrapper to work."); return nothing)
d = init_module(false, kwargs...)[1] # Also checks if the user wants ONLY the HELP mode

cmd, = parse_common_opts(d, "", [:G :V_params :f])
Expand Down
55 changes: 39 additions & 16 deletions src/potential/gmtgravmag3d.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
gmtgravmag3d(cmd0::String=""; kwargs...)
gravmag3d(cmd0::String=""; kwargs...)
Compute the gravity/magnetic anomaly of a 3-D body by the method of Okabe.
Expand Down Expand Up @@ -32,9 +32,9 @@ Parameters
- **S** | **radius** :: [Type => Number]
Set search radius in km (valid only in the two grids mode OR when `thickness`) [Default = 30 km].
- **Tv** | **index** :: [Type => Str]
- **Tr** | **raw_triang** :: [Type => Str]
- **Ts** | **stl** :: [Type => Str]
- **T+v** | **index** :: [Type => Str]
- **T+r** | **raw_triang** :: [Type => Str]
- **T+s** | **stl** :: [Type => Str]
Gives names of a xyz and vertex (ndex="vert_file") files defining a close surface.
- **Z** | **z_level** | **reference_level** :: [Type => Number]
Expand All @@ -47,39 +47,62 @@ Parameters
imshow(G)
```
"""
function gmtgravmag3d(cmd0::String="", arg1=nothing; kwargs...)
#
gravmag3d(cmd0::String; kwargs...) = gravmag3d_helper(cmd0, nothing; kwargs...)
gravmag3d(arg1; kwargs...) = gravmag3d_helper("", arg1; kwargs...)
gravmag3d(; kwargs...) = gravmag3d_helper("", nothing; kwargs...)

function gravmag3d_helper(cmd0::String, arg1; kwargs...)
d = init_module(false, kwargs...)[1] # Also checks if the user wants ONLY the HELP mode
gravmag3d_helper(cmd0, arg1, d)
end
#

# ---------------------------------------------------------------------------------------------------
function gravmag3d_helper(cmd0::String, arg1, d::Dict{Symbol,Any})

cmd = parse_common_opts(d, "", [:G :RIr :V_params :bi :f])[1]
cmd = parse_these_opts(cmd, d, [[:C :density], [:E :thickness], [:F :track], [:G :grid :outgrid],
cmd = parse_these_opts(cmd, d, [[:C :density], [:E :thickness], [:G :grid :outgrid],
[:L :z_obs :observation_level], [:S :radius], [:Z :reference_level]])
cmd = add_opt(d, cmd, "H", [:H :mag_params], (field_dec="", field_dip="", mag="", mag_dec="", mag_dip=""))

arg2 = nothing;
arg2, arg3 = nothing, nothing;
if ((val = find_in_dict(d, [:Tv :index], true, "GMTdataset | Mx3 array | String")[1]) !== nothing && val != "") opt = " -Tv"
elseif ((val = find_in_dict(d, [:Tr :raw_triang], true, "GMTdataset | Mx3 array | String")[1]) !== nothing && val != "") opt = " -Tr"
elseif ((val = find_in_dict(d, [:Ts :stl :STL], true, "String (file name)")[1]) !== nothing && val != "") opt = " -Ts"
elseif ((val = find_in_dict(d, [:M :body], false, "String | NamedTuple")[1]) !== nothing && val != "")
cmd = add_opt(d, cmd, "M", [:M :body], (shape="+s", params=","))
if (isa(val, Tuple))
cmd = add_opt(Dict(:body => val[1]), cmd, "M", [:M :body], (shape="+s", params=","))
for n = 2:numel(val)
cmd = add_opt(Dict(:body => val[n]), cmd, "", [:M :body], (shape="+s", params=","))
end
delete!(d, [:M :body])
else
cmd = add_opt(d, cmd, "M", [:M :body], (shape="+s", params=","))
end
opt = " -Ts"
elseif (SHOW_KWARGS[1]) opt = ""
else error("Missing one of 'index', 'raw_triang' or 'str' data")
end

if (opt != " -Ts") # The STL format can only be requested via file
cmd *= opt
if (opt != " -Ts") # The STL format can only be requested via file (so we can keep testing that with old syntax)
if (isa(val, Array{<:Real}) || isa(val, GDtype))
(arg1 === nothing) ? arg1 = val : arg2 = val # Find the free slot
cmd *= " -T+" * opt[end] # New -T syntax
else
cmd *= arg2str(val)
end
if (opt == " -Tv")
(!occursin(" -I", cmd)) && error("For grid output MUST specify grid increment ('I' or 'inc')")
cmd *= opt # Old syntax. This should be removed when we rebuild GMT_jll
end
end

if ((val = find_in_dict(d, [:F :track], true, "GMTdataset | Mx2 array | String")[1]) !== nothing && val != "")
cmd *= " -F"
isa(val, String) ? (cmd *= val) : (arg1 === nothing) ? arg1 = val : (arg2 === nothing) ? arg2 = val : arg3 = val
end

(!occursin(" -F", cmd) && !occursin(" -G", cmd)) && (cmd *= " -G")

(opt != " -Tv") && return finish_PS_module(d, "gmtgravmag3d " * cmd, "", true, false, false, arg1, arg2)
common_grd(d, cmd0, cmd, "gmtgravmag3d ", arg1, arg2) # Finish build cmd and run it
end
finish_PS_module(d, "gmtgravmag3d " * cmd, "", true, false, false, arg1, arg2, arg3)
end

const gmtgravmag3d = gravmag3d # Alias
4 changes: 2 additions & 2 deletions src/potential/grdgravmag3d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ Parameters
- **S** | **radius** :: [Type => Number]
Set search radius in km (valid only in the two grids mode OR when `thickness`) [Default = 30 km].
- **Z** | **z_level** | **reference_level** :: [Type => Number]
- **Z** | **level** | **reference_level** :: [Type => Number]
Level of reference plane [Default = 0].
- $(opt_V)
Expand All @@ -57,7 +57,7 @@ function grdgravmag3d(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

cmd::String = parse_common_opts(d, "", [:I :G :R :V_params :f :x])[1]
cmd = parse_these_opts(cmd, d, [[:E :thickness], [:Q :pad], [:L :z_obs :observation_level], [:S :radius]])
opt_Z = add_opt(d, "", "Z", [:Z :z_level :reference_level])
opt_Z = add_opt(d, "", "Z", [:Z :level :reference_level], (bottom="_b", top="_t"))
if (opt_Z != "")
if (opt_Z[4] == 't') cmd *= " -Zt"
elseif (opt_Z[4] == 'b') cmd *= " -Zb"
Expand Down
15 changes: 12 additions & 3 deletions test/test_B-GMTs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,9 +66,18 @@
gmtconvert([1.1 2; 3 4], o=0)

println(" GMTGRAVMAG3D")
gmtgravmag3d(M=(shape=:prism, params=(1,1,1,5)), I=1.0, R="-15/15/-15/15", H="10/60/10/-10/40", Vd=dbg2);
@test_throws ErrorException("Missing one of 'index', 'raw_triang' or 'str' data") gmtgravmag3d(I=1.0);
@test_throws ErrorException("For grid output MUST specify grid increment ('I' or 'inc')") gmtgravmag3d(Tv=true);
#gmtgravmag3d(M=(shape=:prism, params=(1,1,1,5)), I=1.0, R="-15/15/-15/15", H="10/60/10/-10/40", Vd=dbg2);
#@test isa(G, GMTgrid)
#D = gravmag3d(region="-15/15/-15/15", I=0.1, mag_params="10/60/10/-10/40", body=(shape=:prism, params="1/1/1/-5/-10/1"), F=[-14 -14; 14 14]);
#@test isa(D, GMTdataset)
#D = grd2xyz("@earth_relief_10m_g", R="-41:50/-41:20/47:30/47:50", f=:g);
#Dtri = triangulate(D);
#D = gmtgravmag3d(D, index=Dtri, C=1700, R="-41.833/-41.333/47.5/47.833", f=:g, Z=-4300, F=[-41.793 47.552; -41.567 47.672; -41.357 47.809]);
#@test isa(D, GMTdataset)
#G = gmtgravmag3d(D, index=Dtri, C=1700, R="-41.833/-41.333/47.5/47.833", f=:g, Z=-4300, I="1m");
#@test isa(G, GMTgrid)
#@test_throws ErrorException("Missing one of 'index', 'raw_triang' or 'str' data") gmtgravmag3d(I=1.0);
#@test_throws ErrorException("For grid output MUST specify grid increment ('I' or 'inc')") gmtgravmag3d(Tv=true);

println(" GRDGRAVMAG3D")
grdgravmag3d("@earth_relief_10m", region=(-12.5,-10,35.5,37.5), density=1700, inc=0.05, pad=0.5, z_level=:b, f=:g, Vd=dbg2)
Expand Down
3 changes: 2 additions & 1 deletion test/test_misc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,7 @@
image_alpha!(img, alpha_band=mask, burn=:red)
mask[1] = 100; # Force variable mask
image_alpha!(img, alpha_band=mask, burn=(0,255,0))
@info "after image_alpha!"

I1 = mat2img(rand(UInt8, 16,16)); I2 = mat2img(rand(UInt8, 16,16)); I3 = mat2img(rand(UInt8, 16,16));
grays2rgb(I1,I2,I3);
Expand All @@ -234,8 +235,8 @@
GMT.transpcmap!(I, false)
@info "before image_cpt!"
image_cpt!(I, clear=true)

@info "after image_cpt!"

GMT.linspace(1,1,100);
GMT.logspace(1,5);
GMT.fakedata(50,1);
Expand Down

0 comments on commit e4a935b

Please sign in to comment.