Skip to content

Commit 492f31e

Browse files
haakon-ehaakon-e
committed
migrate plot_timeheight to Makie
1 parent afc4fe3 commit 492f31e

File tree

1 file changed

+90
-114
lines changed

1 file changed

+90
-114
lines changed

‎test/plotting_utils.jl‎

Lines changed: 90 additions & 114 deletions
Original file line numberDiff line numberDiff line change
@@ -447,80 +447,87 @@ function plot_timeheight_p3(nc_data_file, precip; output = "output")
447447
Plots.png(p, joinpath(path, "timeheight.png"))
448448
end
449449

450+
function make_tz_hm(i, j, t, z, var; imax = 3, title = "", kw...)
451+
fig = current_figure()
452+
ax = Axis(fig[i, j];
453+
xlabel = i == imax ? "time [s]" : "",
454+
ylabel = j == 1 ? "z [m]" : "",
455+
yticklabelsvisible = j == 1,
456+
xticklabelsvisible = i == 3,
457+
title,
458+
)
459+
hm = CairoMakie.heatmap!(ax, t, z, var; colormap = :BuPu, kw...)
460+
Makie.Colorbar(fig[i, j], hm;
461+
halign = 1, valign = 0,
462+
height = Makie.Relative(0.3), width = 10, tellwidth = false,
463+
ticks = WilkinsonTicks(3; k_min = 3, k_max = 6),
464+
)
465+
end
466+
450467
function plot_timeheight(nc_data_file; output = "output", mixed_phase = true, pysdm = false)
451468
path = joinpath(@__DIR__, output)
452469
mkpath(path)
453470

