Skip to content

Commit c00070a

Browse files
briochemcbriochemc
committed
Add MLD plot
1 parent fa8df9a commit c00070a

File tree

2 files changed

+301
-2
lines changed

2 files changed

+301
-2
lines changed

‎src/plot_ACCESS_MLDs.jl‎

Lines changed: 250 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,250 @@
1+
# # qsub -I -P xv83 -q express -l mem=47GB -l storage=scratch/gh0+scratch/xv83 -l walltime=01:00:00 -l ncpus=12
2+
# # This is Fig. 1 in Pasquier et al. (GRL, 2025)
3+
4+
# using Pkg
5+
# Pkg.activate(".")
6+
# Pkg.instantiate()
7+
8+
# using OceanTransportMatrixBuilder
9+
# using NetCDF
10+
# using YAXArrays
11+
# using DataFrames
12+
# using DimensionalData
13+
# # using SparseArrays
14+
# # using LinearAlgebra
15+
# using Unitful
16+
# using Unitful: s, yr
17+
# try
18+
# using CairoMakie
19+
# catch
20+
# using CairoMakie
21+
# end
22+
# using GeoMakie
23+
# using Interpolations
24+
# using OceanBasins
25+
# using Statistics
26+
# using NaNStatistics
27+
# using StatsBase
28+
# using FileIO
29+
# using Contour
30+
# using GeometryBasics
31+
# using GeometryOps
32+
# using LibGEOS
33+
# using Format
34+
35+
# model = "ACCESS-ESM1-5"
36+
37+
# time_window_1850s = "Jan1850-Dec1859"
38+
# time_window_2030s = "Jan2030-Dec2039"
39+
# time_window_2090s = "Jan2090-Dec2099"
40+
# time_windows = ["1850s", "2030s", "2090s"]
41+
# experiments = ["historical", "ssp370", "ssp370"]
42+
# experiments2 = ["historical", "SSP3-7.0", "SSP3-7.0"]
43+
44+
# members = ["r$(r)i1p1f1" for r in 1:40]
45+
# # members = ["r$(r)i1p1f1" for r in 1:3]
46+
47+
# # Load areacello and volcello for grid geometry
48+
# fixedvarsinputdir = "/scratch/xv83/TMIP/data/$model"
49+
# volcello_ds = open_dataset(joinpath(fixedvarsinputdir, "volcello.nc"))
50+
# areacello_ds = open_dataset(joinpath(fixedvarsinputdir, "areacello.nc"))
51+
52+
# # Load fixed variables in memory
53+
# areacello = readcubedata(areacello_ds.areacello)
54+
# volcello = readcubedata(volcello_ds.volcello)
55+
# lon = readcubedata(volcello_ds.lon)
56+
# lat = readcubedata(volcello_ds.lat)
57+
# lev = volcello_ds.lev
58+
# # Identify the vertices keys (vary across CMIPs / models)
59+
# volcello_keys = propertynames(volcello_ds)
60+
# lon_vertices_key = volcello_keys[findfirst(x -> occursin("lon", x) & occursin("vert", x), string.(volcello_keys))]
61+
# lat_vertices_key = volcello_keys[findfirst(x -> occursin("lat", x) & occursin("vert", x), string.(volcello_keys))]
62+
# lon_vertices = readcubedata(getproperty(volcello_ds, lon_vertices_key))
63+
# lat_vertices = readcubedata(getproperty(volcello_ds, lat_vertices_key))
64+
# # Make makegridmetrics
65+
# gridmetrics = makegridmetrics(; areacello, volcello, lon, lat, lev, lon_vertices, lat_vertices)
66+
# (; lon_vertices, lat_vertices, lon, lat, zt, v3D, thkcello, Z3D) = gridmetrics
67+
# lev = zt
68+
# # Make indices
69+
# indices = makeindices(gridmetrics.v3D)
70+
# (; wet3D, N) = indices
71+
72+
# experiment_dir = "/scratch/xv83/TMIP/data/$model/$(experiments[1])"
73+
# mlotst_files_1850s = [joinpath(experiment_dir, member, time_window_1850s, "cyclomonth", "mlotst.nc") for member in members]
74+
# mlotst_1850s_ds = open_mfdataset(DimArray(mlotst_files_1850s, Dim{:member}(members)))
75+
# mlotst_1850s = readcubedata(mlotst_1850s_ds.mlotst)
76+
77+
# experiment_dir = "/scratch/xv83/TMIP/data/$model/$(experiments[2])"
78+
# mlotst_files_2030s = [joinpath(experiment_dir, member, time_window_2030s, "cyclomonth", "mlotst.nc") for member in members]
79+
# mlotst_2030s_ds = open_mfdataset(DimArray(mlotst_files_2030s, Dim{:member}(members)))
80+
# mlotst_2030s = readcubedata(mlotst_2030s_ds.mlotst)
81+
82+
# experiment_dir = "/scratch/xv83/TMIP/data/$model/$(experiments[3])"
83+
# mlotst_files_2090s = [joinpath(experiment_dir, member, time_window_2090s, "cyclomonth", "mlotst.nc") for member in members]
84+
# mlotst_2090s_ds = open_mfdataset(DimArray(mlotst_files_2090s, Dim{:member}(members)))
85+
# mlotst_2090s = readcubedata(mlotst_2090s_ds.mlotst)
86+
87+
# mlotst_1850s_yearlymax = dropdims(maximum(mlotst_1850s, dims = :month), dims = :month)
88+
# mlotst_2030s_yearlymax = dropdims(maximum(mlotst_2030s, dims = :month), dims = :month)
89+
# mlotst_2090s_yearlymax = dropdims(maximum(mlotst_2090s, dims = :month), dims = :month)
90+
91+
# mlotst_1850s_yearlymax_ensemblemean = dropdims(mean(mlotst_1850s_yearlymax, dims = :member), dims = :member)
92+
# mlotst_1850s_yearlymax_ensemblemax = dropdims(maximum(mlotst_1850s_yearlymax, dims = :member), dims = :member)
93+
# mlotst_1850s_yearlymax_ensemblemin = dropdims(minimum(mlotst_1850s_yearlymax, dims = :member), dims = :member)
94+
# mlotst_1850s_yearlymax_ensemblerange = mlotst_1850s_yearlymax_ensemblemax - mlotst_1850s_yearlymax_ensemblemin
95+
96+
# mlotst_2030s_yearlymax_ensemblemean = dropdims(mean(mlotst_2030s_yearlymax, dims = :member), dims = :member)
97+
# mlotst_2030s_yearlymax_ensemblemax = dropdims(maximum(mlotst_2030s_yearlymax, dims = :member), dims = :member)
98+
# mlotst_2030s_yearlymax_ensemblemin = dropdims(minimum(mlotst_2030s_yearlymax, dims = :member), dims = :member)
99+
# mlotst_2030s_yearlymax_ensemblerange = mlotst_2030s_yearlymax_ensemblemax - mlotst_2030s_yearlymax_ensemblemin
100+
101+
# mlotst_2090s_yearlymax_ensemblemean = dropdims(mean(mlotst_2090s_yearlymax, dims = :member), dims = :member)
102+
# mlotst_2090s_yearlymax_ensemblemax = dropdims(maximum(mlotst_2090s_yearlymax, dims = :member), dims = :member)
103+
# mlotst_2090s_yearlymax_ensemblemin = dropdims(minimum(mlotst_2090s_yearlymax, dims = :member), dims = :member)
104+
# mlotst_2090s_yearlymax_ensemblerange = mlotst_2090s_yearlymax_ensemblemax - mlotst_2090s_yearlymax_ensemblemin
105+
106+
# MLD_ensemble_means = [mlotst_1850s_yearlymax_ensemblemean, mlotst_2030s_yearlymax_ensemblemean, mlotst_2090s_yearlymax_ensemblemean]
107+
# MLD_ensemble_ranges = [mlotst_1850s_yearlymax_ensemblerange, mlotst_2030s_yearlymax_ensemblerange, mlotst_2090s_yearlymax_ensemblerange]
108+
109+
110+
include("plotting_functions.jl")
111+
112+
usecontourf = false
113+
114+
axs = Array{Any,2}(undef, (3, 2))
115+
contours = Array{Any,2}(undef, (3, 2))
116+
nrows, ncols = size(axs)
117+
118+
fig = Figure(size = (ncols * 500, nrows * 250 + 100), fontsize = 18)
119+
120+
yticks = -60:30:60
121+
xticks = -120:60:120 + 360
122+
123+
# myscale = ReversibleScale(
124+
# x -> sign(x) * log10(abs(x / 5) + 1),
125+
# x -> sign(x) * (exp10(abs(x)) - 1) * 5;
126+
# # x -> x,
127+
# # x -> x;
128+
# limits=(0f0, 3f0),
129+
# name=:myscale
130+
# )
131+
132+
# levels = [0, 50, 70, 100, 140, 200, 300, 500, 700, 1000, 1400, 2000, 3000]
133+
levels = [0, 50, 100, 200, 500, 1000, 2000]
134+
colorscale = mk_piecewise_linear(levels)
135+
136+
colorrange = extrema(levels)
137+
# pseudocolorrange = myscale.(colorrange)
138+
colormap = cgrad(:thermal, length(levels); categorical=true)
139+
extendlow = lowclip = nothing
140+
extendhigh = highclip = colormap[end]
141+
colormap = cgrad(colormap[1:end-1], categorical=true)
142+
143+
# colormap = cgrad(:tol_ylorbr, length(levels); categorical=true)
144+
# lowclip = nothing
145+
# highclip = colormap[end]
146+
# colormap = cgrad(colormap[1:end-1], categorical=true)
147+
# # pseudologlevels = myscale.(levels)
148+
149+
for (irow, (MLD_ensemble_mean, MLD_ensemble_range)) in enumerate(zip(MLD_ensemble_means, MLD_ensemble_ranges))
150+
# Plot ensemble mean
151+
icol = 1
152+
axs[irow, icol] = ax = Axis(fig[irow, icol]; yticks, xticks, xtickformat, ytickformat, aspect = DataAspect())
153+
contours[irow, icol] = if usecontourf
154+
plotcontourfmap!(ax, MLD_ensemble_mean, gridmetrics; levels, colormap, extendhigh, colorscale)
155+
else
156+
plotmap!(ax, MLD_ensemble_mean, gridmetrics; colorrange, colormap, highclip, colorscale) # <- need to fix wrapping longitude for contour levels
157+
end
158+
myhidexdecorations!(ax, irow < nrows)
159+
myhideydecorations!(ax, icol > 1)
160+
161+
# Plot ensemble range
162+
icol = 2
163+
axs[irow, icol] = ax = Axis(fig[irow, icol]; yticks, xticks, xtickformat, ytickformat, aspect = DataAspect())
164+
contours[irow, icol] = if usecontourf
165+
plotcontourfmap!(ax, MLD_ensemble_range, gridmetrics; levels, colormap, colorscale, extendhigh)
166+
else
167+
plotmap!(ax, MLD_ensemble_range, gridmetrics; colorrange, colormap, highclip, colorscale) # <- need to fix wrapping longitude for contour levels
168+
end
169+
myhidexdecorations!(ax, irow < nrows)
170+
myhideydecorations!(ax, icol > 1)
171+
172+
Label(fig[irow, 0]; text = "$(experiments2[irow]) $(time_windows[irow])", rotation = π/2, tellheight = false)
173+
end
174+
175+
176+
label = "ensemble mean MLD (m)"
177+
# cb = Colorbar(fig[nrows + 1, 1], contours[1, 1]; label, vertical = false, flipaxis = false, ticks = levels)
178+
# cb.width = Relative(2/3)
179+
180+
cb = Colorbar(fig[nrows + 1, 1];
181+
limits = (1, length(levels)),
182+
label,
183+
vertical = false,
184+
flipaxis = false,
185+
colormap,
186+
highclip,
187+
ticks = (1:length(levels), string.(levels)),
188+
)
189+
cb.width = Relative(2/3)
190+
191+
label = "ensemble range MLD (m)"
192+
cb = Colorbar(fig[nrows + 1, 2];
193+
limits = (1, length(levels)),
194+
label,
195+
vertical = false,
196+
flipaxis = false,
197+
colormap,
198+
highclip,
199+
ticks = (1:length(levels), string.(levels)),
200+
)
201+
cb.width = Relative(2/3)
202+
# cb = Colorbar(fig[nrows + 1, 2], contours[1, 2]; label, vertical = false, flipaxis = false, ticks = levels, scale = colorscale)
203+
# cb.width = Relative(2/3)
204+
205+
# column labels
206+
# Label(fig[0, 1]; text = "ensemble mean", tellwidth = false)
207+
# Label(fig[0, 2]; text = "ensemble range (internal variability)", tellwidth = false)
208+
209+
labels = [
210+
"a" "b"
211+
"c" "d"
212+
"e" "f"
213+
]
214+
215+
labeloptions = (
216+
font = :bold,
217+
align = (:left, :bottom),
218+
offset = (5, 2),
219+
space = :relative,
220+
fontsize = 24
221+
)
222+
223+
for (ax, label) in zip(axs, labels)
224+
txt = text!(ax, 0, 0; text = "$label", labeloptions..., strokecolor = :white, strokewidth = 3)
225+
translate!(txt, 0, 0, 100)
226+
txt = text!(ax, 0, 0; text = "$label", labeloptions...)
227+
translate!(txt, 0, 0, 100)
228+
end
229+
230+
# Label(fig[0, 1:2]; text = "$(time_window[4:7])s Seafloor Sequestration Efficiency ($(length(members)) members)$(yearly_str2)", fontsize = 24, tellwidth = false)
231+
rowgap!(fig.layout, 10)
232+
colgap!(fig.layout, 10)
233+
234+
colsize!(fig.layout, 1, Aspect(1, 2.0))
235+
colsize!(fig.layout, 2, Aspect(1, 2.0))
236+
237+
# save plot
238+
suffix = usecontourf ? "_ctrf" : ""
239+
240+
resize_to_layout!(fig)
241+
242+
outputdir = "/scratch/xv83/TMIP/data/$model/$(experiments[2])/all_members"
243+
244+
245+
outputfile = joinpath(outputdir, "MLDs.png")
246+
@info "Saving MLD image file:\n $(outputfile)"
247+
save(outputfile, fig)
248+
outputfile = joinpath(outputdir, "MLDs.pdf")
249+
@info "Saving MLD image file:\n $(outputfile)"
250+
save(outputfile, fig)

