diff --git a/doc/Diurnal Cycle Diagram.pdf b/doc/Diurnal Cycle Diagram.pdf new file mode 100644 index 000000000..cb6fa79b1 Binary files /dev/null and b/doc/Diurnal Cycle Diagram.pdf differ diff --git a/run_tests.py b/run_tests.py index ba19c4453..b5beb5212 100755 --- a/run_tests.py +++ b/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: diff --git a/setup.py b/setup.py index 9ea1696df..c0c0ffc3a 100755 --- a/setup.py +++ b/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 = [], # ] - ) \ No newline at end of file + ) diff --git a/share/default_regions.py b/share/default_regions.py index f681cb8c7..12f2b22b7 100755 --- a/share/default_regions.py +++ b/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'] diff --git a/share/test_data_files.txt b/share/test_data_files.txt new file mode 100644 index 000000000..02c0c4873 --- /dev/null +++ b/share/test_data_files.txt @@ -0,0 +1,2 @@ +http://uvcdat.llnl.gov/cdat/pmp +60582c6986451ffd4925a303c4f0040a sample_data_pr_CMCC.nc diff --git a/src/python/diurnal/__init__.py b/src/python/diurnal/__init__.py new file mode 100644 index 000000000..04d4c2971 --- /dev/null +++ b/src/python/diurnal/__init__.py @@ -0,0 +1,2 @@ +from . import common # noqa +from . import fourierFFT # noqa diff --git a/src/python/diurnal/common.py b/src/python/diurnal/common.py new file mode 100644 index 000000000..f24aa3190 --- /dev/null +++ b/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)") diff --git a/src/python/diurnal/fourierFFT.py b/src/python/diurnal/fourierFFT.py new file mode 100755 index 000000000..9549d1e69 --- /dev/null +++ b/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 diff --git a/src/python/diurnal/scripts/compositeDiurnalStatisticsWrapped.py b/src/python/diurnal/scripts/compositeDiurnalStatisticsWrapped.py new file mode 100755 index 000000000..b2b275974 --- /dev/null +++ b/src/python/diurnal/scripts/compositeDiurnalStatisticsWrapped.py @@ -0,0 +1,225 @@ +#!/usr/bin/env python + +# /export/covey1/CMIP5/Precipitation/DiurnalCycle/HistoricalRuns/compositeDiurnalStatisticsWrapped.py + +# This modifiction of ./compositeDiurnalStatistics.py will have the PMP Parser "wrapped" around it, +# so that it can be executed with input parameters in the Unix command line, for example: +# ---> python compositeDiurnalStatisticsWrapped.py -t "sample_data_%(variable)_%(model).nc" -m 7 + +# These are the models with CMIP5 historical run output at 3h frequency, which this script is designed to process: +# 'ACCESS1-0', 'ACCESS1-3', 'bcc-csm1-1', 'bcc-csm1-1-m', 'BNU-ESM', +# 'CCSM4', 'CMCC-CM', 'CNRM-CM5', 'EC-EARTH', +# 'FGOALS-g2', 'GFDL-CM3', 'GFDL-ESM2M', 'GISS-E2-R', +# 'GISS-E2-H', 'inmcm4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', +# 'MIROC4h', 'MIROC5', 'MIROC-ESM', 'MIROC-ESM-CHEM' + +import cdms2 +import genutil +import MV2 +import os +import sys +import glob +import cdtime +import cdp + +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor, INPUT + + +def compute(params): + fileName = params.fileName + month = params.args.month + monthname = params.monthname + varbname = params.varname + template = populateStringConstructor(args.filename_template, args) + template.variable = varbname + # Units on output (*may be converted below from the units of input*) + outunits = 'mm/d' + startime = 1.5 # GMT value for starting time-of-day + + reverted = template.reverse(os.path.basename(fileName)) + dataname = reverted["model"] + if dataname not in args.skip: + try: + print 'Data source:', dataname + print 'Opening %s ...' % fileName + f = cdms2.open(fileName) + +# Composite-mean and composite-s.d diurnal cycle for month and year(s): + iYear = 0 + for year in range(args.firstyear, args.lastyear + 1): + print 'Year %s:' % year + startTime = cdtime.comptime(year, month, 1, 1, 30) + # Last possible second to get all tpoints + finishtime = startTime.add( + 1, cdtime.Month).add(-1.5, cdtime.Hour).add(.1, cdtime.Second) + print 'Reading %s from %s for time interval %s to %s ...' % (varbname, fileName, startTime, finishtime) + # Transient variable stores data for current year's month. + tvarb = f(varbname, time=(startTime, finishtime)) + # *HARD-CODES conversion from kg/m2/sec to mm/day. + tvarb *= 86400 + print 'Shape:', tvarb.shape + # The following tasks need to be done only once, extracting + # metadata from first-year file: + if year == args.firstyear: + tc = tvarb.getTime().asComponentTime() + day1 = cdtime.comptime(tc[0].year, tc[0].month) + firstday = tvarb( + time=( + day1, + day1.add( + 1, + cdtime.Day), + "con")) + dimensions = firstday.shape + # print ' Shape = ', dimensions + # Number of time points in the selected month for one year + N = dimensions[0] + nlats = dimensions[1] + nlons = dimensions[2] + deltaH = 24. / N + dayspermo = tvarb.shape[0] / N + print ' %d timepoints per day, %d hr intervals between timepoints' % (N, deltaH) + comptime = firstday.getTime() + modellons = tvarb.getLongitude() + modellats = tvarb.getLatitude() + # Longitude values are needed later to compute Local Solar + # Times. + lons = modellons[:] + print ' Creating temporary storage and output fields ...' + # Sorts tvarb into separate GMTs for one year + tvslice = MV2.zeros((N, dayspermo, nlats, nlons)) + # Concatenates tvslice over all years + concatenation = MV2.zeros( + (N, dayspermo * nYears, nlats, nlons)) + LSTs = MV2.zeros((N, nlats, nlons)) + for iGMT in range(N): + hour = iGMT * deltaH + startime + print ' Computing Local Standard Times for GMT %5.2f ...' % hour + for j in range(nlats): + for k in range(nlons): + LSTs[iGMT, j, k] = (hour + lons[k] / 15) % 24 + for iGMT in range(N): + hour = iGMT * deltaH + startime + print ' Choosing timepoints with GMT %5.2f ...' % hour + # Transient-variable slice: every Nth tpoint gets all of + # the current GMT's tpoints for current year: + tvslice[iGMT] = tvarb[iGMT:tvarb.shape[0]:N] + concatenation[iGMT, iYear * + dayspermo: (iYear + + 1) * + dayspermo] = tvslice[iGMT] + iYear += 1 + f.close() + + # For each GMT, take mean and standard deviation over all years for + # the chosen month: + avgvalues = MV2.zeros((N, nlats, nlons)) + stdvalues = MV2.zeros((N, nlats, nlons)) + for iGMT in range(N): + hour = iGMT * deltaH + startime + print 'Computing mean and standard deviation over all GMT %5.2f timepoints ...' % hour + # Assumes first dimension of input ("axis#0") is time + avgvalues[iGMT] = MV2.average(concatenation[iGMT], axis=0) + stdvalues[iGMT] = genutil.statistics.std(concatenation[iGMT]) + avgvalues.id = 'diurnalmean' + stdvalues.id = 'diurnalstd' + LSTs.id = 'LST' + avgvalues.units = outunits + # Standard deviation has same units as mean (not so for + # higher-moment stats). + stdvalues.units = outunits + LSTs.units = 'hr' + LSTs.longname = 'Local Solar Time' + avgvalues.setAxis(0, comptime) + avgvalues.setAxis(1, modellats) + avgvalues.setAxis(2, modellons) + stdvalues.setAxis(0, comptime) + stdvalues.setAxis(1, modellats) + stdvalues.setAxis(2, modellons) + LSTs.setAxis(0, comptime) + LSTs.setAxis(1, modellats) + LSTs.setAxis(2, modellons) + avgoutfile = ('%s_%s_%s_%s-%s_diurnal_avg.nc') % (varbname, + dataname, monthname, + str(args.firstyear), str(args.lastyear)) + stdoutfile = ('%s_%s_%s_%s-%s_diurnal_std.nc') % (varbname, + dataname, monthname, str( + args.firstyear), + str(args.lastyear)) + LSToutfile = ('%s_%s_LocalSolarTimes.nc' % (varbname, dataname)) + if not os.path.exists(args.output_directory): + os.makedirs(args.output_directory) + f = cdms2.open( + os.path.join( + args.output_directory, + avgoutfile), + 'w') + g = cdms2.open( + os.path.join( + args.output_directory, + stdoutfile), + 'w') + h = cdms2.open( + os.path.join( + args.output_directory, + LSToutfile), + 'w') + f.write(avgvalues) + g.write(stdvalues) + h.write(LSTs) + f.close() + g.close() + h.close() + except Exception as err: + print "Failed for model %s with erro: %s" % (dataname, err) + + +print 'done' +args = P.parse_args(sys.argv[1:]) + +month = args.month +monthname = monthname_d[args.month] + +# -------------------------------------HARD-CODED INPUT (add to command line later?): + +# These models have been processed already (or tried and found wanting, +# e.g. problematic time coordinates): +skipMe = args.skip + +# Choose only one ensemble member per model, with the following ensemble-member code (for definitions, see +# http://cmip-pcmdi.llnl.gov/cmip5/docs/cmip5_data_reference_syntax.pdf): + + +# NOTE--These models do not supply 3hr data from the 'r1i1p1' ensemble member, +# but do supply it from other ensemble members: +# bcc-csm1-1 (3hr data is from r2i1p1) +# CCSM4 (3hr data is from r6i1p1) +# GFDL-CM3 (3hr data is from r2i1p1, r3i1p1, r4i1p1, r5i1p1) +# GISS-E2-H (3hr data is from r6i1p1, r6i1p3) +# GISS-E2-R (3hr data is from r6i1p2) + +varbname = 'pr' + +# Note that CMIP5 specifications designate (01:30, 04:30, 07:30, ..., 22:30) GMT for 3hr flux fields, but +# *WARNING* some GMT timepoints are actually (0, 3, 6,..., 21) in submitted CMIP5 data, despite character strings in +# file names (and time axis metadata) to the contrary. See CMIP5 documentation and errata! Overrides to +# correct these problems are given below: +# startGMT = '0:0:0.0' # Include 00Z as a possible starting time, to accomodate (0, 3, 6,..., 21)GMT in the input data. +# startime = -1.5 # Subtract 1.5h from (0, 3, 6,..., 21)GMT input data. This is needed for BNU-ESM, CCSM4 and CNRM-CM5. +# startime = -3.0 # Subtract 1.5h from (0, 3, 6,..., 21)GMT input +# data. This is needed for CMCC-CM. + +# ------------------------------------- + +nYears = args.lastyear - args.firstyear + 1 + +template = populateStringConstructor(args.filename_template, args) +template.variable = varbname + +print "TEMPLATE:", template() +fileList = glob.glob(os.path.join(args.modroot, template())) +print "FILES:", fileList +params = [INPUT(args, name, template) for name in fileList] +print "PARAMS:", params + +cdp.cdp_run.multiprocess(compute, params, num_workers=args.num_workers) diff --git a/src/python/diurnal/scripts/computeStdDailyMeansWrapped.py b/src/python/diurnal/scripts/computeStdDailyMeansWrapped.py new file mode 100755 index 000000000..35aa934c4 --- /dev/null +++ b/src/python/diurnal/scripts/computeStdDailyMeansWrapped.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python + +# From a high frequency (hourly, 3-hourly, ...) time series, compute daily means for one month over one or more years at +# each gridpoint, then their standard deviation. This version processes +# CMIP5 historical precipitation for years 1999-2005. + +# Charles Doutriaux September 2017 +# Curt Covey July 2017 +# (from ./old_computeDailyMeans.py) + +# This version has the PMP Parser "wrapped" around it, so it can be executed with input parameters in the command line +# ---> computeStdDailyMeansWrapped.py -i data -m7 --realization="r1i1p1" -t "sample_data_%(variable)_%(model).nc" + +import cdms2 +import cdtime +import genutil +import numpy.ma +import os +import glob +import sys +import cdp + +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor, INPUT + + +def compute(params): + fileName = params.fileName + startyear = params.args.firstyear + finalyear = params.args.lastyear + month = params.args.month + monthname = params.monthname + varbname = params.varname + template = populateStringConstructor(args.filename_template, args) + template.variable = varbname + + reverted = template.reverse(os.path.basename(fileName)) + dataname = reverted["model"] + if dataname not in args.skip: + try: + print 'Data source:', dataname + print 'Opening %s ...' % fileName + f = cdms2.open(fileName) + iYear = 0 + dmean = None + for year in range(startyear, finalyear + 1): + print 'Year %s:' % year + startTime = cdtime.comptime(year, month) + # Last possible second to get all tpoints + finishtime = startTime.add( + 1, cdtime.Month).add(-1, cdtime.Minute) + print 'Reading %s from %s for time interval %s to %s ...' % (varbname, fileName, startTime, finishtime) + # Transient variable stores data for current year's month. + tvarb = f(varbname, time=(startTime, finishtime, "ccn")) + # *HARD-CODES conversion from kg/m2/sec to mm/day. + tvarb *= 86400 + # The following tasks need to be done only once, extracting + # metadata from first-year file: + tc = tvarb.getTime().asComponentTime() + current = tc[0] + while current.month == month: + end = cdtime.comptime( + current.year, + current.month, + current.day).add( + 1, + cdtime.Day) + sub = tvarb(time=(current, end, "con")) + # Assumes first dimension of input ("axis#0") is time + tmp = numpy.ma.average(sub, axis=0) + sh = list(tmp.shape) + sh.insert(0, 1) + if dmean is None: + dmean = tmp.reshape(sh) + else: + dmean = numpy.ma.concatenate( + (dmean, tmp.reshape(sh)), axis=0) + current = end + iYear += 1 + f.close() + stdvalues = cdms2.MV2.array(genutil.statistics.std(dmean)) + stdvalues.setAxis(0, tvarb.getLatitude()) + stdvalues.setAxis(1, tvarb.getLongitude()) + stdvalues.id = 'dailySD' + # Standard deviation has same units as mean. + stdvalues.units = "mm/d" + stdoutfile = ('%s_%s_%s_%s-%s_std_of_dailymeans.nc') % (varbname, dataname, + monthname, str(startyear), str(finalyear)) + except Exception as err: + print "Failed for model: %s with error: %s" % (dataname, err) + if not os.path.exists(args.output_directory): + os.makedirs(args.output_directory) + g = cdms2.open(os.path.join(args.output_directory, stdoutfile), 'w') + g.write(stdvalues) + g.close() + + +args = P.parse_args(sys.argv[1:]) +month = args.month +startyear = args.firstyear +finalyear = args.lastyear +directory = args.modroot # Input directory for model data +# These models have been processed already (or tried and found wanting, +# e.g. problematic time coordinates): +skipMe = args.skip +print "SKIPPING:", skipMe + +# Choose only one ensemble member per model, with the following ensemble-member code (for definitions, see +# http://cmip-pcmdi.llnl.gov/cmip5/docs/cmip5_data_reference_syntax.pdf): + +# NOTE--These models do not supply 3hr data from the 'r1i1p1' ensemble member, +# but do supply it from other ensemble members: +# bcc-csm1-1 (3hr data is from r2i1p1) +# CCSM4 (3hr data is from r6i1p1) +# GFDL-CM3 (3hr data is from r2i1p1, r3i1p1, r4i1p1, r5i1p1) +# GISS-E2-H (3hr data is from r6i1p1, r6i1p3) +# GISS-E2-R (3hr data is from r6i1p2) + +varbname = "pr" + +# Note that CMIP5 specifications designate (01:30, 04:30, 07:30, ..., 22:30) GMT for 3hr flux fields, but +# *WARNING* some GMT timepoints are actually (0, 3, 6,..., 21) in submitted CMIP5 data, despite character strings in +# file names (and time axis metadata) to the contrary. See CMIP5 documentation and errata! Overrides to +# correct these problems are given below: +# Include 00Z as a possible starting time, to accomodate (0, 3, 6,..., +# 21)GMT in the input data. +# startime = -1.5 # Subtract 1.5h from (0, 3, 6,..., 21)GMT input +# data. This is needed for BNU-ESM, CCSM4 and CNRM-CM5. +# Subtract 1.5h from (0, 3, 6,..., 21)GMT input data. This is needed for +# CMCC-CM. + +# ------------------------------------- + +monthname = monthname_d[month] +nYears = finalyear - startyear + 1 +# Character strings for starting and ending day/GMT (*HARD-CODES +# particular GMT timepoints*): +# *WARNING* GMT timepoints are actually (0, 3, 6,..., 21) in the original TRMM/Obs4MIPs data, despite character strings +# in file names (and time axis metadata). See CMIP5 documentation and +# errata! + + +template = populateStringConstructor(args.filename_template, args) +template.variable = varbname + +fileList = glob.glob(os.path.join(directory, template())) +print "FILES:", fileList + +params = [INPUT(args, name, template) for name in fileList] +print "PARAMS:", params + +cdp.cdp_run.multiprocess(compute, params, num_workers=args.num_workers) diff --git a/src/python/diurnal/scripts/fourierDiurnalAllGridWrapped.py b/src/python/diurnal/scripts/fourierDiurnalAllGridWrapped.py new file mode 100755 index 000000000..2646d91a8 --- /dev/null +++ b/src/python/diurnal/scripts/fourierDiurnalAllGridWrapped.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python + +# This modifiction of ./fourierDiurnalAllGrid.py has the PMP Parser "wrapped" around it, so it's executed +# with input parameters in the Unix command line. For example: + +# ---> python fourierDiurnalAllGridWrapped.py -m7 + +# Charles Doutriaux September 2017 +# Curt Covey January 2017 + +# ------------------------------------------------------------------------- + +import cdms2 +import MV2 +from pcmdi_metrics.diurnal.fourierFFT import fastAllGridFT +import sys +import glob +import os + +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor + +P.add_argument("-t", "--filename_template", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_diurnal_avg.nc", + help="template for file names containing diurnal average") +P.add_argument("--model", default="*") +P.add_argument("--filename_template_LST", + default="pr_%(model)_LocalSolarTimes.nc", + help="template for file names point to Local Solar Time Files") + +args = P.parse_args(sys.argv[1:]) +month = args.month +monthname = monthname_d[month] +startyear = args.firstyear +finalyear = args.lastyear +yearrange = "%s-%s" % (startyear, finalyear) + +template = populateStringConstructor(args.filename_template, args) +template.month = monthname +template_LST = populateStringConstructor(args.filename_template_LST, args) +template_LST.month = monthname + +LSTfiles = glob.glob(os.path.join(args.modroot, template_LST())) +print "LSTFILES:", LSTfiles +print "TMPL", template_LST() +for LSTfile in LSTfiles: + print 'Reading %s ...' % LSTfile, os.path.basename(LSTfile) + reverted = template_LST.reverse(os.path.basename(LSTfile)) + model = reverted["model"] + print '====================' + print model + print '====================' + template.model = model + avgfile = template() + print 'Reading time series of mean diurnal cycle ...' + f = cdms2.open(LSTfile) + g = cdms2.open(os.path.join(args.modroot, avgfile)) + LSTs = f('LST') + avgs = g('diurnalmean') + print 'Input shapes: ', LSTs.shape, avgs.shape + + print 'Getting latitude and longitude coordinates.' + # Any file with grid info will do, so use Local Standard Times file: + modellats = LSTs.getLatitude() + modellons = LSTs.getLongitude() + + f.close() + g.close() + + print 'Taking fast Fourier transform of the mean diurnal cycle ...' + cycmean, maxvalue, tmax = fastAllGridFT(avgs, LSTs) + print ' Output:' + print ' cycmean', cycmean.shape + print ' maxvalue', maxvalue.shape + print ' tmax', tmax.shape + + print '"Re-decorating" Fourier harmonics with grid info, etc., ...' + cycmean = MV2.array(cycmean) + maxvalue = MV2.array(maxvalue) + tmax = MV2.array(tmax) + + cycmean.setAxis(0, modellats) + cycmean.setAxis(1, modellons) + cycmean.id = 'tmean' + cycmean.units = 'mm / day' + + maxvalue.setAxis(1, modellats) + maxvalue.setAxis(2, modellons) + maxvalue.id = 'S' + maxvalue.units = 'mm / day' + + tmax.setAxis(1, modellats) + tmax.setAxis(2, modellons) + tmax.id = 'tS' + tmax.units = 'GMT' + + print '... and writing to netCDF.' + f = cdms2.open( + os.path.join( + args.output_directory, + 'pr_' + + model + + '_' + + monthname + + '_' + + yearrange + + '_tmean.nc'), + 'w') + g = cdms2.open( + os.path.join( + args.output_directory, + 'pr_' + + model + + '_' + + monthname + + '_' + + yearrange + + '_S.nc'), + 'w') + h = cdms2.open( + os.path.join( + args.output_directory, + 'pr_' + + model + + '_' + + monthname + + '_' + + yearrange + + '_tS.nc'), + 'w') + f.write(cycmean) + g.write(maxvalue) + h.write(tmax) + f.close() + g.close() + h.close() diff --git a/src/python/diurnal/scripts/fourierDiurnalGridpoints.py b/src/python/diurnal/scripts/fourierDiurnalGridpoints.py new file mode 100755 index 000000000..a7e604e63 --- /dev/null +++ b/src/python/diurnal/scripts/fourierDiurnalGridpoints.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python + +# For a few selected gridpoints, read previously computed composite-mean and standard deviation of the diurnal +# cycle, compute Fourier harmonics and write output in a form suitable for Mathematica. The Fourier output of +# this script should match the output fields of +# ./fourierDiurnalAllGrid*.py at the selected gridpoints + +# Curt Covey June 2017 + +# (from ~/CMIP5/Tides/OtherFields/Models/CMCC-CM_etal/old_fourierDiurnalGridpoints.py) +# ------------------------------------------------------------------------- + + +from __future__ import print_function +import cdms2 +import MV2 +import sys +import glob +import os + +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor +from pcmdi_metrics.diurnal.fourierFFT import fastFT + +P.add_argument("-t", "--filename_template", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_diurnal_avg.nc", + help="template for file names containing diurnal average") +P.add_argument("--model", default="*") +P.add_argument("--filename_template_LST", + default="pr_%(model)_LocalSolarTimes.nc", + help="template for file names point to Local Solar Time Files") +P.add_argument("--filename_template_std", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_diurnal_std.nc", + help="template for file names containing diurnal std") +P.add_argument( + "-l", + "--lats", + nargs="*", + default=[ + 31.125, + 31.125, + 36.4, + 5.125, + 45.125, + 45.125], + help="latitudes") +P.add_argument("-L", "--lons", nargs="*", + default=[-83.125, 111.145, -97.5, 147.145, -169.145, -35.145], help="longitudes") +P.add_argument("-A", "--outnameasc", + type=str, + dest='outnameasc', + default='pr_%(month)_%(firstyear)-%(lastyear)_fourierDiurnalGridPoints.asc', + help="Output name for ascs") +args = P.parse_args(sys.argv[1:]) +month = args.month +monthname = monthname_d[month] +startyear = args.firstyear +finalyear = args.lastyear +yearrange = "%s-%s" % (startyear, finalyear) + +template = populateStringConstructor(args.filename_template, args) +template.month = monthname +template_std = populateStringConstructor(args.filename_template_std, args) +template_std.month = monthname +template_LST = populateStringConstructor(args.filename_template_LST, args) +template_LST.month = monthname + +LSTfiles = glob.glob(os.path.join(args.modroot, template_LST())) +print("LSTFILES:", LSTfiles) +print("TMPL", template_LST()) + +ascFile = populateStringConstructor(args.outnameasc, args) +ascFile.month = monthname +ascname = os.path.join(os.path.abspath(args.output_directory), ascFile()) + +if not os.path.exists(os.path.dirname(ascname)): + os.makedirs(os.path.dirname(ascname)) +fasc = open(ascname, "w") + + +gridptlats = [float(x) for x in args.lats] +gridptlons = [float(x) for x in args.lons] +nGridPoints = len(gridptlats) +assert(len(gridptlons) == nGridPoints) + +# gridptlats = [-29.125, -5.125, 45.125, 45.125] +# gridptlons = [-57.125, 75.125, -169.145, -35.145] +# Gridpoints for JULY samples in Figure 4 of Covey et al., JClimate 29: 4461 (2016): +# nGridPoints = 6 +# gridptlats = [ 31.125, 31.125, 36.4, 5.125, 45.125, 45.125] +# gridptlons = [-83.125, 111.145, -97.5, 147.145, -169.145, -35.145] + +N = 8 # Number of timepoints in a 24-hour cycle +for LSTfile in LSTfiles: + print('Reading %s ...' % LSTfile, os.path.basename(LSTfile), file=fasc) + print('Reading %s ...' % LSTfile, os.path.basename(LSTfile), file=fasc) + reverted = template_LST.reverse(os.path.basename(LSTfile)) + model = reverted["model"] + print('====================', file=fasc) + print(model, file=fasc) + print('====================', file=fasc) + template.model = model + avgfile = template() + template_std.model = model + stdfile = template_std() + print('Reading time series of mean diurnal cycle ...', file=fasc) + f = cdms2.open(LSTfile) + g = cdms2.open(os.path.join(args.modroot, avgfile)) + h = cdms2.open(os.path.join(args.modroot, stdfile)) + LSTs = f('LST') + print('Input shapes: ', LSTs.shape, file=fasc) + + modellats = LSTs.getLatitude() + modellons = LSTs.getLongitude() + latbounds = modellats.getBounds() + lonbounds = modellons.getBounds() + + # Gridpoints selected above may be offset slightly from points in full + # grid ... + closestlats = MV2.zeros(nGridPoints) + closestlons = MV2.zeros(nGridPoints) + pointLSTs = MV2.zeros((nGridPoints, N)) + avgvalues = MV2.zeros((nGridPoints, N)) + stdvalues = MV2.zeros((nGridPoints, N)) + # ... in which case, just pick the closest full-grid point: + for i in range(nGridPoints): + print(' (lat, lon) = (%8.3f, %8.3f)' % + (gridptlats[i], gridptlons[i]), file=fasc) + closestlats[i] = gridptlats[i] + closestlons[i] = gridptlons[i] % 360 + print(' Closest (lat, lon) for gridpoint = (%8.3f, %8.3f)' % + (closestlats[i], closestlons[i]), file=fasc) + # Time series for selected grid point: + avgvalues[i] = g( + 'diurnalmean', lat=( + closestlats[i], closestlats[i], "cob"), lon=( + closestlons[i], closestlons[i], "cob"), squeeze=1) + stdvalues[i] = h( + 'diurnalstd', lat=( + closestlats[i], closestlats[i], "cob"), lon=( + closestlons[i], closestlons[i], "cob"), squeeze=1) + pointLSTs[i] = f( + 'LST', lat=( + closestlats[i], closestlats[i], "cob"), lon=( + closestlons[i], closestlons[i], "cob"), squeeze=1) + print(' ', file=fasc) + f.close() + g.close() + h.close() + # Print results for input to Mathematica. + if monthname == 'Jan': + # In printed output, numbers for January data follow 0-5 for July data, + # hence begin with 6. + deltaI = 6 + else: + deltaI = 0 + prefix = args.modroot + for i in range(nGridPoints): + print( + 'For gridpoint %d at %5.1f deg latitude, %6.1f deg longitude ...' % + (i, gridptlats[i], gridptlons[i]), file=fasc) + print(' Local Solar Times are:', file=fasc) + print((prefix + 'LST%d = {') % (i + deltaI), file=fasc) + print(N * '%5.3f, ' % tuple(pointLSTs[i]), end="", file=fasc) + print('};', file=fasc) + print(' Mean values for each time-of-day are:', file=fasc) + print((prefix + 'mean%d = {') % (i + deltaI), file=fasc) + print(N * '%5.3f, ' % tuple(avgvalues[i]), end="", file=fasc) + print('};', file=fasc) + print(' Standard deviations for each time-of-day are:', file=fasc) + print((prefix + 'std%d = {') % (i + deltaI), file=fasc) + print(N * '%6.4f, ' % tuple(stdvalues[i]), end="", file=fasc) + print('};', file=fasc) + print(' ', file=fasc) + + # Take fast Fourier transform of the overall multi-year mean diurnal cycle. + print('************** ', avgvalues[0][0], file=fasc) + cycmean, maxvalue, tmax = fastFT(avgvalues, pointLSTs) + print('************** ', avgvalues[0][0], file=fasc) + # Print Fourier harmonics: + for i in range(nGridPoints): + print( + 'For gridpoint %d at %5.1f deg latitude, %6.1f deg longitude ...' % + (i, gridptlats[i], gridptlons[i]), file=fasc) + print(' Mean value over cycle = %6.2f' % cycmean[i], file=fasc) + print(' Diurnal maximum = %6.2f at %6.2f hr Local Solar Time.' % + (maxvalue[i, 0], tmax[i, 0] % 24), file=fasc) + print(' Semidiurnal maximum = %6.2f at %6.2f hr Local Solar Time.' % + (maxvalue[i, 1], tmax[i, 1] % 24), file=fasc) + print(' Terdiurnal maximum = %6.2f at %6.2f hr Local Solar Time.' % + (maxvalue[i, 2], tmax[i, 2] % 24), file=fasc) + +print("Results sent to:", ascname) diff --git a/src/python/diurnal/scripts/savg_fourierWrappedInOut.py b/src/python/diurnal/scripts/savg_fourierWrappedInOut.py new file mode 100755 index 000000000..8a14be16a --- /dev/null +++ b/src/python/diurnal/scripts/savg_fourierWrappedInOut.py @@ -0,0 +1,291 @@ +#!/usr/bin/env python + +# Spatially vector-average pairs of lat/lon fields that store Fourier (amplitude, phase) values for the diurnal +# cycle of preciprecipitation. This version does not explicitly mask out areas where the diurnal cycle is weak, +# but such areas should count for little since the vector amplitudes there +# are small. + + +# This modifiction of ./savg_fourier.py has the PMP Parser "wrapped" around it, so it's executed with both input and +# output parameters specified in the Unix command line. +# ... but this version is designed to get land-sea masks for any model, not just MIROC5 and CCSM4. + +# Charles Doutriaux September 2017 +# Curt Covey (from ./savg_fourierWrappedInOutCCSMandMIROC.py) +# June 2017 + +import cdms2 +import cdutil +import MV2 +import os +import sys +import glob +import pcmdi_metrics +import collections +import json +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor + + +P.add_argument("-j", "--outnamejson", + type=str, + dest='outnamejson', + default='pr_%(month)_%(firstyear)-%(lastyear)_savg_DiurnalFourier.json', + help="Output name for jsons") + +P.add_argument("--lat1", type=float, default=-50., help="First latitude") +P.add_argument("--lat2", type=float, default=50., help="Last latitude") +P.add_argument("--lon1", type=float, default=0., help="First longitude") +P.add_argument("--lon2", type=float, default=360., help="Last longitude") +P.add_argument("--region_name", type=str, default="TRMM", + help="name for the region of interest") + +P.add_argument("-t", "--filename_template", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_S.nc", + help="template for getting at amplitude files") +P.add_argument("--filename_template_tS", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_tS.nc", + help="template for phase files") +P.add_argument("--filename_template_sftlf", + default="cmip5.%(model).%(experiment).r0i0p0.fx.atm.fx.sftlf.%(version).latestX.xml", + help="template for sftlf file names") +P.add_argument("--model", default="*") + +args = P.parse_args(sys.argv[1:]) +month = args.month +monthname = monthname_d[month] +startyear = args.firstyear +finalyear = args.lastyear +years = "%s-%s" % (startyear, finalyear) + + +print 'Specifying latitude / longitude domain of interest ...' +# TRMM (observed) domain: +latrange = (args.lat1, args.lat2) +lonrange = (args.lon1, args.lon2) + +region = cdutil.region.domain(latitude=latrange, longitude=lonrange) + +if args.region_name == "": + region_name = "{:g}_{:g}&{:g}_{:g}".format(*(latrange + lonrange)) +else: + region_name = args.region_name + +# Amazon basin: +# latrange = (-15.0, -5.0) +# lonrange = (285.0, 295.0) + +# Functions to convert phase between angle-in-radians and hours, for +# either a 12- or 24-hour clock, i.e. for clocktype = 12 or 24: + + +def hrs_to_rad(hours, clocktype): + import MV2 + return 2 * MV2.pi * hours / clocktype + + +def rad_to_hrs(phase, clocktype): + import MV2 + return phase * clocktype / 2 / MV2.pi + + +def vectoravg(hr1, hr2, clocktype): + 'Function to test vector-averaging of two time values:' + import MV2 + + sin_avg = (MV2.sin(hrs_to_rad(hr1, clocktype)) + + MV2.sin(hrs_to_rad(hr2, clocktype))) / 2 + cos_avg = (MV2.cos(hrs_to_rad(hr1, clocktype)) + + MV2.cos(hrs_to_rad(hr2, clocktype))) / 2 + return rad_to_hrs(MV2.arctan2(sin_avg, cos_avg), clocktype) + + +def spacevavg(tvarb1, tvarb2, sftlf, model): + ''' + Given a "root filename" and month/year specifications, vector-average lat/lon arrays in an (amplitude, phase) + pair of input data files. Each input data file contains diurnal (24h), semidiurnal (12h) and terdiurnal (8h) + Fourier harmonic components of the composite mean day/night cycle. + + Vector-averaging means we consider the input data to be readings on an 8-, 12- or 24-hour clock and separately + average the Cartesian components (called "cosine" and "sine" below). Then the averaged components are combined + back into amplitude and phase values and returned. + + Space-averaging is done globally, as well as separately for land and ocean areas. + ''' + + glolf = cdutil.averager(sftlf, axis='xy') + print ' Global mean land fraction = %5.3f' % glolf + outD = {} # Output dictionary to be returned by this function + harmonics = [1, 2, 3] + for harmonic in harmonics: + ampl = tvarb1[harmonic - 1] + tmax = tvarb2[harmonic - 1] + # print ampl[:, :] + # print tmax[:, :] + clocktype = 24 / harmonic + cosine = MV2.cos(hrs_to_rad(tmax, clocktype)) * ampl # X-component + sine = MV2.sin(hrs_to_rad(tmax, clocktype)) * ampl # Y-component + + print 'Area-averaging globally, over land only, and over ocean only ...' + # Average Cartesian components ... + cos_avg_glo = cdutil.averager(cosine, axis='xy') + sin_avg_glo = cdutil.averager(sine, axis='xy') + cos_avg_lnd = cdutil.averager(cosine * sftlf, axis='xy') + sin_avg_lnd = cdutil.averager(sine * sftlf, axis='xy') + cos_avg_ocn = cos_avg_glo - cos_avg_lnd + sin_avg_ocn = sin_avg_glo - sin_avg_lnd + # ... normalized by land-sea fraction: + cos_avg_lnd /= glolf + sin_avg_lnd /= glolf + cos_avg_ocn /= (1 - glolf) + sin_avg_ocn /= (1 - glolf) + # Amplitude and phase: + # * 86400 Convert kg/m2/s -> mm/d? + amp_avg_glo = MV2.sqrt(sin_avg_glo**2 + cos_avg_glo**2) + # * 86400 Convert kg/m2/s -> mm/d? + amp_avg_lnd = MV2.sqrt(sin_avg_lnd**2 + cos_avg_lnd**2) + # * 86400 Convert kg/m2/s -> mm/d? + amp_avg_ocn = MV2.sqrt(sin_avg_ocn**2 + cos_avg_ocn**2) + pha_avg_glo = MV2.remainder( + rad_to_hrs( + MV2.arctan2( + sin_avg_glo, + cos_avg_glo), + clocktype), + clocktype) + pha_avg_lnd = MV2.remainder( + rad_to_hrs( + MV2.arctan2( + sin_avg_lnd, + cos_avg_lnd), + clocktype), + clocktype) + pha_avg_ocn = MV2.remainder( + rad_to_hrs( + MV2.arctan2( + sin_avg_ocn, + cos_avg_ocn), + clocktype), + clocktype) + if 'CMCC-CM' in model: + # print '** Correcting erroneous time recording in ', rootfname + pha_avg_lnd -= 3.0 + pha_avg_lnd = MV2.remainder(pha_avg_lnd, clocktype) + elif 'BNU-ESM' in model or 'CCSM4' in model or 'CNRM-CM5' in model: + # print '** Correcting erroneous time recording in ', rootfname + pha_avg_lnd -= 1.5 + pha_avg_lnd = MV2.remainder(pha_avg_lnd, clocktype) + print 'Converting singleton transient variables to plain floating-point numbers ...' + amp_avg_glo = float(amp_avg_glo) + pha_avg_glo = float(pha_avg_glo) + amp_avg_lnd = float(amp_avg_lnd) + pha_avg_lnd = float(pha_avg_lnd) + amp_avg_ocn = float(amp_avg_ocn) + pha_avg_ocn = float(pha_avg_ocn) + print '%s %s-harmonic amplitude, phase = %7.3f mm/d, %7.3f hrsLST averaged globally' \ + % (monthname, harmonic, amp_avg_glo, pha_avg_glo) + print '%s %s-harmonic amplitude, phase = %7.3f mm/d, %7.3f hrsLST averaged over land'\ + % (monthname, harmonic, amp_avg_lnd, pha_avg_lnd) + print '%s %s-harmonic amplitude, phase = %7.3f mm/d, %7.3f hrsLST averaged over ocean'\ + % (monthname, harmonic, amp_avg_ocn, pha_avg_ocn) + # Sub-dictionaries, one for each harmonic component: + outD['harmonic' + str(harmonic)] = {} + outD['harmonic' + str(harmonic)]['amp_avg_lnd'] = amp_avg_lnd + outD['harmonic' + str(harmonic)]['pha_avg_lnd'] = pha_avg_lnd + outD['harmonic' + str(harmonic)]['amp_avg_ocn'] = amp_avg_ocn + outD['harmonic' + str(harmonic)]['pha_avg_ocn'] = pha_avg_ocn + return outD + + +print 'Preparing to write output to JSON file ...' +if not os.path.exists(args.output_directory): + os.makedirs(args.output_directory) +jsonFile = populateStringConstructor(args.outnamejson, args) +jsonFile.month = monthname + +jsonname = os.path.join(os.path.abspath(args.output_directory), jsonFile()) + +if not os.path.exists(jsonname) or args.append is False: + print 'Initializing dictionary of statistical results ...' + stats_dic = {} + metrics_dictionary = collections.OrderedDict() +else: + with open(jsonname) as f: + metrics_dictionary = json.load(f) + stats_dic = metrics_dictionary["RESULTS"] + +OUT = pcmdi_metrics.io.base.Base( + os.path.abspath( + args.output_directory), + os.path.basename(jsonname)) +disclaimer = open( + os.path.join( + sys.prefix, + "share", + "pmp", + "disclaimer.txt")).read() +metrics_dictionary["DISCLAIMER"] = disclaimer +metrics_dictionary["REFERENCE"] = "The statistics in this file are based on Covey et al., J Climate 2016" + +# Accumulate output from each model (or observed) data source in the +# Python dictionary. +template_S = populateStringConstructor(args.filename_template, args) +template_S.month = monthname +template_tS = populateStringConstructor(args.filename_template_tS, args) +template_tS.month = monthname +template_sftlf = populateStringConstructor(args.filename_template_sftlf, args) +template_sftlf.month = monthname + +print "TEMPLATE:", template_S() +files_S = glob.glob(os.path.join(args.modroot, template_S())) +print files_S +for file_S in files_S: + print 'Reading Amplitude from %s ...' % file_S + reverted = template_S.reverse(os.path.basename(file_S)) + model = reverted["model"] + try: + template_tS.model = model + template_sftlf.model = model + S = cdms2.open(file_S)("S", region) + print 'Reading Phase from %s ...' % os.path.join(args.modroot, template_tS()) + tS = cdms2.open( + os.path.join( + args.modroot, + template_tS()))( + "tS", + region) + print 'Reading sftlf from %s ...' % os.path.join(args.modroot, template_sftlf()) + try: + sftlf_fnm = glob.glob( + os.path.join( + args.modroot, + template_sftlf()))[0] + sftlf = cdms2.open(sftlf_fnm)("sftlf", region) / 100. + except BaseException as err: + print 'Failed reading sftlf from file (error was: %s)' % err + print 'Creating one for you' + sftlf = cdutil.generateLandSeaMask(S.getGrid()) + + if model not in stats_dic: + stats_dic[model] = {region_name: spacevavg(S, tS, sftlf, model)} + else: + stats_dic[model].update( + {region_name: spacevavg(S, tS, sftlf, model)}) + print stats_dic + except Exception as err: + print "Failed for model %s with error %s" % (model, err) + +# Write output to JSON file. +metrics_dictionary["RESULTS"] = stats_dic +rgmsk = metrics_dictionary.get("RegionalMasking", {}) +nm = region_name +region.id = nm +rgmsk[nm] = {"id": nm, "domain": region} +metrics_dictionary["RegionalMasking"] = rgmsk +OUT.write( + metrics_dictionary, + json_structure=["model", "domain", "harmonic", "statistic"], + indent=4, + separators=( + ',', + ': ')) +print 'done' diff --git a/src/python/diurnal/scripts/std_of_dailymeansWrappedInOut.py b/src/python/diurnal/scripts/std_of_dailymeansWrappedInOut.py new file mode 100755 index 000000000..54830de7b --- /dev/null +++ b/src/python/diurnal/scripts/std_of_dailymeansWrappedInOut.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python + +# Globally average the variance of daily mean precipitation values from a given month spanning a given set of years, +# then take the square root. This is the correct order of operations to get the fraction of total variance, as shown +# by Covey and Gehne 2016 (https://e-reports-ext.llnl.gov/pdf/823104.pdf). + +# For Januarys, taking the square root before globally averaging would be close to the DJF results mapped in the +# middle-right panel of Figure 2 in Trenberth, Zhang & Gehne (2107) "Intermittency in Precipitation," J. Hydromet. +# This is done in ./plot_std_of_dailymeans.py.) + +# This has the PMP Parser "wrapped" around it, so it's executed with both input and output parameters specified in the +# Unix command line. For example: + +# ---> python std_of_dailymeansWrappedInOut.py -i out -o jsons -m 7 + +# Charles Doutriaux September 2017 +# Curt Covey (from ./old_std_of_dailymeansWrappedInOut.py) +# July 2017 + +import cdms2 +import cdutil +import os +import numpy.ma +import pcmdi_metrics +import collections +import glob +import sys +import json + +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor, INPUT +import cdp + + +def compute(param): + template = populateStringConstructor(args.filename_template, args) + template.variable = param.varname + template.month = param.monthname + fnameRoot = param.fileName + reverted = template.reverse(os.path.basename(fnameRoot)) + model = reverted["model"] + print 'Specifying latitude / longitude domain of interest ...' + datanameID = 'dailySD' # Short ID name of output data + latrange = (param.args.lat1, param.args.lat2) + lonrange = (param.args.lon1, param.args.lon2) + region = cdutil.region.domain(latitude=latrange, longitude=lonrange) + if param.args.region_name == "": + region_name = "{:g}_{:g}&{:g}_{:g}".format(*(latrange + lonrange)) + else: + region_name = param.args.region_name + print 'Reading {} ...'.format(fnameRoot) + try: + f = cdms2.open(fnameRoot) + x = f(datanameID, region) + units = x.units + print ' Shape =', x.shape + print 'Finding RMS area-average ...' + x = x * x + x = cdutil.averager(x, axis='xy') + x = numpy.ma.sqrt(x) + print 'For %8s in %s, average variance of daily values = (%5.2f %s)^2' % (model, monthname, x, units) + f.close() + except Exception as err: + print "Failed for model %s with error %s" % (model, err) + x = 1.e20 + return model, region, {region_name: x} + + +P.add_argument("-j", "--outnamejson", + type=str, + dest='outnamejson', + default='pr_%(month)_%(firstyear)_%(lastyear)_std_of_dailymeans.json', + help="Output name for jsons") + +P.add_argument("--lat1", type=float, default=-50., help="First latitude") +P.add_argument("--lat2", type=float, default=50., help="Last latitude") +P.add_argument("--lon1", type=float, default=0., help="First longitude") +P.add_argument("--lon2", type=float, default=360., help="Last longitude") +P.add_argument("--region_name", type=str, default="TRMM", + help="name for the region of interest") + +P.add_argument("-t", "--filename_template", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_std_of_dailymeans.nc") +P.add_argument("--model", default="*") + +args = P.parse_args(sys.argv[1:]) +month = args.month +monthname = monthname_d[month] +startyear = args.firstyear +finalyear = args.lastyear + +template = populateStringConstructor(args.filename_template, args) +template.month = monthname + +print "TEMPLATE NAME:", template() + + +print 'Preparing to write output to JSON file ...' +if not os.path.exists(args.output_directory): + os.makedirs(args.output_directory) +jsonFile = populateStringConstructor(args.outnamejson, args) +jsonFile.month = monthname + +jsonname = os.path.join(os.path.abspath(args.output_directory), jsonFile()) + +if not os.path.exists(jsonname) or args.append is False: + print 'Initializing dictionary of statistical results ...' + stats_dic = {} + metrics_dictionary = collections.OrderedDict() +else: + with open(jsonname) as f: + metrics_dictionary = json.load(f) + stats_dic = metrics_dictionary["RESULTS"] + +OUT = pcmdi_metrics.io.base.Base( + os.path.abspath( + args.output_directory), + jsonFile()) +disclaimer = open( + os.path.join( + sys.prefix, + "share", + "pmp", + "disclaimer.txt")).read() +metrics_dictionary["DISCLAIMER"] = disclaimer +metrics_dictionary["REFERENCE"] = "The statistics in this file are based on Trenberth, Zhang & Gehne, J Hydromet. 2017" + + +files = glob.glob(os.path.join(args.modroot, template())) +print files +params = [INPUT(args, name, template) for name in files] +print "PARAMS:", params + +results = cdp.cdp_run.multiprocess( + compute, params, num_workers=args.num_workers) + +for r in results: + m, region, res = r + if r[0] not in stats_dic: + stats_dic[m] = res + else: + stats_dic[m].update(res) + +print 'Writing output to JSON file ...', stats_dic +metrics_dictionary["RESULTS"] = stats_dic +rgmsk = metrics_dictionary.get("RegionalMasking", {}) +print "REG MASK:", rgmsk +nm = res.keys()[0] +region.id = nm +rgmsk[nm] = {"id": nm, "domain": region} +metrics_dictionary["RegionalMasking"] = rgmsk +OUT.write( + metrics_dictionary, + json_structure=["model", "domain"], + indent=4, + separators=( + ',', + ': ')) +print 'done' diff --git a/src/python/diurnal/scripts/std_of_hourlyvaluesWrappedInOut.py b/src/python/diurnal/scripts/std_of_hourlyvaluesWrappedInOut.py new file mode 100755 index 000000000..e17101394 --- /dev/null +++ b/src/python/diurnal/scripts/std_of_hourlyvaluesWrappedInOut.py @@ -0,0 +1,170 @@ +#!/usr/bin/env python + +# Globally average the day-to-day variance of CMIP5/AMIP composite-diurnal-cycle precipitation--"Var(Pi)" in Trenberth, +# Zhang & Gehne 2107: "Intermittency in Precipitation," J. Hydrometeorology-from a given month spanning years 1999-2005. +# Square the input standard deviations before averaging over the hours of the day and then area, and then finally take +# the square root. This is the correct order of operations to get the fraction of total variance, as shown by Covey and +# Gehne 2016 (https://e-reports-ext.llnl.gov/pdf/823104.pdf). + +# This has the PMP Parser "wrapped" around it, so it's executed with both input and output parameters specified in the +# Unix command line. For example: + +# ---> python std_of_hourlyvaluesWrappedInOut.py -m7 + +# Curt Covey (from ./old_std_of_hourlyvaluesWrappedInOut.py) July 2017 + + +import cdms2 +import cdutil +import os +import numpy.ma +import pcmdi_metrics +import collections +import glob +import sys +import cdp +import json +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor, INPUT + + +def compute(param): + template = populateStringConstructor(args.filename_template, args) + template.variable = param.varname + template.month = param.monthname + fnameRoot = param.fileName + reverted = template.reverse(os.path.basename(fnameRoot)) + model = reverted["model"] + print 'Specifying latitude / longitude domain of interest ...' + datanameID = 'diurnalstd' # Short ID name of output data + latrange = (param.args.lat1, param.args.lat2) + lonrange = (param.args.lon1, param.args.lon2) + region = cdutil.region.domain(latitude=latrange, longitude=lonrange) + if param.args.region_name == "": + region_name = "{:g}_{:g}&{:g}_{:g}".format(*(latrange + lonrange)) + else: + region_name = param.args.region_name + print 'Reading %s ...' % fnameRoot + reverted = template.reverse(os.path.basename(fnameRoot)) + model = reverted["model"] + try: + f = cdms2.open(fnameRoot) + x = f(datanameID, region) + units = x.units + print ' Shape =', x.shape + print 'Finding RMS area-average ...' + x = x * x + x = cdutil.averager(x, weights='unweighted') + x = cdutil.averager(x, axis='xy') + x = numpy.ma.sqrt(x) + print 'For %8s in %s, average variance of hourly values = (%5.2f %s)^2' % (model, monthname, x, units) + f.close() + except Exception as err: + print "Failed model %s with error: %s" % (model, err) + x = 1.e20 + return model, region, {region_name: x} + + +P.add_argument("-j", "--outnamejson", + type=str, + dest='outnamejson', + default='pr_%(month)_%(firstyear)-%(lastyear)_std_of_hourlymeans.json', + help="Output name for jsons") + +P.add_argument("--lat1", type=float, default=-50., help="First latitude") +P.add_argument("--lat2", type=float, default=50., help="Last latitude") +P.add_argument("--lon1", type=float, default=0., help="First longitude") +P.add_argument("--lon2", type=float, default=360., help="Last longitude") +P.add_argument("--region_name", type=str, default="TRMM", + help="name for the region of interest") + +P.add_argument("-t", "--filename_template", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_diurnal_std.nc") +P.add_argument("--model", default="*") + +args = P.parse_args(sys.argv[1:]) +month = args.month +monthname = monthname_d[month] +startyear = args.firstyear +finalyear = args.lastyear + +template = populateStringConstructor(args.filename_template, args) +template.month = monthname + +print "TEMPLATE NAME:", template() + +print 'Specifying latitude / longitude domain of interest ...' +# TRMM (observed) domain: +latrange = (args.lat1, args.lat2) +lonrange = (args.lon1, args.lon2) + +region = cdutil.region.domain(latitude=latrange, longitude=lonrange) + +# Amazon basin: +# latrange = (-15.0, -5.0) +# lonrange = (285.0, 295.0) + + +print 'Preparing to write output to JSON file ...' +if not os.path.exists(args.output_directory): + os.makedirs(args.output_directory) +jsonFile = populateStringConstructor(args.outnamejson, args) +jsonFile.month = monthname + +jsonname = os.path.join(os.path.abspath(args.output_directory), jsonFile()) + +if not os.path.exists(jsonname) or args.append is False: + print 'Initializing dictionary of statistical results ...' + stats_dic = {} + metrics_dictionary = collections.OrderedDict() +else: + with open(jsonname) as f: + metrics_dictionary = json.load(f) + stats_dic = metrics_dictionary["RESULTS"] + +OUT = pcmdi_metrics.io.base.Base( + os.path.abspath( + args.output_directory), + jsonFile()) + +disclaimer = open( + os.path.join( + sys.prefix, + "share", + "pmp", + "disclaimer.txt")).read() +metrics_dictionary["DISCLAIMER"] = disclaimer +metrics_dictionary["REFERENCE"] = "The statistics in this file are based on Trenberth, Zhang & Gehne, J Hydromet. 2017" + + +files = glob.glob(os.path.join(args.modroot, template())) +print files + + +params = [INPUT(args, name, template) for name in files] +print "PARAMS:", params + +results = cdp.cdp_run.multiprocess( + compute, params, num_workers=args.num_workers) + +for r in results: + m, region, res = r + if r[0] not in stats_dic: + stats_dic[m] = res + else: + stats_dic[m].update(res) + +print 'Writing output to JSON file ...' +metrics_dictionary["RESULTS"] = stats_dic +rgmsk = metrics_dictionary.get("RegionalMasking", {}) +nm = res.keys()[0] +region.id = nm +rgmsk[nm] = {"id": nm, "domain": region} +metrics_dictionary["RegionalMasking"] = rgmsk +OUT.write( + metrics_dictionary, + json_structure=["model", "domain"], + indent=4, + separators=( + ',', + ': ')) +print 'done' diff --git a/src/python/diurnal/scripts/std_of_meandiurnalcycWrappedInOut.py b/src/python/diurnal/scripts/std_of_meandiurnalcycWrappedInOut.py new file mode 100755 index 000000000..eed34d30d --- /dev/null +++ b/src/python/diurnal/scripts/std_of_meandiurnalcycWrappedInOut.py @@ -0,0 +1,172 @@ +#!/usr/bin/env python + +# Globally average the variance of CMIP5/AMIP mean-diurnal-cycle precipitation from a given month and range of years. +# Square the standard deviatios before averaging over area, and then take the square root. This is the correct order of +# operations to get the fraction of total variance, per Covey and Gehne +# 2016 (https://e-reports-ext.llnl.gov/pdf/823104.pdf). + +# This has the PMP Parser "wrapped" around it, so it's executed with both input and output parameters specified in the +# Unix command line. + +# Charles Doutriaux September 2017 +# Curt Covey (from ./old_meandiurnalcycWrappedInOut.py) July 2017 + + +import cdms2 +import cdutil +import os +import pcmdi_metrics +import collections +import glob +import sys +import genutil +import json +import cdp +from pcmdi_metrics.diurnal.common import monthname_d, P, populateStringConstructor, INPUT + + +def compute(param): + template = populateStringConstructor(args.filename_template, args) + template.variable = param.varname + template.month = param.monthname + fnameRoot = param.fileName + reverted = template.reverse(os.path.basename(fnameRoot)) + model = reverted["model"] + print 'Specifying latitude / longitude domain of interest ...' + datanameID = 'diurnalmean' # Short ID name of output data + latrange = (param.args.lat1, param.args.lat2) + lonrange = (param.args.lon1, param.args.lon2) + region = cdutil.region.domain(latitude=latrange, longitude=lonrange) + if param.args.region_name == "": + region_name = "{:g}_{:g}&{:g}_{:g}".format(*(latrange + lonrange)) + else: + region_name = param.args.region_name + print 'Reading %s ...' % fnameRoot + try: + f = cdms2.open(fnameRoot) + x = f(datanameID, region) + units = x.units + print ' Shape =', x.shape + + print 'Finding standard deviation over first dimension (time of day) ...' + x = genutil.statistics.std(x) + print ' Shape =', x.shape + + print 'Finding r.m.s. average over 2nd-3rd dimensions (area) ...' + x = x * x + x = cdutil.averager(x, axis='xy') + x = cdms2.MV2.sqrt(x) + + print 'For %8s in %s, average variance of hourly values = (%5.2f %s)^2' % (model, monthname, x, units) + f.close() + except Exception as err: + print "Failed model %s with error" % (err) + x = 1.e20 + return model, region, {region_name: float(x)} + + +P.add_argument("-j", "--outnamejson", + type=str, + dest='outnamejson', + default='pr_%(month)_%(firstyear)-%(lastyear)_std_of_meandiurnalcyc.json', + help="Output name for jsons") + +P.add_argument("--lat1", type=float, default=-50., help="First latitude") +P.add_argument("--lat2", type=float, default=50., help="Last latitude") +P.add_argument("--lon1", type=float, default=0., help="First longitude") +P.add_argument("--lon2", type=float, default=360., help="Last longitude") +P.add_argument("--region_name", type=str, default="TRMM", + help="name for the region of interest") + +P.add_argument("-t", "--filename_template", + default="pr_%(model)_%(month)_%(firstyear)-%(lastyear)_diurnal_avg.nc") +P.add_argument("--model", default="*") + +args = P.parse_args(sys.argv[1:]) +month = args.month +monthname = monthname_d[month] +startyear = args.firstyear +finalyear = args.lastyear + +template = populateStringConstructor(args.filename_template, args) +template.month = monthname + +print "TEMPLATE NAME:", template() + +print 'Specifying latitude / longitude domain of interest ...' +# TRMM (observed) domain: +latrange = (args.lat1, args.lat2) +lonrange = (args.lon1, args.lon2) + +region = cdutil.region.domain(latitude=latrange, longitude=lonrange) + +# Amazon basin: +# latrange = (-15.0, -5.0) +# lonrange = (285.0, 295.0) + + +print 'Preparing to write output to JSON file ...' +if not os.path.exists(args.output_directory): + os.makedirs(args.output_directory) +jsonFile = populateStringConstructor(args.outnamejson, args) +jsonFile.month = monthname + +jsonname = os.path.join(os.path.abspath(args.output_directory), jsonFile()) + +if not os.path.exists(jsonname) or args.append is False: + print 'Initializing dictionary of statistical results ...' + stats_dic = {} + metrics_dictionary = collections.OrderedDict() +else: + with open(jsonname) as f: + metrics_dictionary = json.load(f) + print "LOADE WITH KEYS:", metrics_dictionary.keys() + stats_dic = metrics_dictionary["RESULTS"] + +OUT = pcmdi_metrics.io.base.Base( + os.path.abspath( + args.output_directory), + jsonFile()) +disclaimer = open( + os.path.join( + sys.prefix, + "share", + "pmp", + "disclaimer.txt")).read() +metrics_dictionary["DISCLAIMER"] = disclaimer +metrics_dictionary["REFERENCE"] = "The statistics in this file are based on Trenberth, Zhang & Gehne, J Hydromet. 2017" + + +files = glob.glob(os.path.join(args.modroot, template())) +print files + +params = [INPUT(args, name, template) for name in files] +print "PARAMS:", params + +results = cdp.cdp_run.multiprocess( + compute, params, num_workers=args.num_workers) + +for r in results: + m, region, res = r + if r[0] not in stats_dic: + stats_dic[m] = res + else: + stats_dic[m].update(res) + +print 'Writing output to JSON file ...' +metrics_dictionary["RESULTS"] = stats_dic +print "KEYS AT END:", metrics_dictionary.keys() +rgmsk = metrics_dictionary.get("RegionalMasking", {}) +print "REG MASK:", rgmsk +nm = res.keys()[0] +region.id = nm +rgmsk[nm] = {"id": nm, "domain": region} +metrics_dictionary["RegionalMasking"] = rgmsk +OUT.write( + metrics_dictionary, + json_structure=["model", "domain"], + indent=4, + separators=( + ',', + ': ')) +print 'done' diff --git a/tests/diurnal/results/json/pr_Jul_1999-2005_savg_DiurnalFourier.json b/tests/diurnal/results/json/pr_Jul_1999-2005_savg_DiurnalFourier.json new file mode 100644 index 000000000..2d303a7a5 --- /dev/null +++ b/tests/diurnal/results/json/pr_Jul_1999-2005_savg_DiurnalFourier.json @@ -0,0 +1,127 @@ +{ + "REFERENCE": "The statistics in this file are based on Covey et al., J Climate 2016", + "json_structure": [ + "model", + "domain", + "harmonic", + "statistic" + ], + "provenance": { + "platform": { + "OS": "Linux", + "Version": "4.13.0-19-generic", + "Name": "drdoom" + }, + "userId": "doutriaux1", + "osAccess": false, + "commandLine": "/home/doutriaux1/anaconda2/envs/nightly2/bin/savg_fourierWrappedInOut.py --append -i tests/diurnal/results/nc -o test_data/results/jsons -m7", + "date": "2017-12-21 17:20:06", + "conda": { + "Platform": "linux-64", + "Version": "4.3.30", + "IsPrivate": "False", + "envVersion": "4.3.30", + "buildVersion": "3.0.28", + "PythonVersion": "2.7.14.final.0", + "RootEnvironment": "/home/doutriaux1/anaconda2 (writable)", + "DefaultEnvironment": "/home/doutriaux1/anaconda2/envs/nightly2" + }, + "packages": { + "blas": "0.2.20", + "cdms": "2.12.2017.11.30.06.06.14b9b9380c14680498a0ba8e9aefb4c84d2756be", + "CDP": "1.1.0", + "cdtime": "2017.11.15", + "cdutil": "2.12.2017.11.16.21.46.e717f782e07d0a187f2b616a408103af20ceb57f", + "clapack": "3.2.1", + "lapack": "3.6.1", + "esmf": "7.0.0", + "esmpy": "7.0.0", + "genutil": "2017.11.15", + "python": "2.7.14", + "matplotlib": "2.1.1", + "mesalib": "17.2.0", + "numpy": "1.13.3", + "vcs": "2.12.2017.12.14.17.01.99491e457b26e4e5fbb875e5ba75e7265f5e5b1b", + "vtk": "7.1.0.2.12", + "PMP": "g6ca4d58", + "PMPObs": null + }, + "openGL": { + "vendor": "X.Org", + "renderer": "AMD TURKS (DRM 2.50.0 / 4.13.0-19-generic, LLVM 5.0.0)", + "version": "3.0 Mesa 17.2.2", + "shading language version": "1.30", + "GLX": { + "client": { + "vendor": "Mesa Project and SGI", + "version": "1.4" + }, + "version": "1.4", + "server": { + "vendor": "SGI", + "version": "1.4" + } + } + } + }, + "RESULTS": { + "CMCC": { + "TROPICS": { + "harmonic2": { + "pha_avg_lnd": 4.630479623191498, + "amp_avg_ocn": 0.1013535443192936, + "pha_avg_ocn": 6.216779605786119, + "amp_avg_lnd": 0.4446492127610328 + }, + "harmonic3": { + "pha_avg_lnd": 0.6424494812695326, + "amp_avg_ocn": 0.018649923583445832, + "pha_avg_ocn": 2.9074929649270276, + "amp_avg_lnd": 0.05872192253917548 + }, + "harmonic1": { + "pha_avg_lnd": 16.571398252926915, + "amp_avg_ocn": 0.36913808560762545, + "pha_avg_ocn": 8.261227664939977, + "amp_avg_lnd": 0.8902230200805754 + } + }, + "TRMM": { + "harmonic2": { + "pha_avg_lnd": 4.527929393121403, + "amp_avg_ocn": 0.08501471624047707, + "pha_avg_ocn": 6.358694740438092, + "amp_avg_lnd": 0.3772614118017368 + }, + "harmonic3": { + "pha_avg_lnd": 0.44504280439914284, + "amp_avg_ocn": 0.012034623053092585, + "pha_avg_ocn": 2.3086496727977948, + "amp_avg_lnd": 0.04654577252551838 + }, + "harmonic1": { + "pha_avg_lnd": 16.275695610410192, + "amp_avg_ocn": 0.31771654735858623, + "pha_avg_ocn": 7.969509128258713, + "amp_avg_lnd": 0.77801266028589 + } + } + } + }, + "RegionalMasking": { + "TROPICS": { + "domain": { + "TROPICS": "cdutil.region.domain(latitude=(-30.0, 30.0),longitude=(0.0, 360.0))" + }, + "id": "TROPICS" + }, + "TRMM": { + "domain": { + "TRMM": "cdutil.region.domain(latitude=(-50.0, 50.0),longitude=(0.0, 360.0))" + }, + "id": "TRMM" + } + }, + "DISCLAIMER": "USER-NOTICE: The results in this file were produced with the PMP v1.1 (https://github.com/PCMDI/pcmdi_metrics). They are for research purposes only. They are subject to ongoing quality control and change as the PMP software advances, interpolation methods are modified, observational data sets are updated, problems with model data are corrected, etc. Use of these results for research (presentation, publications, etc.) should reference: Gleckler, P. J., C. Doutriaux, P. J. Durack, K. E. Taylor, Y. Zhang, and D. N. Williams, E. Mason, and J. Servonnat (2016), A more powerful reality test for climate models, Eos, 97, doi:10.1029/2016EO051663. If any problems are uncovered in using these results please contact the PMP development team at pcmdi-metrics@llnl.gov\n", + "json_version": 3.0 +} \ No newline at end of file diff --git a/tests/diurnal/results/json/pr_Jul_1999-2005_std_of_hourlymeans.json b/tests/diurnal/results/json/pr_Jul_1999-2005_std_of_hourlymeans.json new file mode 100644 index 000000000..eb88ac4e6 --- /dev/null +++ b/tests/diurnal/results/json/pr_Jul_1999-2005_std_of_hourlymeans.json @@ -0,0 +1,87 @@ +{ + "REFERENCE": "The statistics in this file are based on Trenberth, Zhang & Gehne, J Hydromet. 2017", + "json_structure": [ + "model", + "domain" + ], + "provenance": { + "platform": { + "OS": "Darwin", + "Version": "16.7.0", + "Name": "loki" + }, + "userId": "doutriaux1", + "osAccess": false, + "commandLine": "/Users/doutriaux1/anaconda2/envs/dev-nox/bin/std_of_hourlyvaluesWrappedInOut.py --append -i tests/diurnal/results/nc -o test_data/results/jsons -m7", + "date": "2017-11-10 12:27:00", + "conda": { + "Platform": "osx-64", + "Version": "4.3.29", + "IsPrivate": "False", + "envVersion": "4.3.29", + "buildVersion": "3.0.9", + "PythonVersion": "2.7.13.final.0", + "RootEnvironment": "/Users/doutriaux1/anaconda2 (writable)", + "DefaultEnvironment": "/Users/doutriaux1/anaconda2/envs/dev-nox" + }, + "packages": { + "cdms": "2.12", + "CDP": "1.1.0", + "cdtime": "2.12", + "cdutil": "2.12", + "clapack": "3.2.1", + "lapack": "3.6.1", + "esmf": "7.0.0", + "esmpy": "7.0.0", + "genutil": "2.12", + "python": "2.7.14", + "mesalib": "17.2.0", + "numpy": "1.13.1", + "vcs": "2.12", + "vtk": "7.1.0.2.12", + "blas": null, + "matplotlib": null, + "PMP": "g9a0d4e8", + "PMPObs": null + }, + "openGL": { + "vendor": "ATI Technologies Inc.", + "renderer": "AMD Radeon Pro 455 OpenGL Engine", + "version": "2.1 ATI-1.51.8", + "shading language version": "1.20", + "GLX": { + "client": { + "vendor": "Mesa Project and SGI", + "version": "1.4" + }, + "version": "1.4", + "server": { + "vendor": "SGI", + "version": "1.4" + } + } + } + }, + "RESULTS": { + "CMCC": { + "TROPICS": 10.966332640748378, + "TRMM": 10.77210397230895 + } + }, + "RegionalMasking": { + "TROPICS": { + "domain": { + "TROPICS": "cdutil.region.domain(latitude=(-30.0, 30.0),longitude=(0.0, 360.0))" + }, + "id": "TROPICS" + }, + "TRMM": { + "domain": { + "TRMM": "cdutil.region.domain(latitude=(-50.0, 50.0),longitude=(0.0, 360.0))" + }, + "id": "TRMM" + } + }, + "DISCLAIMER": "USER-NOTICE: The results in this file were produced with the PMP v1.1 (https://github.com/PCMDI/pcmdi_metrics). They are for research purposes only. They are subject to ongoing quality control and change as the PMP software advances, interpolation methods are modified, observational data sets are updated, problems with model data are corrected, etc. Use of these results for research (presentation, publications, etc.) should reference: Gleckler, P. J., C. Doutriaux, P. J. Durack, K. E. Taylor, Y. Zhang, and D. N. Williams, E. Mason, and J. Servonnat (2016), A more powerful reality test for climate models, Eos, 97, doi:10.1029/2016EO051663. If any problems are uncovered in using these results please contact the PMP development team at pcmdi-metrics@llnl.gov\n", + "json_version": 3.0 +} \ No newline at end of file diff --git a/tests/diurnal/results/json/pr_Jul_1999-2005_std_of_meandiurnalcyc.json b/tests/diurnal/results/json/pr_Jul_1999-2005_std_of_meandiurnalcyc.json new file mode 100644 index 000000000..906f17e00 --- /dev/null +++ b/tests/diurnal/results/json/pr_Jul_1999-2005_std_of_meandiurnalcyc.json @@ -0,0 +1,87 @@ +{ + "REFERENCE": "The statistics in this file are based on Trenberth, Zhang & Gehne, J Hydromet. 2017", + "json_structure": [ + "model", + "domain" + ], + "provenance": { + "platform": { + "OS": "Darwin", + "Version": "16.7.0", + "Name": "loki" + }, + "userId": "doutriaux1", + "osAccess": false, + "commandLine": "/Users/doutriaux1/anaconda2/envs/dev-nox/bin/std_of_meandiurnalcycWrappedInOut.py --append -i tests/diurnal/results/nc -o test_data/results/jsons -m7", + "date": "2017-11-10 12:27:09", + "conda": { + "Platform": "osx-64", + "Version": "4.3.29", + "IsPrivate": "False", + "envVersion": "4.3.29", + "buildVersion": "3.0.9", + "PythonVersion": "2.7.13.final.0", + "RootEnvironment": "/Users/doutriaux1/anaconda2 (writable)", + "DefaultEnvironment": "/Users/doutriaux1/anaconda2/envs/dev-nox" + }, + "packages": { + "cdms": "2.12", + "CDP": "1.1.0", + "cdtime": "2.12", + "cdutil": "2.12", + "clapack": "3.2.1", + "lapack": "3.6.1", + "esmf": "7.0.0", + "esmpy": "7.0.0", + "genutil": "2.12", + "python": "2.7.14", + "mesalib": "17.2.0", + "numpy": "1.13.1", + "vcs": "2.12", + "vtk": "7.1.0.2.12", + "blas": null, + "matplotlib": null, + "PMP": "g9a0d4e8", + "PMPObs": null + }, + "openGL": { + "vendor": "ATI Technologies Inc.", + "renderer": "AMD Radeon Pro 455 OpenGL Engine", + "version": "2.1 ATI-1.51.8", + "shading language version": "1.20", + "GLX": { + "client": { + "vendor": "Mesa Project and SGI", + "version": "1.4" + }, + "version": "1.4", + "server": { + "vendor": "SGI", + "version": "1.4" + } + } + } + }, + "RESULTS": { + "CMCC": { + "TROPICS": 1.2867014136692358, + "TRMM": 1.1373195943304841 + } + }, + "RegionalMasking": { + "TROPICS": { + "domain": { + "TROPICS": "cdutil.region.domain(latitude=(-30.0, 30.0),longitude=(0.0, 360.0))" + }, + "id": "TROPICS" + }, + "TRMM": { + "domain": { + "TRMM": "cdutil.region.domain(latitude=(-50.0, 50.0),longitude=(0.0, 360.0))" + }, + "id": "TRMM" + } + }, + "DISCLAIMER": "USER-NOTICE: The results in this file were produced with the PMP v1.1 (https://github.com/PCMDI/pcmdi_metrics). They are for research purposes only. They are subject to ongoing quality control and change as the PMP software advances, interpolation methods are modified, observational data sets are updated, problems with model data are corrected, etc. Use of these results for research (presentation, publications, etc.) should reference: Gleckler, P. J., C. Doutriaux, P. J. Durack, K. E. Taylor, Y. Zhang, and D. N. Williams, E. Mason, and J. Servonnat (2016), A more powerful reality test for climate models, Eos, 97, doi:10.1029/2016EO051663. If any problems are uncovered in using these results please contact the PMP development team at pcmdi-metrics@llnl.gov\n", + "json_version": 3.0 +} \ No newline at end of file diff --git a/tests/diurnal/results/json/pr_Jul_1999_2005_std_of_dailymeans.json b/tests/diurnal/results/json/pr_Jul_1999_2005_std_of_dailymeans.json new file mode 100644 index 000000000..f7b6c735a --- /dev/null +++ b/tests/diurnal/results/json/pr_Jul_1999_2005_std_of_dailymeans.json @@ -0,0 +1,87 @@ +{ + "REFERENCE": "The statistics in this file are based on Trenberth, Zhang & Gehne, J Hydromet. 2017", + "json_structure": [ + "model", + "domain" + ], + "provenance": { + "platform": { + "OS": "Darwin", + "Version": "16.7.0", + "Name": "loki" + }, + "userId": "doutriaux1", + "osAccess": false, + "commandLine": "/Users/doutriaux1/anaconda2/envs/dev-nox/bin/std_of_dailymeansWrappedInOut.py --append -i tests/diurnal/results/nc -o test_data/results/jsons -m7", + "date": "2017-11-10 12:26:38", + "conda": { + "Platform": "osx-64", + "Version": "4.3.29", + "IsPrivate": "False", + "envVersion": "4.3.29", + "buildVersion": "3.0.9", + "PythonVersion": "2.7.13.final.0", + "RootEnvironment": "/Users/doutriaux1/anaconda2 (writable)", + "DefaultEnvironment": "/Users/doutriaux1/anaconda2/envs/dev-nox" + }, + "packages": { + "cdms": "2.12", + "CDP": "1.1.0", + "cdtime": "2.12", + "cdutil": "2.12", + "clapack": "3.2.1", + "lapack": "3.6.1", + "esmf": "7.0.0", + "esmpy": "7.0.0", + "genutil": "2.12", + "python": "2.7.14", + "mesalib": "17.2.0", + "numpy": "1.13.1", + "vcs": "2.12", + "vtk": "7.1.0.2.12", + "blas": null, + "matplotlib": null, + "PMP": "g9a0d4e8", + "PMPObs": null + }, + "openGL": { + "vendor": "ATI Technologies Inc.", + "renderer": "AMD Radeon Pro 455 OpenGL Engine", + "version": "2.1 ATI-1.51.8", + "shading language version": "1.20", + "GLX": { + "client": { + "vendor": "Mesa Project and SGI", + "version": "1.4" + }, + "version": "1.4", + "server": { + "vendor": "SGI", + "version": "1.4" + } + } + } + }, + "RESULTS": { + "CMCC": { + "TROPICS": 7.860331218213006, + "TRMM": 7.396705786166231 + } + }, + "RegionalMasking": { + "TROPICS": { + "domain": { + "TROPICS": "cdutil.region.domain(latitude=(-30.0, 30.0),longitude=(0.0, 360.0))" + }, + "id": "TROPICS" + }, + "TRMM": { + "domain": { + "TRMM": "cdutil.region.domain(latitude=(-50.0, 50.0),longitude=(0.0, 360.0))" + }, + "id": "TRMM" + } + }, + "DISCLAIMER": "USER-NOTICE: The results in this file were produced with the PMP v1.1 (https://github.com/PCMDI/pcmdi_metrics). They are for research purposes only. They are subject to ongoing quality control and change as the PMP software advances, interpolation methods are modified, observational data sets are updated, problems with model data are corrected, etc. Use of these results for research (presentation, publications, etc.) should reference: Gleckler, P. J., C. Doutriaux, P. J. Durack, K. E. Taylor, Y. Zhang, and D. N. Williams, E. Mason, and J. Servonnat (2016), A more powerful reality test for climate models, Eos, 97, doi:10.1029/2016EO051663. If any problems are uncovered in using these results please contact the PMP development team at pcmdi-metrics@llnl.gov\n", + "json_version": 3.0 +} \ No newline at end of file diff --git a/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_S.nc b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_S.nc new file mode 100644 index 000000000..a91ffa3b0 Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_S.nc differ diff --git a/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_diurnal_avg.nc b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_diurnal_avg.nc new file mode 100644 index 000000000..65566c3d9 Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_diurnal_avg.nc differ diff --git a/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_diurnal_std.nc b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_diurnal_std.nc new file mode 100644 index 000000000..01acf6e82 Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_diurnal_std.nc differ diff --git a/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_std_of_dailymeans.nc b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_std_of_dailymeans.nc new file mode 100644 index 000000000..361cafc41 Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_std_of_dailymeans.nc differ diff --git a/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_tS.nc b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_tS.nc new file mode 100644 index 000000000..97f4015c3 Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_tS.nc differ diff --git a/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_tmean.nc b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_tmean.nc new file mode 100644 index 000000000..4358e3ab8 Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_Jul_1999-2005_tmean.nc differ diff --git a/tests/diurnal/results/nc/pr_CMCC_LocalSolarTimes.nc b/tests/diurnal/results/nc/pr_CMCC_LocalSolarTimes.nc new file mode 100644 index 000000000..3d5f379ec Binary files /dev/null and b/tests/diurnal/results/nc/pr_CMCC_LocalSolarTimes.nc differ diff --git a/tests/test_pmp_diurnal.py b/tests/test_pmp_diurnal.py new file mode 100644 index 000000000..c97e62425 --- /dev/null +++ b/tests/test_pmp_diurnal.py @@ -0,0 +1,82 @@ +import unittest +import cdat_info +import subprocess +import shlex +import os +import cdms2 +import numpy +import json + + +class DiurnalTest(unittest.TestCase): + + def assertSame(self,a,b): + self.assertTrue(numpy.ma.allclose(a,b)) + + def compare_nc(self,test_name): + print "Checking integrity of",test_name + self.assertTrue(os.path.exists(os.path.join("test_data",test_name))) + good_name = test_name.replace("test_data","tests/diurnal") + + test_out = cdms2.open(os.path.join("test_data",test_name)) + good_out = cdms2.open(os.path.join("tests/diurnal",test_name)) + + self.assertEqual(test_out.variables.keys(),good_out.variables.keys()) + + for v in good_out.variables.keys(): + test = test_out(v) + good = good_out(v) + + self.assertSame(test,good) + + def testDiurnaliComputeStdDailyMean(self): + cmd = 'computeStdDailyMeansWrapped.py -i test_data -o test_data/results/nc -t "sample_data_pr_%(model).nc" -m7' + p = subprocess.Popen(shlex.split(cmd)) + p.communicate() + + self.compare_nc("results/nc/pr_CMCC_Jul_1999-2005_std_of_dailymeans.nc") + + def testFourierDiurnalAllGridWrapped(self): + cmd = 'fourierDiurnalAllGridWrapped.py -i test_data/results/nc -o test_data/results/nc -m7' + p = subprocess.Popen(shlex.split(cmd)) + p.communicate() + self.compare_nc("results/nc/pr_CMCC_Jul_1999-2005_tmean.nc") + self.compare_nc("results/nc/pr_CMCC_Jul_1999-2005_tS.nc") + self.compare_nc("results/nc/pr_CMCC_Jul_1999-2005_S.nc") + + def testDiurnalStdDailyVariance(self): + self.runJsoner("std_of_dailymeansWrappedInOut.py","pr_Jul_1999_2005_std_of_dailymeans.json") + def runJsoner(self,script,json_file): + cmd = '{} --region_name=TROPICS --lat1=-30. --lat2=30. --lon1=0. --lon2=360 -i tests/diurnal/results/nc -o test_data/results/jsons -m7'.format(script) + p = subprocess.Popen(shlex.split(cmd)) + p.communicate() + cmd = '{} --append -i tests/diurnal/results/nc -o test_data/results/jsons -m7'.format(script) + p = subprocess.Popen(shlex.split(cmd)) + p.communicate() + good = open("tests/diurnal/results/json/{}".format(json_file)) + test = open("test_data/results/jsons/{}".format(json_file)) + test = json.load(test) + good = json.load(good) + self.assertEqual(test["RESULTS"],good["RESULTS"]) + def testCompositeDiurnalStatisticsWrapped(self): + cmd = 'compositeDiurnalStatisticsWrapped.py -i test_data -o test_data/results/nc -t "sample_data_pr_%(model).nc" -m7' + p = subprocess.Popen(shlex.split(cmd)) + p.communicate() + self.compare_nc("results/nc/pr_CMCC_Jul_1999-2005_diurnal_avg.nc") + self.compare_nc("results/nc/pr_CMCC_Jul_1999-2005_diurnal_std.nc") + self.compare_nc("results/nc/pr_CMCC_LocalSolarTimes.nc") + + def testStd_of_hourlyvaluesWrappedInOut(self): + self.runJsoner("std_of_hourlyvaluesWrappedInOut.py","pr_Jul_1999-2005_std_of_hourlymeans.json") + + def testStd_of_meandiurnalcycWrappedInOut(self): + self.runJsoner("std_of_meandiurnalcycWrappedInOut.py","pr_Jul_1999-2005_std_of_meandiurnalcyc.json") + + def testSavg_fourierWrappedInOut(self): + self.runJsoner("savg_fourierWrappedInOut.py","pr_Jul_1999-2005_savg_DiurnalFourier.json") + + def testfourierDiurnalGridpoints(self): + cmd = "fourierDiurnalGridpoints.py -i tests/diurnal/results/nc -o test_data/results/ascii" + p = subprocess.Popen(shlex.split(cmd)) + p.communicate() + self.assertTrue(os.path.exists("test_data/results/ascii/pr_Jul_1999-2005_fourierDiurnalGridPoints.asc")) diff --git a/tests/test_pmp_flake8.py b/tests/test_pmp_flake8.py index 83bed0667..9a73da19b 100644 --- a/tests/test_pmp_flake8.py +++ b/tests/test_pmp_flake8.py @@ -26,8 +26,7 @@ def testFlake8(self): stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - P.wait() - out = P.stdout.read() + out,e = P.communicate() if out != "": print out self.assertEqual(out, "")