jiwoo branch ready?

share/test_data_files.txt

@@ -1,3 +1,4 @@
 https://cdat.llnl.gov/cdat/pmp
-9c55c78adc9cbfba4ad98b1e68c403aa monsoon_sperber_aus.nc
 60582c6986451ffd4925a303c4f0040a sample_data_pr_CMCC.nc
+37acc58e3ba0acd6b35ad3439ee96453 test_monsoon_sperber_input.nc
+3731b19e00226bb6a482c27e19bc1d91 test_monsoon_sperber_input_lf.nc

src/python/monsoon_sperber/doc/myParam.py

@@ -3,7 +3,6 @@
 # Background Information
 # -------------------------------------------------
 mip = 'cmip5'
-#exp = 'piControl'
 exp = 'historical'
 frequency = 'da'
 realm = 'atm'
@@ -12,11 +11,11 @@
 # Observation
 # -------------------------------------------------
 reference_data_name = 'GPCP'
-reference_data_path = '/p/user_pub/pmp/pmp_results/tree_v0.3/pmp_v1.1.2/data/PMPObs/PMPObs_v1.3/atmos/day/pr/GPCP-1-3/gn/v20180816/pr_day_GPCP-1-3_BE_gn_19961002-20170101.nc'
+reference_data_path = '/p/user_pub/pmp/pmp_results/tree_v0.3/pmp_v1.1.2/data/PMPObs/PMPObs_v1.3/atmos/day/pr/GPCP-1-3/gn/v20180816/pr_day_GPCP-1-3_BE_gn_19961002-20170101.nc' # noqa
 reference_data_lf_path = '/work/lee1043/DATA/LandSeaMask_1x1_NCL/NCL_LandSeaMask_rewritten.nc'
 
 varOBS = 'pr'
-ObsUnitsAdjust = (True, 'multiply', 86400.0) # kg m-2 s-1 to mm day-1
+ObsUnitsAdjust = (True, 'multiply', 86400.0)  # kg m-2 s-1 to mm day-1
 
 osyear = 1996
 oeyear = 2016
@@ -29,15 +28,12 @@
 modpath = '/work/lee1043/ESGF/xmls/cmip5/historical/day/pr/cmip5.%(model).%(exp).%(realization).day.pr.xml'
 modpath_lf = '/work/lee1043/ESGF/xmls/cmip5/fx/fx/sftlf/cmip5.%(model).fx.r0i0p0.fx.sftlf.xml'
 
