404 jwl pentmonsoon (#563)

* place holder

* renamed

* place holder

* renamed

* extracting time series for each year

* add getting precip ave from pentad chunk

* clean up conflicts

* add testing plot

* minor bug fix for plot file name

* add for testing purpose archive

* rename plot file

* rename plot file saving

* minor change for test plot name

* add composite of year by year pentad time sereies

* minor clean up

* debug test plot for multiple lines

* Delete cmip5_CNRM-CM5_historical_r9i1p1_ASM_1850.png

* add sample nc output for develop process

* clarify variable name for monsoon regions

* minor fix

* use one testing plot for multiple monsoon regions

* unify coordinate axis between individual year and composite time series

* minor fix

* test with multiple models

* fix for testing plot when analyzing multiple model

* minor fix for plot legend

* minor: fix legend box location to lower right

* add result example for develop process

* start merging pmpparser capability

* download complete for CPC daily ref dataset

* issue list added in classical way

* add comment for issue list note

* add comment on issue note

* minor revise for comment

* comment adding for issue list note

* land only (mask out ocean area)

* d.unit to d.units

* reduce redundancy of type conversion to MV2.array for pentad_time_series

* reduce hardcorded part

* merge debugging part to shorten code

* add issue list

* rename testing output dir

* working on leap year issue

* add Sperber monsoon domains

* move divide chunk to lib and draft leap year handling

* advance chunk dividing, but still need to correct

* bug fix for leap year chunk getting

* bug fix for leap year chunk getting

* clean up

* remove redundant debug print

* mark completed task in issue list

* minor: plot legend location change

* test modpath, update issue list

* reduce debug print

* update model list as available

* Get modnames from myParam.py

* test with multi runs

* add numpy cumsum for cumulative pentad time series

* clean up

* check with r1i1p1 when testing

* add JSON and fractional accumulation for metrics calculation

* add metrics calc

* minor grammer fix

* add copyfile import

* add json import

* add description at top

* advance json structure

* disable try statement when debugging

* minor fix for chunk average

* test run

* test run

* plot option from myParam.py applied for plotting part

* test with multiple regions

* advance print statements

* clean up

* clean up

* code clean up

* add map to visualize geographical location of monsoon domains

* show plot legend for only one panel

* GoG domain monsoon decay at 60%

* start pendtad ts from 7/1 for SH monsoon domains

* shorten code

* advance

* add check time spend for each run

* add check time spend for each run (typo fix)

* advance code efficiency (load lf once per model, close xml files after used)

* minor fix

* bug fix for plotting

* simplify plot and code for it

* BNU-ESM returned axis no match error: make sure array to have consistent axes after get through the function

* plot label fix

* code clean up

* put place holder for 72 to 73 interpolation for 360 day/year model (e.g., HadGEM2 family)

* add 360 day calendar handling: 1d interpolation

* minor bug fix for string joining

* adjust 1d interp process

* minor fix

* minor fix 2

* update issue list

* tweaks to cleanup PR before going to master

* partial fixes to @lee1043 PR

* example of string contructor template as argument

* move domain check plotting files to doc directory

* Delete monsoon_domain_map.png

moved to doc

* Delete plot_monsoon_domain.py

moved to doc

* remove test output files

* few changes: add sperber_metrics for importing, replace StringConstructor to param.process_templated_argument

* some clean up and add inline document

* clean up and typo fix

* Include OBS dataset for metric calculation

* bug fix for when daily file saved as e.g., 1997-1-1 0:0:0.0

* bug fix for when daily file saved as e.g., 1997-1-1 0:0:0.0 (cont)

* obs capability add bug fix

* to work with obs

* turn off debugging mode

* enable option for run across cmip5

* add switch to decide include obs or not

* clean up

* add import sys for case being called

* HadGEM2 model family 360 calendar handling case when bound not properly saved

* clean up

* define modpath_lf for parser

* add one more try level for model, not only for realizations, because there was a case that failure of loading modpath_lf caused entire loop stop

* typo fix

* use relative path for default_regions.py when testing purpose

* add import numpy as np

* add debugging message

* json structure revised

* Improve start and end year getting part

* minor: add inline description

* improve inline description

* minor: remove repeated blank line

* remove importing divide_chunks which no more used

* move f_lf close upward

* simplify if statement for plot and nc_out

* clean up inline description

* jiwoo branch ready?

* mac ffmpeg issue fix

* travis mac ffmpeg fix

* add debug print

* fix little flake8 tweak

* I had FFMPEG env defined at the wrong place for circleci and not used for travis...

* magic conda pkg combination?

* travis needs scipy

* correct bash...

.circleci/config.yml

@@ -22,8 +22,8 @@ aliases:
  conda config --set always_yes yes --set changeps1 no
  conda update -y -q conda
  conda config --set anaconda_upload no
- conda create -q -n py2 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib cia testsrunner "proj4<5" "python<3"
- conda create -q -n py3 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib cia testsrunner "proj4<5" "python>3"
+ conda create -q -n py2 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib matplotlib scipy cia testsrunner "proj4<5" "python<3" $FFMPEG
+ conda create -q -n py3 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib matplotlib scipy cia testsrunner "proj4<5" "python>3" $FFMPEG
 
  - &setup_pmp
  name: setup_pmp
@@ -101,6 +101,7 @@ jobs:
  environment:
  WORKDIR: "workspace/test_macos_pmp"
  OS: "osx-64"
+ FFMPEG: "'ffmpeg>4' 'libpng>1.6.34'"
  steps:
  - checkout
  - run: *setup_miniconda
@@ -120,6 +121,7 @@ jobs:
  environment:
  WORKDIR: "workspace/test_linux_pmp"
  OS: "linux-64"
+ FFMPEG: "ffmpeg"
  steps:
  - checkout
  - run: *setup_miniconda

.travis.yml

@@ -16,9 +16,11 @@ before_install:
  - conda update -y -q conda
 
 script:
- - export UVCDAT_ANONYMOUS_LOG="False"
- - conda create -q -n py2 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib testsrunner cia "proj4<5" "python<3"
- - conda create -q -n py3 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib testsrunner cia "proj4<5" "python>3"
+ - export CDAT_ANONYMOUS_LOG="False"
+ - export FFMPEG="ffmpeg"
+ - if [ "$TRAVIS_OS_NAME" = "osx" ]; then export FFMPEG="'ffmpeg>4' 'libpng>1.6.34'"; fi
+ - conda create -q -n py2 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib testsrunner cia "proj4<5" "python<3" scipy $FFMPEG
+ - conda create -q -n py3 -c cdat/label/nightly -c conda-forge -c cdat -c pcmdi vcs vcsaddons mesalib testsrunner cia "proj4<5" "python>3" scipy $FFMPEG
  # Useful for debugging any issues with conda
  - conda info -a
  - source activate py2

setup.py

@@ -55,6 +55,7 @@
  'pcmdi_metrics.graphics': 'src/python/graphics',
  'pcmdi_metrics.driver': 'src/python/pcmdi/scripts/driver',
  'pcmdi_metrics.monsoon_wang': 'src/python/monsoon_wang/lib',
+ 'pcmdi_metrics.monsoon_sperber': 'src/python/monsoon_sperber/lib',
  }
 scripts = ['src/python/pcmdi/scripts/pcmdi_metrics_driver.py',
  'src/python/pcmdi/scripts/pcmdi_metrics_driver_legacy.py',
@@ -64,6 +65,7 @@
  'demo/pmp_demo_1.py',
  'demo/pmp_demo.py',
  'src/python/monsoon_wang/scripts/mpindex_compute.py',
+ 'src/python/monsoon_sperber/scripts/driver_monsoon_sperber.py',
  ]
 scripts += glob.glob("src/python/diurnal/scripts/*.py")
 

