commit 1c7dcced9c0a9c9b8915593268159a4d4490e322
parent 92cec776e88d9a786c967f6d5e5f3bc0faf8d246
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Tue, 7 Dec 2021 22:38:28 +0100
add plots for AGU
Diffstat:
3 files changed, 143 insertions(+), 9 deletions(-)
diff --git a/Makefile b/Makefile
@@ -30,14 +30,22 @@ OUT_COMMON =\
-ys 0 -ye ${T_END}\
all: \
+ ex_deltaSL-linear_1d-till-evol.pdf\
ex_deltaSL-linear_1d-profile.pdf\
+ ex_constant-linear_1d-till-evol.pdf\
ex_constant-linear_1d-profile.pdf\
ex_nt-deltaSL-linear_1d-profile.pdf\
ex_nt-constant-linear_1d-profile.pdf\
+ex_deltaSL-linear_1d-till-evol.pdf: ex_deltaSL-linear_1d.nc plot-till-evolution.py
+ ./plot-till-evolution.py ex_deltaSL-linear_1d.nc
+
ex_deltaSL-linear_1d-profile.pdf: ex_deltaSL-linear_1d.nc plot-time-series.py
./plot-time-series.py ex_deltaSL-linear_1d.nc
+ex_constant-linear_1d-till-evol.pdf: ex_constant-linear_1d.nc plot-till-evolution.py
+ ./plot-till-evolution.py ex_constant-linear_1d.nc
+
ex_constant-linear_1d-profile.pdf: ex_constant-linear_1d.nc plot-time-series.py
./plot-time-series.py ex_constant-linear_1d.nc
diff --git a/plot-till-evolution.py b/plot-till-evolution.py
@@ -0,0 +1,123 @@
+#!/usr/bin/env python3
+
+from pylab import figure, subplots, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout, cm, legend
+from sys import exit
+
+import MISMIP
+
+import numpy as np
+from optparse import OptionParser
+import os.path
+
+try:
+ from netCDF4 import Dataset as NC
+except:
+ print("netCDF4 is not installed!")
+ sys.exit(1)
+
+def process_options():
+ "Process command-line options and arguments."
+ parser = OptionParser()
+ parser.usage = "%prog <input files> [options]"
+ parser.description = "Plots the ice flux as a function of the distance from the divide."
+ parser.add_option("-o", "--output", dest="output", type="string",
+ help="Output image file name (e.g. -o foo.png)")
+
+ opts, args = parser.parse_args()
+
+ if len(args) == 0:
+ print("ERROR: An input file is requied.")
+ exit(0)
+
+ if len(args) > 1 and opts.output:
+ print("More than one input file given. Ignoring the -o option...\n")
+ opts.output = None
+
+ return args, opts.output, opts
+
+
+def read(filename, name):
+ "Read a variable and extract the middle row."
+ nc = NC(filename)
+
+ try:
+ var = nc.variables[name][:]
+ except:
+ print("ERROR: Variable '%s' not present in '%s'" % (name, filename))
+ exit(1)
+
+ return var
+
+
+def find_grounding_line(x, topg, thk, mask):
+ "Find the modeled grounding line position."
+ # "positive" parts of x, topg, thk, mask
+ topg = topg[x > 0]
+ thk = thk[x > 0]
+ mask = mask[x > 0]
+ x = x[x > 0] # this should go last
+
+ def f(j):
+ "See equation (7) in Pattyn et al, 'Role of transition zones in marine ice sheet dynamics', 2005."
+ z_sl = 0
+ return (z_sl - topg[j]) * MISMIP.rho_w() / (MISMIP.rho_i() * thk[j])
+
+ for j in range(x.size):
+ if mask[j] == 2 and mask[j + 1] == 3: # j is grounded, j+1 floating
+ nabla_f = (f(j + 1) - f(j)) / (x[j + 1] - x[j])
+
+ # See equation (8) in Pattyn et al
+ return (1.0 - f(j) + nabla_f * x[j]) / nabla_f
+
+ raise Exception("Can't find the grounding line")
+
+
+def plot_profile(in_file, out_file):
+
+ if out_file is None:
+ out_file = os.path.splitext(in_file)[0] + "-till-evol.pdf"
+
+ steps = read(in_file, 'thk').shape[0]
+
+ fig, ax = subplots(1, 1, sharex=True, figsize=[6, 3])
+ for i in range(1, steps):
+
+ mask = read(in_file, 'mask')[i]
+ uvelbase = read(in_file, 'uvelbase')[i]
+ utillflux = read(in_file, 'utillflux')[i]
+ till_deposit = read(in_file, 'tilldeposit')[i]
+ x = read(in_file, 'x')
+
+ # convert x to kilometers
+ x /= 1e3
+
+ # modeled grounding line position
+ #xg_PISM = find_grounding_line(x, lsrf, thk, mask)
+ #plot(x, np.zeros_like(x), ls='dotted', color='red')
+ icecolor = cm.cool(i / steps)
+ #ax[0].plot(x, uvelbase, color=icecolor)
+ #ax[1].plot(x, utillflux, color=icecolor) #, label='{}'.format(i))
+ ax.plot(x, till_deposit, color=icecolor)
+
+ ax.set_xlabel('distance from the divide, km')
+ #ax[0].set_ylabel('$v_{SSA}$, m/a')
+ #ax[1].set_ylabel('$q_{t,x}$, m$^2$/a')
+ ax.set_ylabel('$\Delta b$, m')
+
+ #_, _, ymin, ymax = axis(xmin=0, xmax=x.max())
+ _, _, ymin, ymax = axis(xmin=950, xmax=1150)
+ #_, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max())
+
+ #vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')
+ #vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red')
+
+ #legend()
+ fig.tight_layout()
+ print("Saving '%s'...\n" % out_file)
+ savefig(out_file)
+
+if __name__ == "__main__":
+ args, out_file, opts = process_options()
+
+ for in_file in args:
+ plot_profile(in_file, out_file)+
\ No newline at end of file
diff --git a/plot-time-series.py b/plot-time-series.py
@@ -1,6 +1,6 @@
#!/usr/bin/env python3
-from pylab import figure, subplot, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout, cm, legend
+from pylab import figure, subplots, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout, cm, legend
from sys import exit
import MISMIP
@@ -79,8 +79,7 @@ def plot_profile(in_file, out_file):
steps = read(in_file, 'thk').shape[0]
- figure(1)
- ax = subplot(111)
+ fig, ax = subplots(1, 1, figsize=[4, 3])
for i in range(1, steps):
mask = read(in_file, 'mask')[i]
@@ -89,6 +88,7 @@ def plot_profile(in_file, out_file):
lsrf = topg.copy()
thk = read(in_file, 'thk')[i]
sea_lvl = read(in_file, 'sea_level')[i]
+ #till_deposit = read(in_file, 'tilldeposit')[i]
x = read(in_file, 'x')
# mask out ice-free areas
@@ -104,20 +104,22 @@ def plot_profile(in_file, out_file):
#xg_PISM = find_grounding_line(x, lsrf, thk, mask)
#plot(x, np.zeros_like(x), ls='dotted', color='red')
icecolor = cm.cool(i / steps)
- plot(x, topg, color='black')
- plot(x, usurf, color=icecolor, label='{}'.format(i))
- plot(x, lsrf, color=icecolor)
+ ax.plot(x, topg, color='black')
+ ax.plot(x, usurf, color=icecolor, label='{}'.format(i))
+ ax.plot(x, lsrf, color=icecolor)
- xlabel('distance from the divide, km')
- ylabel('elevation, m')
+ ax.set_xlabel('distance from the divide, km')
+ ax.set_ylabel('elevation, m')
#_, _, ymin, ymax = axis(xmin=0, xmax=x.max())
+ _, _, ymin, ymax = axis(xmin=950, xmax=1150, ymin=-500, ymax=1150)
#_, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max())
#vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')
#vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red')
- legend()
+ #legend()
+ fig.tight_layout()
print("Saving '%s'...\n" % out_file)
savefig(out_file)