-modnames = ['ACCESS1-0', 'ACCESS1-3', 'BCC-CSM1-1', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-BGC', 'CESM1-CAM5', 'CESM1-FASTCHEM', 'CMCC-CESM', 'CMCC-CM', 'CMCC-CMS', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'EC-EARTH', 'FGOALS-g2', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadGEM2-AO', 'HadGEM2-CC', 'HadGEM2-ES', 'INMCM4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC-ESM', 'MIROC-ESM-CHEM', 'MIROC4h', 'MIROC5', 'MPI-ESM-MR', 'MPI-ESM-P', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M']
+modnames = ['ACCESS1-0', 'ACCESS1-3', 'BCC-CSM1-1', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-BGC', 'CESM1-CAM5', 'CESM1-FASTCHEM', 'CMCC-CESM', 'CMCC-CM', 'CMCC-CMS', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'EC-EARTH', 'FGOALS-g2', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadGEM2-AO', 'HadGEM2-CC', 'HadGEM2-ES', 'INMCM4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC-ESM', 'MIROC-ESM-CHEM', 'MIROC4h', 'MIROC5', 'MPI-ESM-MR', 'MPI-ESM-P', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M'] # noqa
 
-#modnames = ['CanESM2']
-
-#realization = '*' # realizations
 realization = 'r1i1p1'
 
 varModel = 'pr'
-ModUnitsAdjust = (True, 'multiply', 86400.0) # kg m-2 s-1 to mm day-1
+ModUnitsAdjust = (True, 'multiply', 86400.0)  # kg m-2 s-1 to mm day-1
 units = 'mm/d'
 
 msyear = 1961
@@ -53,6 +49,5 @@
 # =================================================
 # Miscellaneous
 # -------------------------------------------------
-update_json = True # False
-#debug = True # False
-debug = False # False
+update_json = True
+debug = False

src/python/monsoon_sperber/doc/plot_monsoon_domain.py

@@ -2,43 +2,46 @@
 import matplotlib.pyplot as plt
 from matplotlib.patches import Polygon
 
-def draw_screen_poly( lats, lons, m):
- x, y = m( lons, lats )
- xy = zip(x,y)
- poly = Polygon( xy, facecolor='red', alpha=0.4, edgecolor='red', linewidth=2)
+
+def draw_screen_poly(lats, lons, m):
+ x, y = m(lons, lats)
+ xy = zip(x, y)
+ poly = Polygon(xy, facecolor='red', alpha=0.4,
+ edgecolor='red', linewidth=2)
  plt.gca().add_patch(poly)
 
+
 lats = {}
 lons = {}
 
-lats['AIR'] = [ 7, 25, 25, 7 ]
-lons['AIR'] = [ 65, 65, 85, 85 ]
+lats['AIR'] = [7, 25, 25, 7]
+lons['AIR'] = [65, 65, 85, 85]
 
-lats['AUS'] = [ -20, -10, -10, -20 ]
-lons['AUS'] = [ 120, 120, 150, 150 ]
+lats['AUS'] = [-20, -10, -10, -20]
+lons['AUS'] = [120, 120, 150, 150]
 
-lats['Sahel'] = [ 13, 18, 18, 13 ]
-lons['Sahel'] = [ -10, -10, 10, 10 ]
+lats['Sahel'] = [13, 18, 18, 13]
+lons['Sahel'] = [-10, -10, 10, 10]
 
-lats['GoG'] = [ 0, 5, 5, 0 ]
-lons['GoG'] = [ -10, -10, 10, 10 ]
+lats['GoG'] = [0, 5, 5, 0]
+lons['GoG'] = [-10, -10, 10, 10]
 
-lats['NAM'] = [ 20, 37, 37, 20 ]
-lons['NAM'] = [ -112, -112, -103, -103 ]
+lats['NAM'] = [20, 37, 37, 20]
+lons['NAM'] = [-112, -112, -103, -103]
 
-lats['SAM'] = [ -20, -2.5, -2.5, -20 ]
-lons['SAM'] = [ -65, -65, -40, -40 ]
+lats['SAM'] = [-20, -2.5, -2.5, -20]
+lons['SAM'] = [-65, -65, -40, -40]
 
 regions = lats.keys()
 
-plt.figure(figsize=(12,6))
+plt.figure(figsize=(12, 6))
 
 m = Basemap(lon_0=0)
 m.drawmapboundary(fill_color='aqua')
-m.fillcontinents(color='grey',lake_color='aqua')
+m.fillcontinents(color='grey', lake_color='aqua')
 m.drawcoastlines()
 
 for region in regions:
- draw_screen_poly( lats[region], lons[region], m )
+ draw_screen_poly(lats[region], lons[region], m)
 
 plt.savefig('monsoon_domain_map.png')

src/python/monsoon_sperber/lib/__init__.py

@@ -1,4 +1,4 @@
 from .argparse_functions import AddParserArgument, YearCheck # noqa
 from .calc_metrics import sperber_metrics # noqa
 from .model_land_only import model_land_only # noqa
-from .divide_chunks import divide_chunks, divide_chunks, divide_chunks_advanced, interp1d #noqa
+from .divide_chunks import divide_chunks, divide_chunks_advanced, interp1d # noqa

src/python/monsoon_sperber/lib/argparse_functions.py

@@ -17,21 +17,28 @@ def AddParserArgument(P):
  P.add_argument("--modpath_lf",
  type=str,
  help="Path of landsea mask for model data set")
- P.add_argument("--varobs", dest="varOBS", type=str, help="Variable name in reference data set")
- P.add_argument("--varmod", dest="varMOD", type=str, help="Variable name in model data set")
- P.add_argument("--obs_units_adjust", dest="ObsUnitsAdjust", type=tuple, 
+ P.add_argument("--varobs", dest="varOBS", type=str,
+ help="Variable name in reference data set")
+ P.add_argument("--varmod", dest="varMOD", type=str,
+ help="Variable name in model data set")
+ P.add_argument("--obs_units_adjust", dest="ObsUnitsAdjust", type=tuple,
  help="Unit conversion\n"
  "- if needed: e.g., (True, 'multiply', 86400.),\n"
  "- no needed: (False, 0, 0)")
- P.add_argument("--mod_units_adjust", dest="ModUnitsAdjust", type=tuple, 
+ P.add_argument("--mod_units_adjust", dest="ModUnitsAdjust", type=tuple,
  help="Unit conversion\n"
  "- if needed: e.g., (True, 'multiply', 86400.),\n"
  "- no needed: (False, 0, 0)")
- P.add_argument("--units", dest="units", type=str, help="Final units for the variable")
- P.add_argument("--osyear", dest="osyear", type=int, help="Start year for reference data set")
- P.add_argument("--msyear", dest="msyear", type=int, help="Start year for model data set")
- P.add_argument("--oeyear", dest="oeyear", type=int, help="End year for reference data set")
- P.add_argument("--meyear", dest="meyear", type=int, help="End year for model data set")
+ P.add_argument("--units", dest="units", type=str,
+ help="Final units for the variable")
+ P.add_argument("--osyear", dest="osyear", type=int,
+ help="Start year for reference data set")
+ P.add_argument("--msyear", dest="msyear", type=int,
+ help="Start year for model data set")
+ P.add_argument("--oeyear", dest="oeyear", type=int,
+ help="End year for reference data set")
+ P.add_argument("--meyear", dest="meyear", type=int,
+ help="End year for model data set")
  P.add_argument("--modnames",
  type=list,
  default=None,
@@ -48,17 +55,17 @@ def AddParserArgument(P):
  type=bool,
  default=False,
  help="Option for debug: False (defualt) or True")
- P.add_argument("--nc_out", dest="nc_out", help="record netcdf output", action="store_true", default=False)
- P.add_argument("--plot", dest="plot", help="produce plots", action="store_true", default=False)
- P.add_argument("--include_obs", dest="includeOBS", help="include observation", action="store_true", default=False)
+ P.add_argument("--nc_out", dest="nc_out",
+ help="record netcdf output", action="store_true", default=False)
+ P.add_argument("--plot", dest="plot", help="produce plots",
+ action="store_true", default=False)
+ P.add_argument("--include_obs", dest="includeOBS",
+ help="include observation", action="store_true", default=False)
  return P
 
 
 def YearCheck(syear, eyear, P):
  if syear >= eyear:
- P.error('Given starting year, '
- + str(syear)
- + ', is later than given ending year, '
- + str(eyear))
+ P.error('Given starting year {} is later than given ending year,\ {}'.format(syear, eyear))
  else:
  pass

src/python/monsoon_sperber/lib/calc_metrics.py

@@ -7,25 +7,27 @@
 
 Jiwoo Lee, 2018-07
 
-Note: Code for picking onset/decay index inspired by 
+Note: Code for picking onset/decay index inspired by
 https://stackoverflow.com/questions/2236906/first-python-list-index-greater-than-x
 """
 
 import MV2
 
+
 def sperber_metrics(d, region, debug=False):
  """ d: input, 1d array of cumulative pentad time series """
  # Convert accumulation to fractional accumulation; normalize by sum
  d_sum = d[-1]
  frac_accum = MV2.divide(d, d_sum)
- onset_index = next(i for i,v in enumerate(frac_accum) if v >= 0.2)
+ onset_index = next(i for i, v in enumerate(frac_accum) if v >= 0.2)
  if region == 'GoG':
  decay_threshold = 0.6
  else:
  decay_threshold = 0.8
- decay_index = next(i for i,v in enumerate(frac_accum) if v >= decay_threshold)
- slope = ( frac_accum[decay_index] - frac_accum[onset_index] ) \
- / float(decay_index - onset_index)
+ decay_index = next(i for i, v in enumerate(
+ frac_accum) if v >= decay_threshold)
+ slope = (frac_accum[decay_index] - frac_accum[onset_index]) \
+ / float(decay_index - onset_index)
  return {'frac_accum': frac_accum,
  'onset_index': onset_index,
  'decay_index': decay_index,

src/python/monsoon_sperber/lib/divide_chunks.py

@@ -7,80 +7,82 @@
 Code taken from https://www.geeksforgeeks.org/break-list-chunks-size-n-python/
 """
 # Yield successive n-sized
-# chunks from l.
-def divide_chunks(l, n):
- # looping till length l
- for i in range(0, len(l), n):
- yield l[i:i+n]
+# chunks from data.
+
+
+def divide_chunks(data, n):
+ # looping till length data
+ for i in range(0, len(data), n):
+ yield data[i:i+n]
+
 
 """ Above code advanced considering leap year
 """
-def divide_chunks_advanced(l, n, debug=False):
+
+
+def divide_chunks_advanced(data, n, debug=False):
  # Double check first date should be Jan 1 (except for SH monsoon)
- tim = l.getTime()
+ tim = data.getTime()
  calendar = tim.calendar
  month = tim.asComponentTime()[0].month
  day = tim.asComponentTime()[0].day
- if debug: print('debug: first day of year is '+str(month)+'/'+str(day))
- if month not in [1,7] or day != 1:
- sys.exit('error: first day of year time series is '+str(month)+'/'+str(day))
+ if debug:
+ print('debug: first day of year is '+str(month)+'/'+str(day))
+ if month not in [1, 7] or day != 1:
+ sys.exit('error: first day of year time series is ' +
+ str(month)+'/'+str(day))
 
  # Check number of days in given year
- nday = len(l)
+ nday = len(data)
 
  if nday in [365, 360]:
- # looping till length l
- for i in range(0, len(l), n):
- yield l[i:i+n]
+ # looping till length data
+ for i in range(0, nday, n):
+ yield data[i:i+n]
  elif nday == 366:
  # until leap year day detected
- for i in range(0, len(l), n):
+ for i in range(0, nday, n):
  # Check if leap year date included
  leap_detect = False
  for ii in range(i, i+n):
- date = l.getTime().asComponentTime()[ii]
+ date = data.getTime().asComponentTime()[ii]
  month = date.month
  day = date.day
  if month == 2 and day > 28:
- if debug: 
+ if debug:
  print('debug: leap year detected:', month, '/', day)
  leap_detect = True
  if leap_detect:
- yield l[i:i+n+1]
+ yield data[i:i+n+1]
  tmp = i+n+1
  break
  else:
- yield l[i:i+n]
+ yield data[i:i+n]
  # after leap year day passed
  if leap_detect:
- for i in range(tmp, len(l), n):
- yield l[i:i+n]
+ for i in range(tmp, nday, n):
+ yield data[i:i+n]
  elif nday == 361 and calendar == '360_day':
- # Speacial case handling for HadGEM2 family where time bounds was not 
+ # Speacial case handling for HadGEM2 family where time bounds was not
  # properly saved, so include next year's first day in time series
  if debug:
- print('debug: 361 to 360 revise: l.shape, before:', l.shape)
- l = l[0:360]
+ print('debug: 361 to 360 revise: l.shape, before:', data.shape)
+ data = data[0:360]
  if debug:
- print('debug: 361 to 360 revise: l.shape, after:', l.shape)
- # looping till length l
- for i in range(0, len(l), n):
- yield l[i:i+n]
- #elif nday == 360:
- # l2 = interp1d(l, 365, debug=debug)
- # # looping till length l
- # for i in range(0, len(l2), n):
- # yield l2[i:i+n]
+ print('debug: 361 to 360 revise: l.shape, after:', data.shape)
+ # looping till length data
+ for i in range(0, nday, n):
+ yield data[i:i+n]
  else:
  sys.exit('error: number of days in year is '+str(nday))
 
 
-def interp1d(l, ref_length, debug=False):
- l = np.array(l)
- l_interp = interp.interp1d(np.arange(l.size), l)
- l2 = l_interp(np.linspace(0, l.size-1, ref_length))
+def interp1d(data, ref_length, debug=False):
+ data = np.array(data)
+ data_interp = interp.interp1d(np.arange(data.size), data)
+ data2 = data_interp(np.linspace(0, data.size-1, ref_length))
  if debug:
  print('debug: 1d interpolation')
- print('debug: length before interp: ', len(l))
- print('debug: length after interp: ', len(l2))
- return l2
+ print('debug: length before interp: ', len(data))
+ print('debug: length after interp: ', len(data2))
+ return data2