Merge pull request #121 from jbisits/distributionfixes

Fixes for distribution plotting

Author: jbisits

ext/TLDNSMakieExt/plotrecipes.jl

@@ -65,7 +65,7 @@ function TLDNS.animate_2D_field(field_timeseries::FieldTimeSeries, field_name::A
 end
 function TLDNS.animate_2D_field(rs::Raster, xslice::Int64, yslice::Int64; colormap = :thermal,
                                 colorrange = nothing, highclip = nothing, lowclip = nothing,
-                                aspect_ratio = 1, profile_xlims = nothing)
+                                aspect_ratio = 1, vline = nothing)
 
     x, z, t = lookup(rs, :xC), lookup(rs, :zC), lookup(rs, :Ti)
     field_name = string(rs.name)
@@ -94,8 +94,11 @@ function TLDNS.animate_2D_field(rs::Raster, xslice::Int64, yslice::Int64; colorm
     ax[2].ylabel = "z"
     ax[2].aspect = aspect_ratio
     ax[2].xaxisposition = :top
-    if !isnothing(profile_xlims)
-        xlims!(ax[2], profile_xlims)
+    ax[2].xticklabelrotation = π / 4
+    if !isnothing(vline)
+        vlines!(ax[2], vline, linestyle = :dash, color = :red,
+                label = "Predicted maximum density")
+        axislegend(ax[2])
     end
 
     linkyaxes!(ax[1], ax[2])
@@ -114,23 +117,35 @@ end
 """
     function animate_volume_distributions(rs::Raster)
 Animate the volume distribution from each saved snapshot of data for the variable saved as
-a `Raster`.
+a `Raster`. Optional arguments:
+- `edges` for the bins, in `nothing` the edges for the first fitted `Histogram` are used for
+the subsequent histograms;
+- `unit` for the variable that is being plotted - must be a `String`.
 """
-function TLDNS.animate_volume_distributions(rs::Raster, edges::AbstractVector; unit = nothing)
-
-    rs_series = slice(rs, Ti)
-    t = lookup(rs_series, Ti)
-    rs_series_hist = RasterLayerHistogram.(rs_series, edges)
+function TLDNS.animate_volume_distributions(rs::Raster; edges = nothing, unit = nothing)
+
+    t = lookup(rs, Ti)
+    rs_series_hist = Vector{RasterLayerHistogram}(undef, length(t))
+    for i ∈ eachindex(t)
+        if i == 1
+            rs_series_hist[i] = isnothing(edges) ? RasterLayerHistogram(rs[:, :, :, i]) :
+                                                   RasterLayerHistogram(rs[:, :, :, i], edges)
+        else
+            rs_series_hist[i] = RasterLayerHistogram(rs[:, :, :, i],
+                                                     rs_series_hist[1].histogram.edges[1])
+        end
+    end
     n = Observable(1)
-    @lift rs_series_hist[$n]
+    dₜ = @lift rs_series_hist[$n]
     time_title = @lift @sprintf("t=%1.2f minutes", t[$n] / 60)
 
     fig = Figure(size = (500, 500))
     xlabel = isnothing(unit) ? string(rs.name) : string(rs.name) * unit
+    xlimits = isnothing(edges) ? rs_series_hist[1].histogram.edges[1] : edges
     ylabel = "Volume (m³)"
     ax = Axis(fig[1, 1], title = time_title; xlabel, ylabel)
-
-    plot!(ax, rs_series_hist[1], color = :steelblue)
+    xlims!(ax, xlimits[1], xlimits[end])
+    plot!(ax, dₜ, color = :steelblue)
 
     frames = eachindex(t)
     record(fig, joinpath(pwd(), string(rs.name) * "_vd.mp4"), frames, framerate=8) do i