Commit a2f0cce6 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Implemented the longterm flag.

parent ca4f8b40
......@@ -71,6 +71,12 @@ def main():
base_parser.add_argument('-o', '--output', metavar='<path/to/output>',
type=str, help='where to store the outputs. if '
f'not given, no output is stored')
base_parser.add_argument('--longterm', action='store_true',
help=f'this flag is intended for longterm models.'
f' if given, all times will be interpreted as '
f'ages before the year 1950 in kyrs., i.e. the '
f'epoch 10 refers to the year -8050 without the'
f'flag')
# parser for arguments that are shared between time series commands
series_base_parser = argparse.ArgumentParser(add_help=False)
......@@ -147,21 +153,29 @@ def main():
default=res_default)
# the model to be used
parser.add_argument('model', type=str, help='a model')
parser.add_argument('-i', '--input', type=str, help=f'input path to a '
f'model file. may be used to evaluate your own model. '
f'if given, this path will be used instead over the '
f'included files')
parser.add_argument('-i', '--input', metavar='<path/to/input>', type=str,
help=f'input path to a model file. may be used to '
f'evaluate your own model. if given, this path will '
f'be used instead over the included files')
# parse the command line arguments
args = parser.parse_args()
# check whether model and input are consistent
if args.input is None and args.model in models:
args.input = models[args.model]
if args.input is None and args.model not in models:
raise parser.error(f'Model {args.model} is not include in pycals. '
f'To process your own model, specify an input file '
f'using \n'
raise parser.error(f'Model {args.model} is not include in pymagglobal.'
f' To process your own model, specify an input '
f'file using \n'
f' $ pycals --input <path/to/{args.model}> '
f'... {args.model}')
# get the right time units
if args.longterm:
args.t_unit = 'kyrs.'
args.t_label = r'age [kyrs.]'
else:
args.t_unit = 'yrs.'
args.t_label = r't [yrs.]'
# call the respective function, specified using .set_defaults(func=...)
args.func(args)
......
......@@ -49,18 +49,24 @@ def valid_epoch(args):
----------
args : object
the object returned by ArgumentParser.parse_args()
Returns
-------
float
the epoch in yrs
'''
# get the range of the model from the model file
with open(args.input, 'r') as fh:
t_min, t_max = map(float, fh.readline().split()[0:2])
if args.epoch < t_min or t_max < args.epoch:
warn(f'Epoch {args.epoch} is outside of the models range [{t_min}, '
if args.longterm:
epoch = lt2yr(args.epoch)
else:
epoch = args.epoch
if epoch < t_min or t_max < epoch:
warn(f'Epoch {epoch} is outside of the models range [{t_min}, '
f'{t_max}]. The result will be extrapolated and not robust.',
category=OutOfRangeWarning)
return False
else:
return True
return epoch
def args2times(args):
......@@ -83,6 +89,17 @@ def args2times(args):
return times
def yr2lt(times):
''' translate times given as years to the format / age often used
for longterm models '''
return -(times - 1950) / 1000
def lt2yr(times):
''' translate times given in the longterm format to years '''
return -times*1000 + 1950
def master_curve(args):
''' handle the master command and create a master curve at a given
location.
......@@ -100,9 +117,13 @@ def master_curve(args):
curves = pymagglobal.master_curve(times, (args.lat, args.lon), splines,
ser_type=args.type)
# output formats for dif and nez components
fmts = {'dif': ('%.1f', '%2.6f', '%2.6f', '%1.7e'),
'nez': ('%.1f', '%1.7e', '%1.7e', '%1.7e')}
fmts = {'dif': ('%.2f', '%2.6f', '%2.6f', '%1.7e'),
'nez': ('%.2f', '%1.7e', '%1.7e', '%1.7e')}
# if the --longterm flag is set, translate the years accordingly
if args.longterm:
times = yr2lt(times)
# if an output is specified, save the result to text
if args.output is not None:
np.savetxt(args.output,
......@@ -114,7 +135,7 @@ def master_curve(args):
delimiter=',',
header=f'Master curves at ({args.lat}°, {args.lon}°) '
f'for the model {args.model}\n'
f't [yrs.],'
f'{args.t_label}'
f'{utils.labels[args.type][0]},'
f'{utils.labels[args.type][1]},'
f'{utils.labels[args.type][2]}')
......@@ -132,8 +153,12 @@ def master_curve(args):
for it in range(3):
axs[it].plot(times, curves[it])
axs[it].set_title(utils.names[args.type][it])
axs[it].set_xlabel(r'$t$ [yrs.]')
axs[it].set_ylabel(utils.nicelabel(utils.labels[args.type][it]))
if args.longterm:
axs[it].invert_xaxis()
axs[it].set_xlabel(args.t_label)
else:
axs[it].set_xlabel(utils.nicelabel(args.t_label))
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
......@@ -153,17 +178,20 @@ def dipole_series(args):
splines = pymagglobal.file2splines(args.input)
# create the dipole-moment time series for the given model
dip_ser = pymagglobal.dipole_series(times, splines)
# if the --longterm flag is set, translate the years accordingly
if args.longterm:
times = yr2lt(times)
# if an output is specified, save the result to text
if args.output is not None:
np.savetxt(args.output,
np.array([times,
dip_ser]).T,
fmt=('%.1f', '%1.7e'),
fmt=('%.2f', '%1.7e'),
delimiter=',',
header=f'Dipole moment series for for the model '
f'{args.model}\n'
f't [yrs.],m [m^2 A]')
f'{args.t_label},m [m^2 A]')
# if the --no-show flag is not set, plot the time series
if not args.no_show:
......@@ -171,8 +199,12 @@ def dipole_series(args):
ax = fig.add_subplot(111)
ax.plot(times, dip_ser)
ax.set_title(f'Dipole moment series for for the model {args.model}')
ax.set_xlabel(r'$t$ [yrs.]')
ax.set_ylabel(r'$m$ [m$^2$ A]')
if args.longterm:
ax.invert_xaxis()
ax.set_xlabel(args.t_label)
else:
ax.set_xlabel(utils.nicelabel(args.t_label))
fig.tight_layout()
plt.show()
......@@ -187,12 +219,11 @@ def coefficients(args):
the object returned by ArgumentParser.parse_args()
'''
# check whether the epoch is in the range of the model
valid_epoch(args)
epoch = valid_epoch(args)
# get the splines from the input file
splines = pymagglobal.file2splines(args.input)
# evaluate the splines at the epoch
ls, ms, coeffs = pymagglobal.coefficients(args.epoch, splines)
ls, ms, coeffs = pymagglobal.coefficients(epoch, splines)
# if an output is specified, save the result to text
if args.output is not None:
np.savetxt(args.output,
......@@ -200,13 +231,13 @@ def coefficients(args):
fmt=('%.d', '%d', '%1.7e'),
delimiter=',',
header=f'Coefficients for the model {args.model} at epoch '
f'{args.epoch}\n'
f'{args.epoch} {args.t_unit}\n'
f'l,m,g_l^m')
# if no output is given and the --no-show flag is not set, print the
# coefficients to the command line
elif not args.no_show:
print(f'Coefficients for the model {args.model} at epoch '
f'{args.epoch}')
f'{args.epoch} {args.t_unit}')
right = len(str(ms[-1]))
print(f'l\t' + ' '*(right-1) + 'm\t g_l^m')
for l, m, g in zip(ls, ms, coeffs):
......@@ -223,7 +254,7 @@ def maps(args):
the object returned by ArgumentParser.parse_args()
'''
# check whether the epoch is in the range of the model
valid_epoch(args)
epoch = valid_epoch(args)
# set up a grid of equidistant points on the sphere
tps = utils.equi_sph(args.res)
# set upt the points as expected by pyfield
......@@ -232,7 +263,7 @@ def maps(args):
# BUGFIX: white edges in cartopy
z_at[1] = np.rad2deg(tps[1]) - 180
z_at[2] = REARTH
z_at[3] = args.epoch
z_at[3] = epoch
# get the splines from the input file
splines = pymagglobal.file2splines(args.input)
# use the core function to get the field
......@@ -255,7 +286,7 @@ def maps(args):
fmt=fmts[args.type],
delimiter=',',
header=f'Field map for the model {args.model} at epoch '
f'{args.epoch}\n'
f'{args.epoch} {args.t_unit}\n'
f'lat,lon,'
f'{utils.labels[args.type][0]},'
f'{utils.labels[args.type][1]},'
......@@ -288,7 +319,7 @@ def maps(args):
fig, axs = plt.subplots(1, 3, figsize=(13, 3),
subplot_kw={'projection': proj})
fig.suptitle(f'Field maps for the model {args.model} at epoch '
f'{args.epoch}')
f'{args.epoch} {args.t_unit}')
for it in range(3):
axs[it].tripcolor(plt_lat, plt_lon, field[it],
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment