Merge pull request #510 from PCMDI/diurnal

Green lights!

run_tests.py

@@ -186,7 +186,9 @@ def run_nose(test_name):
  sys.exit(0)
 
 # Make sure we have sample data
-#cdat_info.download_sample_data_files(os.path.join(sys.prefix,"share","vcs","test_data_files.txt"),cdat_info.get_sampledata_path())
+if not os.path.exists("test_data"):
+ os.makedirs("test_data")
+cdat_info.download_sample_data_files(os.path.join(sys.prefix,"share","pmp","test_data_files.txt"),"test_data")
 if args.update:
  os.environ["UPDATE_TESTS"]="True"
 if args.traceback:

setup.py

@@ -50,6 +50,7 @@
 packages = {'pcmdi_metrics': 'src/python',
  'pcmdi_metrics.io': 'src/python/io',
  'pcmdi_metrics.pcmdi': 'src/python/pcmdi',
+ 'pcmdi_metrics.diurnal': 'src/python/diurnal',
  'pcmdi_metrics.graphics': 'src/python/graphics',
  'pcmdi_metrics.driver': 'src/python/pcmdi/scripts/driver',
  }
@@ -60,6 +61,7 @@
  'demo/pmp_demo_1.py',
  'demo/pmp_demo.py',
  ]
+scripts += glob.glob("src/python/diurnal/scripts/*.py")
 
 demo_files = glob.glob("demo/*/*")
 print "demo files"
@@ -87,6 +89,7 @@
  ('share/pmp', ('doc/obs_info_dictionary.json',
  'share/pcmdi_metrics_table',
  'share/disclaimer.txt',
+ 'share/test_data_files.txt',
  'share/default_regions.py'
  )),
  ('share/pmp/demo', demo_files),
@@ -139,4 +142,4 @@
  # extra_compile_args = [],
  # extra_link_args = [],
  # ]
- )
+ )

share/default_regions.py