‎src/plotting_functions.jl‎

Lines changed: 51 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,55 @@ ytickformat(y) = latticklabel.(y)
4040

4141
loninsamewindow(l1, l2) = mod(l1 - l2 + 180, 360) + l2 - 180
4242

43+
"""
44+
To be used as `colorscale` of Makie's `contourf` and `heatmap`.
45+
Picewise linear mapping such that
46+
func.([v1,v2, . . . , vn]) == [0, 1, . . . , n-1].
47+
Outside the range [v1, vn], it's the simple lienar extrapolation.
48+
vs = [v1,v2, . . . , vn] must be strictly increasing.
49+
from https://discourse.julialang.org/t/makie-nonlinear-color-levels-in-colorbar/118056/5
50+
"""
51+
function mk_piecewise_linear(vs)
52+
@assert length(vs) > 1
53+
function is_increasing(ss)
54+
prev = ss[1]
55+
for s in ss[2:end]
56+
(s ≤ prev) && return false
57+
prev = s
58+
end
59+
return true
60+
end
61+
@assert is_increasing(vs)
62+
d1 = vs[2] - vs[1]
63+
d2 = vs[end] - vs[end-1]
64+
un = size(vs, 1) - 1
65+
function piecewise_linear(v)
66+
if v <= vs[1]
67+
(v - vs[1]) / d1
68+
elseif v ≥ vs[end]
69+
(v - vs[end]) / d2 + un
70+
else
71+
i = findfirst(q -> v < q, vs) - 1
72+
d = vs[i + 1] - vs[i]
73+
(v - vs[i]) / d + i - 1
74+
end
75+
end
76+
function its_inverse(u)
77+
if u ≤ 0
78+
u * d1 + vs[1]
79+
elseif u ≥ un
80+
(u - un) * d2 + vs[end]
81+
else
82+
iu = floor(Int, u)
83+
i = iu + 1
84+
d = vs[i + 1] - vs[i]
85+
(u - iu) * d + vs[i]
86+
end
87+
end
88+
return ReversibleScale(piecewise_linear, its_inverse)
89+
end
90+
91+
4392
land1 = GeoMakie.land()
4493
mapwindow = GeometryBasics.Polygon([
4594
Point2f(20, -90),
@@ -54,7 +103,7 @@ land2 = GeometryOps.transform(P -> P + Point2f(360, 0), land1)
54103
land2cut = [LibGEOS.intersection(p, mapwindow) for p in land2]
55104
land2cut = [p for p in land2cut if !LibGEOS.isEmpty(p)]
56105

57-
function plotmap!(ax, x2D, gridmetrics; colorrange, colormap, levels=nothing, highclip = automatic, lowclip = automatic)
106+
function plotmap!(ax, x2D, gridmetrics; colorrange, colormap, levels=nothing, highclip = automatic, lowclip = automatic, colorscale = identity)
58107

59108
# unpack gridmetrics
60109
lonv = gridmetrics.lon_vertices
@@ -71,7 +120,7 @@ function plotmap!(ax, x2D, gridmetrics; colorrange, colormap, levels=nothing, hi
71120
colors_per_point = vcat(fill.(vec(x2D), 4)...)
72121

73122
# create plot
74-
plt = mesh!(ax, quad_points, quad_faces; color = colors_per_point, shading = NoShading, colormap, colorrange, rasterize = 2, highclip, lowclip)
123+
plt = mesh!(ax, quad_points, quad_faces; color = colors_per_point, shading = NoShading, colormap, colorrange, rasterize = 2, highclip, lowclip, colorscale)
75124
xlims!(ax, (20, 20 + 360))
76125
ylims!(ax, (-90, 90))
77126

0 commit comments

Comments
 (0)