tcreate_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
       ---
       tcreate_timeline.py (5580B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 # @package create_timeline
            4 # \author Andy Aschwanden, University of Alaska Fairbanks, USA
            5 # \brief Script creates a timeline for enviornmental forcing.
            6 # \details Script creates a timeline file that can be used with
            7 # the -time_file option and -*_given_file to control the applicaton
            8 # of forcing data, and to control start and end date of a PISM simulation.
            9 #
           10 # Say you have monthly climate forcing from 1980-1-1 through 2001-1-1 in
           11 # the forcing file foo_1980-1999.nc to be used with, e.g. -surface_given_file,
           12 # but you want the model to run from 1991-1-1 through 2001-1-1.
           13 #
           14 # Usage:
           15 #
           16 # \verbatim $ create_timeline.py --start_date '1991-1-1' -end_date '2001-1-1'
           17 # time_1991-2000.nc \endverbatim
           18 
           19 import os
           20 from argparse import ArgumentParser
           21 from dateutil import rrule
           22 from dateutil.parser import parse
           23 from datetime import datetime
           24 import time
           25 import numpy as np
           26 
           27 try:
           28     import netCDF4 as netCDF
           29 except:
           30     print("netCDF4 is not installed!")
           31     sys.exit(1)
           32 NC = netCDF.Dataset
           33 from netcdftime import utime
           34 
           35 # Set up the option parser
           36 parser = ArgumentParser()
           37 parser.description = '''Script creates a time file with time and time
           38 bounds that can be used to determine to force PISM via command line
           39 option -time_file'''
           40 parser.add_argument("FILE", nargs='*')
           41 parser.add_argument("-p", "--periodicity", dest="periodicity",
           42                     help='''periodicity, e.g. monthly, daily, etc. Default=monthly''',
           43                     default="monthly")
           44 parser.add_argument("-a", "--start_date", dest="start_date",
           45                     help='''Start date in ISO format. Default=1989-1-1''',
           46                     default='1989-1-1')
           47 parser.add_argument("-e", "--end_date", dest="end_date",
           48                     help='''End date in ISO format. Default=2012-1-1''',
           49                     default='2012-1-1')
           50 parser.add_argument("-i", "--interval_type", dest="interval_type",
           51                     choices=['start', 'mid', 'end'],
           52                     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".''',
           53                     default='mid')
           54 parser.add_argument("-u", "--ref_unit", dest="ref_unit",
           55                     help='''Reference unit. Default=days. Use of months or
           56                     years is NOT recommended.''', default='days')
           57 parser.add_argument("-d", "--ref_date", dest="ref_date",
           58                     help='''Reference date. Default=1960-1-1''',
           59                     default='1960-1-1')
           60 
           61 options = parser.parse_args()
           62 interval_type = options.interval_type
           63 periodicity = options.periodicity.upper()
           64 start_date = parse(options.start_date)
           65 end_date = parse(options.end_date)
           66 ref_unit = options.ref_unit
           67 ref_date = options.ref_date
           68 args = options.FILE
           69 infile = args[0]
           70 
           71 # Check if file exists. If True, append, otherwise create it.
           72 if os.path.isfile(infile):
           73     nc = NC(infile, 'a')
           74 else:
           75     nc = NC(infile, 'w', format='NETCDF3_CLASSIC')
           76 
           77 time_units = ("%s since %s" % (ref_unit, ref_date))
           78 # currently PISM only supports the gregorian standard calendar
           79 # once this changes, calendar should become a command-line option
           80 time_calendar = "standard"
           81 
           82 cdftime = utime(time_units, time_calendar)
           83 
           84 # create a dictionary so that we can supply the periodicity as a
           85 # command-line argument.
           86 pdict = {}
           87 pdict['SECONDLY'] = rrule.SECONDLY
           88 pdict['MINUTELY'] = rrule.MINUTELY
           89 pdict['HOURLY'] = rrule.HOURLY
           90 pdict['DAILY'] = rrule.DAILY
           91 pdict['WEEKLY'] = rrule.WEEKLY
           92 pdict['MONTHLY'] = rrule.MONTHLY
           93 pdict['YEARLY'] = rrule.YEARLY
           94 prule = pdict[periodicity]
           95 
           96 # reference date from command-line argument
           97 r = time_units.split(' ')[2].split('-')
           98 refdate = datetime(int(r[0]), int(r[1]), int(r[2]))
           99 
          100 # create list with dates from start_date until end_date with
          101 # periodicity prule.
          102 bnds_datelist = list(rrule.rrule(prule, dtstart=start_date, until=end_date))
          103 
          104 # calculate the days since refdate, including refdate, with time being the
          105 bnds_interval_since_refdate = cdftime.date2num(bnds_datelist)
          106 if interval_type == 'mid':
          107     # mid-point value:
          108     # time[n] = (bnds[n] + bnds[n+1]) / 2
          109     time_interval_since_refdate = (bnds_interval_since_refdate[0:-1] +
          110                                    np.diff(bnds_interval_since_refdate) / 2)
          111 elif interval_type == 'start':
          112     time_interval_since_refdate = bnds_interval_since_refdate[:-1]
          113 else:
          114     time_interval_since_refdate = bnds_interval_since_refdate[1:]
          115 
          116 # create a new dimension for bounds only if it does not yet exist
          117 time_dim = "time"
          118 if time_dim not in list(nc.dimensions.keys()):
          119     nc.createDimension(time_dim)
          120 
          121 # create a new dimension for bounds only if it does not yet exist
          122 bnds_dim = "nb2"
          123 if bnds_dim not in list(nc.dimensions.keys()):
          124     nc.createDimension(bnds_dim, 2)
          125 
          126 # variable names consistent with PISM
          127 time_var_name = "time"
          128 bnds_var_name = "time_bnds"
          129 
          130 # create time variable
          131 time_var = nc.createVariable(time_var_name, 'd', dimensions=(time_dim))
          132 time_var[:] = time_interval_since_refdate
          133 time_var.bounds = bnds_var_name
          134 time_var.units = time_units
          135 time_var.calendar = time_calendar
          136 time_var.standard_name = time_var_name
          137 time_var.axis = "T"
          138 
          139 # create time bounds variable
          140 time_bnds_var = nc.createVariable(bnds_var_name, 'd', dimensions=(time_dim, bnds_dim))
          141 time_bnds_var[:, 0] = bnds_interval_since_refdate[0:-1]
          142 time_bnds_var[:, 1] = bnds_interval_since_refdate[1::]
          143 
          144 # writing global attributes
          145 script_command = ' '.join([time.ctime(), ':', __file__.split('/')[-1],
          146                            ' '.join([str(x) for x in args])])
          147 nc.history = script_command
          148 nc.Conventions = "CF 1.5"
          149 nc.close()