Skip to content

Commit f408d68

Browse files
haakon-ehaakon-e
committed
migrate plot_timeheight to Makie
1 parent 16eb99c commit f408d68

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

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

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

540547
function plot_timeheight_no_ice_snow(nc_data_file; output = "output", pysdm = false)
@@ -543,58 +550,27 @@ function plot_timeheight_no_ice_snow(nc_data_file; output = "output", pysdm = fa
543550

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

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

0 commit comments

Comments
 (0)