From 41ec0ba7ef69066a70de255fb0e040fd9a0093b2 Mon Sep 17 00:00:00 2001
From: Boris Kaus <kaus@uni-mainz.de>
Date: Tue, 7 Jan 2025 23:16:50 +0100
Subject: [PATCH] fix plotting topography; add option to read topo files

---
 ext/PlotsExt.jl | 43 ++++++++++++++++++++++++++++---------------
 1 file changed, 28 insertions(+), 15 deletions(-)

diff --git a/ext/PlotsExt.jl b/ext/PlotsExt.jl
index 29a1faa..67bc441 100644
--- 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