share/default_regions.py

@@ -18,13 +18,15 @@
  'ocean_SHEX': {'value': 0, 'domain': cdutil.region.domain(latitude=(-90., -30))},
  'ocean_TROPICS': {'value': 0, 'domain': cdutil.region.domain(latitude=(30., 30))},
  "ocean": {'value': 0, },
- # Below is for modes of variability
+
+ # Modes of variability
  "NAM": {'domain': cdutil.region.domain(latitude=(20., 90), longitude=(-180, 180))},
  "NAO": {'domain': cdutil.region.domain(latitude=(20., 80), longitude=(-90, 40))},
  "SAM": {'domain': cdutil.region.domain(latitude=(-20., -90), longitude=(0, 360))},
  "PNA": {'domain': cdutil.region.domain(latitude=(20., 85), longitude=(120, 240))},
  "PDO": {'domain': cdutil.region.domain(latitude=(20., 70), longitude=(110, 260))},
- # Below is for monsoon domains
+
+ # Monsoon domains for Wang metrics
  # All monsoon domains
  'AllMW': {'domain': cdutil.region.domain(latitude=(-40., 45.), longitude=(0., 360.))},
  'AllM': {'domain': cdutil.region.domain(latitude=(-45., 45.), longitude=(0., 360.))},
@@ -40,6 +42,20 @@
  'ASM': {'domain': cdutil.region.domain(latitude=(0., 45.), longitude=(60., 180.))},
  # Australian Monsoon
  'AUSM': {'domain': cdutil.region.domain(latitude=(-45., 0.), longitude=(90., 160.))},
+
+ # Monsoon domains for Sperber metrics
+ # All India rainfall
+ 'AIR': {'domain': cdutil.region.domain(latitude=(7., 25.), longitude=(65., 85.))},
+ # North Australian
+ 'AUS': {'domain': cdutil.region.domain(latitude=(-20., -10.), longitude=(120., 150.))},
+ # Sahel
+ 'Sahel': {'domain': cdutil.region.domain(latitude=(13., 18.), longitude=(-10., 10.))},
+ # Gulf of Guinea 
+ 'GoG': {'domain': cdutil.region.domain(latitude=(0., 5.), longitude=(-10., 10.))},
+ # North American monsoon
+ 'NAmo': {'domain': cdutil.region.domain(latitude=(20., 37.), longitude=(-112., -103.))},
+ # South American monsoon
+ 'SAmo': {'domain': cdutil.region.domain(latitude=(-20., 2.5), longitude=(-65., -40.))},
 }
 
 default_regions = ['global', 'NHEX', 'SHEX', 'TROPICS']

share/test_data_files.txt

@@ -1,2 +1,4 @@
 https://cdat.llnl.gov/cdat/pmp
-60582c6986451ffd4925a303c4f0040a sample_data_pr_CMCC.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

@@ -0,0 +1,53 @@
+
+# =================================================
+# Background Information
+# -------------------------------------------------
+mip = 'cmip5'
+exp = 'historical'
+frequency = 'da'
+realm = 'atm'
+
+# =================================================
+# 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' # 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
+
+osyear = 1996
+oeyear = 2016
+
+includeOBS = True
+
+# =================================================
+# Models
+# -------------------------------------------------
+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'] # noqa
+
+realization = 'r1i1p1'
+
+varModel = 'pr'
+ModUnitsAdjust = (True, 'multiply', 86400.0) # kg m-2 s-1 to mm day-1
+units = 'mm/d'
+
+msyear = 1961
+meyear = 1999
+
+# =================================================
+# Output
+# -------------------------------------------------
+results_dir = '/work/lee1043/imsi/result_test/monsoon_sperber'
+nc_out = True # Write output in NetCDF
+plot = True # Create map graphics
+
+# =================================================
+# Miscellaneous
+# -------------------------------------------------
+update_json = True
+debug = False

src/python/monsoon_sperber/doc/plot_monsoon_domain.py

@@ -0,0 +1,47 @@
+from mpl_toolkits.basemap import Basemap
+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)
+ plt.gca().add_patch(poly)
+
+
+lats = {}
+lons = {}
+
+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['Sahel'] = [13, 18, 18, 13]
+lons['Sahel'] = [-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['SAM'] = [-20, -2.5, -2.5, -20]
+lons['SAM'] = [-65, -65, -40, -40]
+
+regions = lats.keys()
+
+plt.figure(figsize=(12, 6))
+
+m = Basemap(lon_0=0)
+m.drawmapboundary(fill_color='aqua')
+m.fillcontinents(color='grey', lake_color='aqua')
+m.drawcoastlines()
+
+for region in regions:
+ draw_screen_poly(lats[region], lons[region], m)
+
+plt.savefig('monsoon_domain_map.png')

src/python/monsoon_sperber/lib/__init__.py

@@ -0,0 +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_advanced, interp1d # noqa

src/python/monsoon_sperber/lib/argparse_functions.py

@@ -0,0 +1,71 @@
+def AddParserArgument(P):
+ # Load pre-defined parsers
+ P.use("--mip")
+ P.use("--exp")
+ P.use("--results_dir")
+ P.use("--reference_data_path")
+ P.use("--modpath")
+ # Add parsers for options
+ P.add_argument("--frequency", default="da")
+ P.add_argument("--realm", default="atm")
+ P.add_argument("--reference_data_name",
+ type=str,
+ help="Name of reference data set")
+ P.add_argument("--reference_data_lf_path",
+ type=str,
+ help="Path of landsea mask for reference data set")
+ 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,
+ 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,
+ 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("--modnames",
+ type=list,
+ default=None,
+ help="List of models")
+ P.add_argument("-r", "--realization",
+ type=str,
+ default="r1i1p1",
+ help="Consider all accessible realizations as idividual\n"
+ "- r1i1p1: default, consider only 'r1i1p1' member\n"
+ " Or, specify realization, e.g, r3i1p1'\n"
+ "- *: consider all available realizations")
+ # Add parsers as switches
+ P.add_argument("-d", "--debug",
+ 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)
+ return P
+
+
+def YearCheck(syear, eyear, P):
+ if syear >= 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

@@ -0,0 +1,34 @@
+""" Calculate metrics based on Sperber and Annamalai 2014 Clim. Dyn.,
+which are:
+- onset pentad index: when fractional accumulation hit 20%
+- decay pentad index: when fractional accumulation hit 80%
+- slope: slope between onset and decay pentad time step indices
+calculated from cumulative pentad time series
+
+Jiwoo Lee, 2018-07
+
+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)
+ 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)
+ return {'frac_accum': frac_accum,
+ 'onset_index': onset_index,
+ 'decay_index': decay_index,
+ 'slope': slope}