@@ -1,44 +1,44 @@
 import cdutil
 
 regions_specs = {
- 'NHEX': {'domain': cdutil.region.domain(latitude=(30., 90, 'ccb'))},
- 'SHEX': {'domain': cdutil.region.domain(latitude=(-90., -30, 'ccb'))},
- 'TROPICS': {'domain': cdutil.region.domain(latitude=(-30., 30, 'ccb'))},
+ 'NHEX': {'domain': cdutil.region.domain(latitude=(30., 90))},
+ 'SHEX': {'domain': cdutil.region.domain(latitude=(-90., -30))},
+ 'TROPICS': {'domain': cdutil.region.domain(latitude=(-30., 30))},
  "global": {},
- '90S50S': {'domain': cdutil.region.domain(latitude=(-90., -50, 'ccb'))},
- '50S20S': {'domain': cdutil.region.domain(latitude=(-50., -20, 'ccb'))},
- '20S20N': {'domain': cdutil.region.domain(latitude=(-20., 20, 'ccb'))},
- '20N50N': {'domain': cdutil.region.domain(latitude=(20., 50, 'ccb'))},
- '50N90N': {'domain': cdutil.region.domain(latitude=(50., 90, 'ccb'))},
- 'land_NHEX': {'value': 100, 'domain': cdutil.region.domain(latitude=(30., 90, 'ccb'))},
- 'land_SHEX': {'value': 100, 'domain': cdutil.region.domain(latitude=(-90., -30, 'ccb'))},
- 'land_TROPICS': {'value': 100, 'domain': cdutil.region.domain(latitude=(-30., 30, 'ccb'))},
+ '90S50S': {'domain': cdutil.region.domain(latitude=(-90., -50))},
+ '50S20S': {'domain': cdutil.region.domain(latitude=(-50., -20))},
+ '20S20N': {'domain': cdutil.region.domain(latitude=(-20., 20))},
+ '20N50N': {'domain': cdutil.region.domain(latitude=(20., 50))},
+ '50N90N': {'domain': cdutil.region.domain(latitude=(50., 90))},
+ 'land_NHEX': {'value': 100, 'domain': cdutil.region.domain(latitude=(30., 90))},
+ 'land_SHEX': {'value': 100, 'domain': cdutil.region.domain(latitude=(-90., -30))},
+ 'land_TROPICS': {'value': 100, 'domain': cdutil.region.domain(latitude=(-30., 30))},
  "land": {'value': 100, },
- 'ocean_NHEX': {'value': 0, 'domain': cdutil.region.domain(latitude=(30., 90, 'ccb'))},
- 'ocean_SHEX': {'value': 0, 'domain': cdutil.region.domain(latitude=(-90., -30, 'ccb'))},
- 'ocean_TROPICS': {'value': 0, 'domain': cdutil.region.domain(latitude=(30., 30, 'ccb'))},
+ 'ocean_NHEX': {'value': 0, 'domain': cdutil.region.domain(latitude=(30., 90))},
+ '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
- "NAM": {'domain': cdutil.region.domain(latitude=(20., 90, 'ccb'), longitude=(-180, 180, 'ccb'))},
- "NAO": {'domain': cdutil.region.domain(latitude=(20., 80, 'ccb'), longitude=(-90, 40, 'ccb'))},
- "SAM": {'domain': cdutil.region.domain(latitude=(-20., -90, 'ccb'), longitude=(0, 360, 'ccb'))},
- "PNA": {'domain': cdutil.region.domain(latitude=(20., 85, 'ccb'), longitude=(120, 240, 'ccb'))},
- "PDO": {'domain': cdutil.region.domain(latitude=(20., 70, 'ccb'), longitude=(110, 260, 'ccb'))},
+ "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
  # All monsoon domains
- 'AllM': {'domain': cdutil.region.domain(latitude=(-45., 45., 'ccb'), longitude=(0., 360., 'ccb'))},
+ 'AllM': {'domain': cdutil.region.domain(latitude=(-45., 45.), longitude=(0., 360.))},
  # North American Monsoon
- 'NAMM': {'domain': cdutil.region.domain(latitude=(0., 45., 'ccb'), longitude=(210., 310., 'ccb'))},
+ 'NAMM': {'domain': cdutil.region.domain(latitude=(0., 45.), longitude=(210., 310.))},
  # South American Monsoon
- 'SAMM': {'domain': cdutil.region.domain(latitude=(-45., 0., 'ccb'), longitude=(240., 330., 'ccb'))},
+ 'SAMM': {'domain': cdutil.region.domain(latitude=(-45., 0.), longitude=(240., 330.))},
  # North African Monsoon
- 'NAFM': {'domain': cdutil.region.domain(latitude=(0., 45., 'ccb'), longitude=(310., 60., 'ccb'))},
+ 'NAFM': {'domain': cdutil.region.domain(latitude=(0., 45.), longitude=(310., 60.))},
  # South African Monsoon
- 'SAFM': {'domain': cdutil.region.domain(latitude=(-45., 0., 'ccb'), longitude=(0., 90., 'ccb'))},
+ 'SAFM': {'domain': cdutil.region.domain(latitude=(-45., 0.), longitude=(0., 90.))},
  # Asian Summer Monsoon
- 'ASM': {'domain': cdutil.region.domain(latitude=(0., 45., 'ccb'), longitude=(60., 180., 'ccb'))},
+ 'ASM': {'domain': cdutil.region.domain(latitude=(0., 45.), longitude=(60., 180.))},
  # Australian Monsoon
- 'AUSM': {'domain': cdutil.region.domain(latitude=(-45., 0., 'ccb'), longitude=(90., 160., 'ccb'))},
+ 'AUSM': {'domain': cdutil.region.domain(latitude=(-45., 0.), longitude=(90., 160.))},
 }
 
 default_regions = ['global', 'NHEX', 'SHEX', 'TROPICS']

