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))