tadjust_timeline.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tadjust_timeline.py (5374B)
---
1 #!/usr/bin/env python3
2
3 # @package adjust_timeline
4 # \author Andy Aschwanden, University of Alaska Fairbanks, USA
5 # \brief Script adjusts a time axis of a file.
6 # \details Script adjusts the time axis of a file.
7
8 # Say you have monthly climate forcing from 1980-1-1 through 2001-1-1 in
9 # the forcing file foo_1980-1999.nc to be used with, e.g. -surface_given_file,
10 # but you want the model to run from 1991-1-1 through 2001-1-1.
11 #
12 # Usage:
13 #
14 # \verbatim $ adjust_timeline.py --start_date '1991-1-1'
15 # time_1991-2000.nc \endverbatim
16
17 import os
18 from argparse import ArgumentParser
19 from dateutil import rrule
20 from dateutil.parser import parse
21 import time
22 import numpy as np
23
24 try:
25 import netCDF4 as netCDF
26 except:
27 print("netCDF4 is not installed!")
28 sys.exit(1)
29 NC = netCDF.Dataset
30 from netcdftime import utime, datetime
31
32 # Set up the option parser
33 parser = ArgumentParser()
34 parser.description = """Script adjusts the time file with time and time
35 bounds that can be used to determine to force PISM via command line
36 option -time_file or adjust the time axis for postprocessing."""
37 parser.add_argument("FILE", nargs="*")
38 parser.add_argument(
39 "-p",
40 "--periodicity",
41 dest="periodicity",
42 help="""periodicity, e.g. monthly, daily, etc. Default=monthly""",
43 default="monthly",
44 )
45 parser.add_argument(
46 "-a",
47 "--start_date",
48 dest="start_date",
49 help="""Start date in ISO format. Default=1989-1-1""",
50 default="1989-1-1",
51 )
52 parser.add_argument(
53 "-c",
54 "--calendar",
55 dest="calendar",
56 choices=["standard", "gregorian", "no_leap", "365_day", "360_day", "julian"],
57 help="""Sets the calendar. Default="standard".""",
58 default="standard",
59 )
60 parser.add_argument(
61 "-i",
62 "--interval_type",
63 dest="interval_type",
64 choices=["start", "mid", "end"],
65 help="""Defines whether the time values t_k are the end points of the time bounds tb_k or the mid points 1/2*(tb_k -tb_(k-1)). Default="mid".""",
66 default="mid",
67 )
68 parser.add_argument(
69 "-u",
70 "--ref_unit",
71 dest="ref_unit",
72 help="""Reference unit. Default=days. Use of months or
73 years is NOT recommended.""",
74 default="days",
75 )
76 parser.add_argument(
77 "-d",
78 "--ref_date",
79 dest="ref_date",
80 help="""Reference date. Default=1960-1-1""",
81 default="1960-1-1",
82 )
83
84 options = parser.parse_args()
85 interval_type = options.interval_type
86 periodicity = options.periodicity.upper()
87 start_date = parse(options.start_date)
88 ref_unit = options.ref_unit
89 ref_date = options.ref_date
90 args = options.FILE
91 infile = args[0]
92
93 time1 = time.time()
94
95
96 nc = NC(infile, "a")
97 nt = len(nc.variables["time"])
98
99 time_units = "%s since %s" % (ref_unit, ref_date)
100 time_calendar = options.calendar
101
102 cdftime = utime(time_units, time_calendar)
103
104 # create a dictionary so that we can supply the periodicity as a
105 # command-line argument.
106 pdict = {}
107 pdict["SECONDLY"] = rrule.SECONDLY
108 pdict["MINUTELY"] = rrule.MINUTELY
109 pdict["HOURLY"] = rrule.HOURLY
110 pdict["DAILY"] = rrule.DAILY
111 pdict["WEEKLY"] = rrule.WEEKLY
112 pdict["MONTHLY"] = rrule.MONTHLY
113 pdict["YEARLY"] = rrule.YEARLY
114 prule = pdict[periodicity]
115
116 # reference date from command-line argument
117 r = time_units.split(" ")[2].split("-")
118 refdate = datetime(int(r[0]), int(r[1]), int(r[2]))
119
120 # create list with dates from start_date for nt counts
121 # periodicity prule.
122 bnds_datelist = list(rrule.rrule(prule, dtstart=start_date, count=nt + 1))
123
124 # calculate the days since refdate, including refdate, with time being the
125 bnds_interval_since_refdate = cdftime.date2num(bnds_datelist)
126 if interval_type == "mid":
127 # mid-point value:
128 # time[n] = (bnds[n] + bnds[n+1]) / 2
129 time_interval_since_refdate = (
130 bnds_interval_since_refdate[0:-1] + np.diff(bnds_interval_since_refdate) / 2
131 )
132 elif interval_type == "start":
133 time_interval_since_refdate = bnds_interval_since_refdate[:-1]
134 else:
135 time_interval_since_refdate = bnds_interval_since_refdate[1:]
136
137 # create a new dimension for bounds only if it does not yet exist
138 time_dim = "time"
139 if time_dim not in list(nc.dimensions.keys()):
140 nc.createDimension(time_dim)
141
142 # create a new dimension for bounds only if it does not yet exist
143 bnds_dim = "nb2"
144 if bnds_dim not in list(nc.dimensions.keys()):
145 nc.createDimension(bnds_dim, 2)
146
147 # variable names consistent with PISM
148 time_var_name = "time"
149 bnds_var_name = "time_bnds"
150
151 # create time variable
152 if time_var_name not in nc.variables:
153 time_var = nc.createVariable(time_var_name, "d", dimensions=(time_dim))
154 else:
155 time_var = nc.variables[time_var_name]
156 time_var[:] = time_interval_since_refdate
157 time_var.bounds = bnds_var_name
158 time_var.units = time_units
159 time_var.calendar = time_calendar
160 time_var.standard_name = time_var_name
161 time_var.axis = "T"
162
163 # create time bounds variable
164 if bnds_var_name not in nc.variables:
165 time_bnds_var = nc.createVariable(
166 bnds_var_name, "d", dimensions=(time_dim, bnds_dim)
167 )
168 else:
169 time_bnds_var = nc.variables[bnds_var_name]
170 time_bnds_var[:, 0] = bnds_interval_since_refdate[0:-1]
171 time_bnds_var[:, 1] = bnds_interval_since_refdate[1::]
172
173 # writing global attributes
174 script_command = " ".join(
175 [time.ctime(), ":", __file__.split("/")[-1], " ".join([str(x) for x in args])]
176 )
177 nc.history = script_command
178 nc.Conventions = "CF 1.5"
179 nc.close()
180
181 time2 = time.time()
182 print("adjust_timeline.py took {:2.1f}s".format(time2 - time1))