share/test_data_files.txt

@@ -0,0 +1,2 @@
+http:https://uvcdat.llnl.gov/cdat/pmp
+60582c6986451ffd4925a303c4f0040a sample_data_pr_CMCC.nc

src/python/diurnal/__init__.py

@@ -0,0 +1,2 @@
+from . import common # noqa
+from . import fourierFFT # noqa

src/python/diurnal/common.py

@@ -0,0 +1,64 @@
+import genutil
+import argparse
+from pcmdi_metrics.driver.pmp_parser import PMPParser
+
+monthname_d = {1: 'Jan', 2: 'Feb', 3: 'Mar', 4: 'Apr', 5: 'May', 6: 'Jun',
+ 7: 'Jul', 8: 'Aug', 9: 'Sep', 10: 'Oct', 11: 'Nov', 12: 'Dec'}
+
+
+class INPUT(object):
+ def __init__(self, args, filename, filename_template, varname="pr"):
+ self.fileName = filename
+ self.args = args
+ self.monthname = monthname_d[args.month]
+ self.varname = varname
+ self.template = filename_template
+
+
+def populateStringConstructor(template, args):
+ template = genutil.StringConstructor(template)
+ for k in template.keys():
+ if hasattr(args, k):
+ setattr(template, k, str(getattr(args, k)))
+ return template
+
+
+P = PMPParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+P.add_argument("-i", "--modroot",
+ default='data',
+ help="Root directory for model (or observed) 3-hourly data")
+P.add_argument("-m", "--month",
+ type=int,
+ default=7,
+ help="Month to be processed, given as 2-char month number")
+P.add_argument("-f", "--firstyear",
+ type=int,
+ default=1999,
+ help="First year of data processing")
+P.add_argument("-l", "--lastyear",
+ type=int,
+ default=2005,
+ help="Last year of data processing")
+P.add_argument("-o", "--output_directory",
+ default="out",
+ help="output directory")
+P.add_argument("-r", "--realization",
+ default="r1i1p1",
+ help="Realization used")
+P.add_argument("-w", "--num-workers", default=None, type=int,
+ help="number of workers to use in multiprocessing, 0 means auto")
+P.add_argument("--version", default="*")
+P.add_argument("--frequency", default="3hr")
+P.add_argument("--realm", default="atm")
+P.add_argument("--model", default="*")
+P.add_argument("--experiment", default="historical")
+P.add_argument("-t", "--filename_template",
+ default="cmip5.%(model).%(experiment).%(realization).%(frequency).%(realm).%(frequency).%(variable).%(version).latestX.xml") # noqa
+P.add_argument("--skip", default=[],
+ help="models to skip", nargs="*")
+P.add_argument(
+ "-a",
+ "--append",
+ default=False,
+ action="store_true",
+ help="append in json file if json exist (e.g. adding a model to file)")

src/python/diurnal/fourierFFT.py