454471
ds = NC.NCDataset(joinpath(@__DIR__, nc_data_file))
455472
if pysdm
456-
t_plt = Array(ds["time"])
457-
z_plt = Array(ds["height"])
458-
q_liq_plt = transpose(Array(ds["cloud water mixing ratio"]))
459-
q_rai_plt = transpose(Array(ds["rain water mixing ratio"]))
460-
q_ice_plt = transpose(Array(ds["rain water mixing ratio"])) * FT(0)
461-
q_sno_plt = transpose(Array(ds["rain water mixing ratio"])) * FT(0)
462-
q_vap = transpose(Array(ds["water_vapour_mixing_ratio"])) * 1e3
463-
q_tot_plt = q_vap + q_liq_plt
464-
N_aer_plt = transpose(Array(ds["na"]))
465-
N_liq_plt = transpose(Array(ds["nc"]))
466-
N_rai_plt = transpose(Array(ds["nr"]))
467-
SN_liq_act_plt = transpose(Array(ds["activating"]))
473+
t = Array(ds["time"])
474+
z = Array(ds["height"])
475+
q_liq = Array(ds["cloud water mixing ratio"])
476+
q_rai = Array(ds["rain water mixing ratio"])
477+
q_ice = Array(ds["rain water mixing ratio"]) * FT(0)
478+
q_sno = Array(ds["rain water mixing ratio"]) * FT(0)
479+
q_vap = Array(ds["water_vapour_mixing_ratio"]) * 1e3
480+
q_tot = q_vap + q_liq
481+
N_aer = Array(ds["na"])
482+
N_liq = Array(ds["nc"])
483+
N_rai = Array(ds["nr"])
484+
SN_liq_act = transpose(Array(ds["activating"]))
468485
else
469-
t_plt = Array(ds.group["profiles"]["t"])
470-
z_plt = Array(ds.group["profiles"]["zc"])
471-
q_tot_plt = Array(ds.group["profiles"]["q_tot"])
472-
q_liq_plt = Array(ds.group["profiles"]["q_liq"])
473-
q_ice_plt = Array(ds.group["profiles"]["q_ice"])
474-
q_rai_plt = Array(ds.group["profiles"]["q_rai"])
475-
q_sno_plt = Array(ds.group["profiles"]["q_sno"])
476-
N_aer_plt = Array(ds.group["profiles"]["N_aer"])
477-
N_liq_plt = Array(ds.group["profiles"]["N_liq"])
478-
N_rai_plt = Array(ds.group["profiles"]["N_rai"])
479-
SN_liq_act_plt = Array(ds.group["profiles"]["SN_liq_act"])
486+
data = ds.group["profiles"]
487+
t = Array(data["t"])
488+
z = Array(data["zc"])
489+
q_tot = permutedims(Array(data["q_tot"]))
490+
q_liq = permutedims(Array(data["q_liq"]))
491+
q_ice = permutedims(Array(data["q_ice"]))
492+
q_rai = permutedims(Array(data["q_rai"]))
493+
q_sno = permutedims(Array(data["q_sno"]))
494+
N_aer = permutedims(Array(data["N_aer"]))
495+
N_liq = permutedims(Array(data["N_liq"]))
496+
N_rai = permutedims(Array(data["N_rai"]))
497+
SN_liq_act = permutedims(Array(data["SN_liq_act"]))
480498
end
481-
#! format: off
482-
p1 = Plots.heatmap(t_plt, z_plt, q_tot_plt .* 1e3, title = "q_tot [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
483-
p2 = Plots.heatmap(t_plt, z_plt, q_liq_plt .* 1e3, title = "q_liq [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
484-
p3 = Plots.heatmap(t_plt, z_plt, q_ice_plt .* 1e3, title = "q_ice [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
485-
p4 = Plots.heatmap(t_plt, z_plt, q_rai_plt .* 1e3, title = "q_rai [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
486-
p5 = Plots.heatmap(t_plt, z_plt, q_sno_plt .* 1e3, title = "q_sno [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
487-
p6 = Plots.heatmap(t_plt, z_plt, N_aer_plt .* 1e-6, title = "N_aer [1/cm^3]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
488-
p7 = Plots.heatmap(t_plt, z_plt, N_liq_plt .* 1e-6, title = "N_liq [1/cm^3]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
489-
p8 = Plots.heatmap(t_plt, z_plt, N_rai_plt .* 1e-6, title = "N_rai [1/cm^3]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
490-
p9 = Plots.heatmap(t_plt, z_plt, SN_liq_act_plt .* 1e-6, title = "Activation [1/cm^3/s]", xlabel = "time [s]", ylabel = "z [m]", color = :BuPu)
491-
#! format: on
499+
500+
kg_to_g = 1e3
501+
m⁻³_to_cm⁻³ = 1e-6
502+
492503
if mixed_phase
493-
p = Plots.plot(
494-
p1,
495-
p2,
496-
p3,
497-
p4,
498-
p5,
499-
p6,
500-
p7,
501-
p8,
502-
p9,
503-
size = (1200.0, 1200.0),
504-
bottom_margin = 30.0 * Plots.PlotMeasures.px,
505-
left_margin = 30.0 * Plots.PlotMeasures.px,
506-
layout = (3, 3),
507-
)
504+
fig = Figure(size = (1200, 1200))
505+
make_tz_hm(1, 1, t, z, q_tot .* kg_to_g; colorrange = (8, 15), title = "q_tot [g/kg]")
506+
make_tz_hm(1, 2, t, z, q_liq .* kg_to_g; colorrange = (0, 1), title = "q_liq [g/kg]")
507+
make_tz_hm(1, 3, t, z, q_ice .* kg_to_g; colorrange = (0, 0.25), title = "q_ice [g/kg]")
508+
make_tz_hm(2, 1, t, z, q_rai .* kg_to_g; colorrange = (0, 0.25), title = "q_rai [g/kg]")
509+
make_tz_hm(2, 2, t, z, q_sno .* kg_to_g; title = "q_sno [g/kg]")
510+
make_tz_hm(2, 3, t, z, N_aer .* m⁻³_to_cm⁻³; colorrange = (0, 100), title = "N_aer [1/cm³]")
511+
make_tz_hm(3, 1, t, z, N_liq .* m⁻³_to_cm⁻³; colorrange = (0, 50), title = "N_liq [1/cm³]")
512+
make_tz_hm(3, 2, t, z, N_rai .* m⁻³_to_cm⁻³; colorrange = (0, 1), title = "N_rai [1/cm³]")
513+
make_tz_hm(3, 3, t, z, SN_liq_act .* m⁻³_to_cm⁻³; title = "Activation [1/cm³/s]")
514+
508515
else
509-
p = Plots.plot(
510-
p1,
511-
p2,
512-
p4,
513-
p6,
514-
p7,
515-
p8,
516-
p9,
517-
size = (1200.0, 900.0),
518-
bottom_margin = 30.0 * Plots.PlotMeasures.px,
519-
left_margin = 30.0 * Plots.PlotMeasures.px,
520-
layout = (3, 3),
521-
)
516+
fig = Figure(size = (1200, 600))
517+
make_tz_hm(1, 1, t, z, q_tot .* kg_to_g; colorrange = (8, 15), title = "q_tot [g/kg]")
518+
make_tz_hm(1, 2, t, z, q_liq .* kg_to_g; colorrange = (0, 1), title = "q_liq [g/kg]")
519+
make_tz_hm(1, 3, t, z, q_rai .* kg_to_g; colorrange = (0, 0.25), title = "q_rai [g/kg]")
520+
make_tz_hm(2, 1, t, z, N_aer .* m⁻³_to_cm⁻³; colorrange = (0, 100), title = "N_aer [1/cm³]")
521+
make_tz_hm(2, 2, t, z, N_liq .* m⁻³_to_cm⁻³; colorrange = (0, 50), title = "N_liq [1/cm³]")
522+
make_tz_hm(2, 3, t, z, N_rai .* m⁻³_to_cm⁻³; colorrange = (0, 1), title = "N_rai [1/cm³]")
523+
make_tz_hm(3, 1, t, z, SN_liq_act .* m⁻³_to_cm⁻³; title = "Activation [1/cm³/s]")
522524
end
523-
Plots.png(p, joinpath(path, "timeheight.png"))
525+
fig
526+
527+
axs = filter(x -> x isa Axis, contents(fig[:, :]))
528+
linkaxes!(axs...)
529+
save(joinpath(path, "timeheight.png"), fig)
530+
nothing
524531
end
525532

526533
function plot_timeheight_no_ice_snow(nc_data_file; output = "output", pysdm = false)
@@ -529,58 +536,27 @@ function plot_timeheight_no_ice_snow(nc_data_file; output = "output", pysdm = fa
529536

530537
ds = NC.NCDataset(joinpath(@__DIR__, nc_data_file))
531538
if pysdm
532-
t_plt = Array(ds["time"])
533-
z_plt = Array(ds["height"])
534-
q_liq_plt = transpose(Array(ds["cloud water mixing ratio"])) / 1e3
535-
q_rai_plt = transpose(Array(ds["rain water mixing ratio"])) / 1e3
536-
q_vap = transpose(Array(ds["water_vapour_mixing_ratio"]))
537-
q_tot_plt = q_vap + q_liq_plt
539+
t = Array(ds["time"])
540+
z = Array(ds["height"])
541+
q_liq = Array(ds["cloud water mixing ratio"]) / 1e3
542+
q_rai = Array(ds["rain water mixing ratio"]) / 1e3
543+
q_vap = Array(ds["water_vapour_mixing_ratio"])
544+
q_tot = q_vap + q_liq
538545
else
539-
t_plt = Array(ds.group["profiles"]["t"])
540-
z_plt = Array(ds.group["profiles"]["zc"])
541-
q_tot_plt = Array(ds.group["profiles"]["q_tot"])
542-
q_liq_plt = Array(ds.group["profiles"]["q_liq"])
543-
q_rai_plt = Array(ds.group["profiles"]["q_rai"])
546+
t = Array(ds.group["profiles"]["t"])
547+
z = Array(ds.group["profiles"]["zc"])
548+
q_tot = permutedims(Array(ds.group["profiles"]["q_tot"]))
549+
q_liq = permutedims(Array(ds.group["profiles"]["q_liq"]))
550+
q_rai = permutedims(Array(ds.group["profiles"]["q_rai"]))
544551
end
545552

546-
p1 = Plots.heatmap(
547-
t_plt,
548-
z_plt,
549-
q_tot_plt .* 1e3,
550-
title = "q_tot [g/kg]",
551-
xlabel = "time [s]",
552-
ylabel = "z [m]",
553-
color = :BuPu,
554-
clims = (0, 1),
555-
)
556-
p2 = Plots.heatmap(
557-
t_plt,
558-
z_plt,
559-
q_liq_plt .* 1e3,
560-
title = "q_liq [g/kg]",
561-
xlabel = "time [s]",
562-
ylabel = "z [m]",
563-
color = :BuPu,
564-
clims = (0, 1),
565-
)
566-
p3 = Plots.heatmap(
567-
t_plt,
568-
z_plt,
569-
q_rai_plt .* 1e3,
570-
title = "q_rai [g/kg]",
571-
xlabel = "time [s]",
572-
ylabel = "z [m]",
573-
color = :BuPu,
574-
clims = (0, 0.25),
575-
)
576-
p = Plots.plot(
577-
p1,
578-
p2,
579-
p3,
580-
size = (1200.0, 300.0),
581-
bottom_margin = 30.0 * Plots.PlotMeasures.px,
582-
left_margin = 30.0 * Plots.PlotMeasures.px,
583-
layout = (1, 3),
584-
)
585-
Plots.png(p, joinpath(path, "timeheight.png"))
553+
kg_to_g = 1e3
554+
fig = Figure(size = (1200, 300))
555+
make_tz_hm(1, 1, t, z, q_tot .* kg_to_g; title = "q_tot [g/kg]", imax = 1, colorrange = (0, 1))
556+
make_tz_hm(1, 2, t, z, q_liq .* kg_to_g; title = "q_liq [g/kg]", imax = 1, colorrange = (0, 1))
557+
make_tz_hm(1, 3, t, z, q_rai .* kg_to_g; title = "q_rai [g/kg]", imax = 1, colorrange = (0, 0.25))
558+
axs = filter(x -> x isa Axis, contents(fig[:, :]))
559+
linkaxes!(axs...)
560+
save(joinpath(path, "timeheight.png"), fig)
561+
nothing
586562
end

0 commit comments

Comments
 (0)