Skip to content

Commit 489c92c

Browse files
haakon-ehaakon-e
committed
migrate plot_timeheight to Makie
1 parent 2d09ae4 commit 489c92c

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
@@ -460,80 +460,87 @@ function plot_timeheight_p3(nc_data_file, precip; output = "output")
460460
Plots.png(p, joinpath(path, "timeheight.png"))
461461
end
462462

463+
function make_tz_hm(i, j, t, z, var; imax = 3, title = "", kw...)
464+
fig = current_figure()
465+
ax = Axis(fig[i, j];
466+
xlabel = i == imax ? "time [s]" : "",
467+
ylabel = j == 1 ? "z [m]" : "",
468+
yticklabelsvisible = j == 1,
469+
xticklabelsvisible = i == 3,
470+
title,
471+
)
472+
hm = CairoMakie.heatmap!(ax, t, z, var; colormap = :BuPu, kw...)
473+
Makie.Colorbar(fig[i, j], hm;
474+
halign = 1, valign = 0,
475+
height = Makie.Relative(0.3), width = 10, tellwidth = false,
476+
ticks = WilkinsonTicks(3; k_min = 3, k_max = 6),
477+
)
478+
end
479+
463480
function plot_timeheight(nc_data_file; output = "output", mixed_phase = true, pysdm = false)
464481
path = joinpath(@__DIR__, output)
465482
mkpath(path)
466483

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

539546
function plot_timeheight_no_ice_snow(nc_data_file; output = "output", pysdm = false)
@@ -542,58 +549,27 @@ function plot_timeheight_no_ice_snow(nc_data_file; output = "output", pysdm = fa
542549

543550
ds = NC.NCDataset(joinpath(@__DIR__, nc_data_file))
544551
if pysdm
545-
t_plt = Array(ds["time"])
546-
z_plt = Array(ds["height"])
547-
q_liq_plt = transpose(Array(ds["cloud water mixing ratio"])) / 1e3
548-
q_rai_plt = transpose(Array(ds["rain water mixing ratio"])) / 1e3
549-
q_vap = transpose(Array(ds["water_vapour_mixing_ratio"]))
550-
q_tot_plt = q_vap + q_liq_plt
552+
t = Array(ds["time"])
553+
z = Array(ds["height"])
554+
q_liq = Array(ds["cloud water mixing ratio"]) / 1e3
555+
q_rai = Array(ds["rain water mixing ratio"]) / 1e3
556+
q_vap = Array(ds["water_vapour_mixing_ratio"])
557+
q_tot = q_vap + q_liq
551558
else
552-
t_plt = Array(ds.group["profiles"]["t"])
553-
z_plt = Array(ds.group["profiles"]["zc"])
554-
q_tot_plt = Array(ds.group["profiles"]["q_tot"])
555-
q_liq_plt = Array(ds.group["profiles"]["q_liq"])
556-
q_rai_plt = Array(ds.group["profiles"]["q_rai"])
559+
t = Array(ds.group["profiles"]["t"])
560+
z = Array(ds.group["profiles"]["zc"])
561+
q_tot = permutedims(Array(ds.group["profiles"]["q_tot"]))
562+
q_liq = permutedims(Array(ds.group["profiles"]["q_liq"]))
563+
q_rai = permutedims(Array(ds.group["profiles"]["q_rai"]))
557564
end
558565

559-
p1 = Plots.heatmap(
560-
t_plt,
561-
z_plt,
562-
q_tot_plt .* 1e3,
563-
title = "q_tot [g/kg]",
564-
xlabel = "time [s]",
565-
ylabel = "z [m]",
566-
color = :BuPu,
567-
clims = (0, 1),
568-
)
569-
p2 = Plots.heatmap(
570-
t_plt,
571-
z_plt,
572-
q_liq_plt .* 1e3,
573-
title = "q_liq [g/kg]",
574-
xlabel = "time [s]",
575-
ylabel = "z [m]",
576-
color = :BuPu,
577-
clims = (0, 1),
578-
)
579-
p3 = Plots.heatmap(
580-
t_plt,
581-
z_plt,
582-
q_rai_plt .* 1e3,
583-
title = "q_rai [g/kg]",
584-
xlabel = "time [s]",
585-
ylabel = "z [m]",
586-
color = :BuPu,
587-
clims = (0, 0.25),
588-
)
589-
p = Plots.plot(
590-
p1,
591-
p2,
592-
p3,
593-
size = (1200.0, 300.0),
594-
bottom_margin = 30.0 * Plots.PlotMeasures.px,
595-
left_margin = 30.0 * Plots.PlotMeasures.px,
596-
layout = (1, 3),
597-
)
598-
Plots.png(p, joinpath(path, "timeheight.png"))
566+
kg_to_g = 1e3
567+
fig = Figure(size = (1200, 300))
568+
make_tz_hm(1, 1, t, z, q_tot .* kg_to_g; title = "q_tot [g/kg]", imax = 1, colorrange = (0, 1))
569+
make_tz_hm(1, 2, t, z, q_liq .* kg_to_g; title = "q_liq [g/kg]", imax = 1, colorrange = (0, 1))
570+
make_tz_hm(1, 3, t, z, q_rai .* kg_to_g; title = "q_rai [g/kg]", imax = 1, colorrange = (0, 0.25))
571+
axs = filter(x -> x isa Axis, contents(fig[:, :]))
572+
linkaxes!(axs...)
573+
save(joinpath(path, "timeheight.png"), fig)
574+
nothing
599575
end

0 commit comments

Comments
 (0)