@@ -0,0 +1,88 @@
+def fastFT(x, t):
+ '''
+Use a Numerical Python function to compute a FAST Fourier transform -- which should give the same result as a simple
+SLOW Fourier integration via the trapezoidal rule.
+
+Return mean + amplitudes and times-of-maximum of the first three Fourier harmonic components of a time series x(t).
+Do NOT detrend the time series first, in order to retain the "sawtooth" frequency implied by the input length of the
+time series (e.g. the 24-hour period from a composite-diurnal cycle).
+
+On input: x[i, j] = values at each gridpoint (i) for N times (j), e.g. N = 8 for a 3-hr composite-diurnal cycle
+ t[i, j] = timepoints at each gridpoint (i) for N times (j), e.g. Local Standard Times
+
+On output: c[i] = mean value at each gridpoint (i) in the time series (= constant term in Fourier series)
+ maxvalue[i, k] = amplitude at each gridpoint (i) for each Fourier harmonic (k)
+ tmax [i, k] = time of maximum at each gridpoint (i) for each Fourier harmonic (k)
+
+ Curt Covey, PCMDI/LLNL November 2016
+ (from ./fourier.py and ./fourierTestFFT.py)
+ '''
+ import numpy
+ nGridPoints = len(x)
+ X = numpy.fft.ifft(x)
+ a = X.real
+ b = X.imag
+ S = numpy.sqrt(a**2 + b**2)
+ c = S[:, 0]
+ # time of maximum for nth component (n=0 => diurnal, n=1 => semi...)
+ tmax = numpy.zeros((nGridPoints, 3))
+ # value of maximum for nth component (= 1/2 peak-to-peak amplitude)
+ maxvalue = numpy.zeros((nGridPoints, 3))
+ for n in range(3):
+ # Adding first + last terms, second + second-to-last, ...
+ maxvalue[:, n] = S[:, n + 1] + S[:, -n - 1]
+ tmax[:, n] = numpy.arctan2(b[:, n + 1], a[:, n + 1])
+ tmax[:, n] = tmax[:, n] * 12.0 / \
+ (numpy.pi * (n + 1)) # Radians to hours
+ tmax[:, n] = tmax[:, n] + t[:, 0] # GMT to LST
+ tmax[:, n] = tmax[:, n] % (24 / (n + 1))
+ return c, maxvalue, tmax
+
+
+def fastAllGridFT(x, t):
+ '''
+This version of fastFT (see above) does all gridpoints at once.
+
+Use a Numerical Python function to compute a FAST Fourier transform -- which should give the same result as a simple
+SLOW Fourier integration via the trapezoidal rule.
+
+Return mean + amplitudes and times-of-maximum of the first three Fourier harmonic components of a time series x(t).
+Do NOT detrend the time series first, in order to retain the "sawtooth" frequency implied by the input length of the
+time series (e.g. the 24-hour period from a composite-diurnal cycle).
+
+On input: x[k,i,j] = values at each gridpoint (i,j) for N times (k), e.g. N = 8 for a 3-hr composite-diurnal cycle
+ t[k,i,j] = timepoints at each gridpoint (i,j) for N times (k), e.g. Local Standard Times
+
+On output: c[i,j] = mean value at each gridpoint (i,j) in the time series ("zeroth" term in Fourier series)
+ maxvalue[n,i,j] = amplitude at each gridpoint (i,j) for each Fourier harmonic (n)
+ tmax [n,i,j] = time of maximum at each gridpoint (i,j) for each Fourier harmonic (n)
+
+ Curt Covey, PCMDI/LLNL December 2016
+ '''
+ import numpy
+
+ print 'Creating output arrays ...'
+ nx = x.shape[1]
+ ny = x.shape[2]
+ # time of maximum for nth component (n=0 => diurnal, n=1 => semi...)
+ tmax = numpy.zeros((3, nx, ny))
+ # value of maximum for nth component (= 1/2 peak-to-peak amplitude)
+ maxvalue = numpy.zeros((3, nx, ny))
+
+ print 'Calling numpy FFT function ...'
+ X = numpy.fft.ifft(x, axis=0)
+ print X.shape
+
+ print 'Converting from complex-valued FFT to real-valued amplitude and phase ...'
+ a = X.real
+ b = X.imag
+ S = numpy.sqrt(a**2 + b**2)
+ c = S[0] # Zeroth harmonic = mean-value "constant term" in Fourier series.
+ for n in range(3):
+ # Adding first + last terms, second + second-to-last, ...
+ maxvalue[n] = S[n + 1] + S[-n - 1]
+ tmax[n] = numpy.arctan2(b[n + 1], a[n + 1])
+ tmax[n] = tmax[n] * 12.0 / (numpy.pi * (n + 1)) # Radians to hours
+ tmax[n] = tmax[n] + t[0] # GMT to LST
+ tmax[n] = tmax[n] % (24 / (n + 1))
+ return c, maxvalue, tmax