fix plotting topography; add option to read topo files

JuliaGeodynamics · Jan 7, 2025 · 41ec0ba · 41ec0ba
1 parent 4275414
commit 41ec0ba
Showing 1 changed file with 28 additions and 15 deletions.
diff --git a/ext/PlotsExt.jl b/ext/PlotsExt.jl
@@ -12,18 +12,23 @@ println("Adding Plots.jl plotting extensions for LaMEM")
 """
     plot_cross_section(model::Model, args...; field::Symbol=:phase,  
                         timestep::Union{Nothing, Int64}=nothing,
-                        dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
+                        dim=1, x=nothing, y=nothing, z=nothing, 
+                        aspect_ratio::Union{Real, Symbol}=:equal,
+                        surf=false)
    
 This plots a cross-section through the LaMEM `model`, of the field  `field`. If that is a vector or tensor field, specify `dim` to indicate which of the fields you want to see.
 If the keyword `timestep` is specified, it will load a timestep 
 """
 function plot_cross_section(model::Model, args...; field::Symbol=:phase, 
         timestep::Union{Nothing, Int64}=0, dim=1, 
-        x=nothing, y=nothing, z=nothing, 
-        aspect_ratio::Union{Real, Symbol}=:equal)
+        x=nothing, 
+        y=nothing, 
+        z=nothing, 
+        aspect_ratio::Union{Real, Symbol}=:equal,
+        surf=false)
 
     # load a particular timestep
-    data_cart, time_val = read_LaMEM_timestep(model,timestep)
+    data_cart, time_val = read_LaMEM_timestep(model,timestep, surf=surf)
 
     if isnothing(x) && isnothing(y) && isnothing(z)
         x = mean(extrema(model.Grid.Grid.X))
@@ -42,7 +47,12 @@ end
 """
     plot_cross_section(data::CartData, args...; field::Symbol=:phase,  
                         title_str="",
-                        dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
+                        dim=1, 
+                        x=nothing, 
+                        y=nothing, 
+                        z=nothing, 
+                        aspect_ratio::Union{Real, Symbol}=:equal
+                        surf=false)
    
 This plots a cross-section through a `CartData` dataset `data` of the field  `field`, typically read in from a LaMEM simulation.
 If that is a vector or tensor field, specify `dim` to indicate which of the fields you want to see.
@@ -51,14 +61,15 @@ function plot_cross_section(data::CartData, args...; field::Symbol=:phase,
         dim=1, 
         x=nothing, y=nothing, z=nothing, 
         title_str="",
-        aspect_ratio::Union{Real, Symbol}=:equal)
+        aspect_ratio::Union{Real, Symbol}=:equal,
+        surf=false)
 
 
     if isnothing(x) && isnothing(y) && isnothing(z)
         x = mean(extrema(data.x.val))
     end
 
-    data_tuple, axes_str = cross_section(data, field; x=x, y=y, z=z)
+    data_tuple, axes_str = cross_section(data, field; x=x, y=y, z=z, surf=surf)
 
     if isa(data_tuple.data, Array)
         data_field = data_tuple.data
@@ -99,19 +110,21 @@ end
 
 
 """
-    plot_topo(topo::CartData, args...)
+    plot_topo(topo::CartData; kwargs...)
 
 Simple function to plot the topography 
 """
 function plot_topo(topo::CartData; kwargs...)
 
-    hm = Plots.heatmap(topo.x.val[:,1,1], topo.y.val[1,:,1], topo.fields.Topography[:,:,1]'; 
-                aspect_ratio=:equal, 
-                xlabel="x",
-                ylabel="y",
-                title="topography [km]",
-                colormap=:oleron,
-                kwargs...)
+    hm = Plots.heatmap( topo.x.val[:,1,1], 
+                        topo.y.val[1,:,1], 
+                        topo.fields.topography[:,:,1]'; 
+                        aspect_ratio=:equal, 
+                        xlabel="x",
+                        ylabel="y",
+                        title="topography [km]",
+                        colormap=:oleron,
+                        kwargs...)
 
     return hm
 end