diff --git a/TODO.md b/TODO.md index 1956d2e..60f6376 100755 --- a/TODO.md +++ b/TODO.md @@ -5,6 +5,6 @@ - [X] Add Fernando's notebooks - [X] Add Johan's notebooks - [X] Convert all code from Py2 to Py3 +- [X] Select more utilities to add for version 0.1.0 - [ ] Check install works flawlessly -- [ ] Select more utilities to add for version 1 - [ ] Add header info to all codes diff --git a/captoolkit/corrlaser.py b/captoolkit/corrlaser.py index 3c89376..175a643 100755 --- a/captoolkit/corrlaser.py +++ b/captoolkit/corrlaser.py @@ -106,14 +106,14 @@ def main(fname): files = sys.argv[1:] -print("Number of files:", len(files)) +print(("Number of files:", len(files))) if njobs == 1: print("running sequential code ...") [main(f) for f in files] else: - print("running parallel code (%d jobs) ..." % njobs) + print(("running parallel code (%d jobs) ..." % njobs)) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)(delayed(main)(f) for f in files) diff --git a/captoolkit/corrlaser.py.bak b/captoolkit/corrlaser.py.bak index 9f4c77a..3c89376 100755 --- a/captoolkit/corrlaser.py.bak +++ b/captoolkit/corrlaser.py.bak @@ -1,93 +1,121 @@ #!/usr/bin/env python """ -Apply corrections for ICESat Laser 2 and 3. +Compute and apply corrections for ICESat Laser 2 and 3. From Borsa et al. (2019): Subtract 1.7 cm from Laser 2 and add 1.1 cm to Laser 3 -""" -from __future__ import print_function -from __future__ import unicode_literals -from __future__ import division -from __future__ import absolute_import -from future import standard_library -standard_library.install_aliases() +Notes: + Edit some parameters bellow. + +Example: + python corrlaser.py /path/to/files/*.h5 + +Credits: + captoolkit - JPL Cryosphere Altimetry Processing Toolkit + + Fernando Paolo (paolofer@jpl.nasa.gov) + Johan Nilsson (johan.nilsson@jpl.nasa.gov) + Alex Gardner (alex.s.gardner@jpl.nasa.gov) + + Jet Propulsion Laboratory, California Institute of Technology + +""" import os import sys + import h5py import numpy as np -from glob import glob -from joblib import Parallel, delayed from astropy.time import Time +from future import standard_library + +standard_library.install_aliases() + +# === EDIT HERE ====================================== # + +# time variable +tvar = "t_year" -tvar = 't_year' -hvar = 'h_cor' +# height variable +hvar = "h_cor" +# apply or only store the correction for each height apply_ = True + +# number of jobs in parallel njobs = 16 +# === END EDIT ======================================= # # Corrections for lasers/campaigns (in meters) -# Cor will be subtracted from height: h - lx -bias = {'l1': 0., 'l2': 0.017, 'l3': -0.011, None: 0.} +# Cor will be subtracted from height: h - l{1,2,3} +bias = {"l1": 0.0, "l2": 0.017, "l3": -0.011, None: 0.0} # Definition of ICESat campaigns/lasers # https://nsidc.org/data/icesat/laser_op_periods.html campaigns = [ - (Time('2003-02-20').decimalyear, Time('2003-03-21').decimalyear, 'l1a'), - (Time('2003-03-21').decimalyear, Time('2003-03-29').decimalyear, 'l1b'), - (Time('2003-09-25').decimalyear, Time('2003-11-19').decimalyear, 'l2a'), # 8-d + 91-d + 91-d (NSIDC) - (Time('2004-02-17').decimalyear, Time('2004-03-21').decimalyear, 'l2b'), - (Time('2004-05-18').decimalyear, Time('2004-06-21').decimalyear, 'l2c'), - (Time('2004-10-03').decimalyear, Time('2004-11-08').decimalyear, 'l3a'), - (Time('2005-02-17').decimalyear, Time('2005-03-24').decimalyear, 'l3b'), - (Time('2005-05-20').decimalyear, Time('2005-06-23').decimalyear, 'l3c'), - (Time('2005-10-21').decimalyear, Time('2005-11-24').decimalyear, 'l3d'), - (Time('2006-02-22').decimalyear, Time('2006-03-28').decimalyear, 'l3e'), - (Time('2006-05-24').decimalyear, Time('2006-06-26').decimalyear, 'l3f'), - (Time('2006-10-25').decimalyear, Time('2006-11-27').decimalyear, 'l3g'), - (Time('2007-03-12').decimalyear, Time('2007-04-14').decimalyear, 'l3h'), - (Time('2007-10-02').decimalyear, Time('2007-11-05').decimalyear, 'l3i'), - (Time('2008-02-17').decimalyear, Time('2008-03-21').decimalyear, 'l3j'), - (Time('2008-10-04').decimalyear, Time('2008-10-19').decimalyear, 'l3k'), - (Time('2008-11-25').decimalyear, Time('2008-12-17').decimalyear, 'l2d'), - (Time('2009-03-09').decimalyear, Time('2009-04-11').decimalyear, 'l2e'), - (Time('2009-09-30').decimalyear, Time('2009-10-11').decimalyear, 'l2f'), - ] + (Time("2003-02-20").decimalyear, Time("2003-03-21").decimalyear, "l1a"), + (Time("2003-03-21").decimalyear, Time("2003-03-29").decimalyear, "l1b"), + ( + Time("2003-09-25").decimalyear, + Time("2003-11-19").decimalyear, + "l2a", + ), # 8-d + 91-d + 91-d (NSIDC) + (Time("2004-02-17").decimalyear, Time("2004-03-21").decimalyear, "l2b"), + (Time("2004-05-18").decimalyear, Time("2004-06-21").decimalyear, "l2c"), + (Time("2004-10-03").decimalyear, Time("2004-11-08").decimalyear, "l3a"), + (Time("2005-02-17").decimalyear, Time("2005-03-24").decimalyear, "l3b"), + (Time("2005-05-20").decimalyear, Time("2005-06-23").decimalyear, "l3c"), + (Time("2005-10-21").decimalyear, Time("2005-11-24").decimalyear, "l3d"), + (Time("2006-02-22").decimalyear, Time("2006-03-28").decimalyear, "l3e"), + (Time("2006-05-24").decimalyear, Time("2006-06-26").decimalyear, "l3f"), + (Time("2006-10-25").decimalyear, Time("2006-11-27").decimalyear, "l3g"), + (Time("2007-03-12").decimalyear, Time("2007-04-14").decimalyear, "l3h"), + (Time("2007-10-02").decimalyear, Time("2007-11-05").decimalyear, "l3i"), + (Time("2008-02-17").decimalyear, Time("2008-03-21").decimalyear, "l3j"), + (Time("2008-10-04").decimalyear, Time("2008-10-19").decimalyear, "l3k"), + (Time("2008-11-25").decimalyear, Time("2008-12-17").decimalyear, "l2d"), + (Time("2009-03-09").decimalyear, Time("2009-04-11").decimalyear, "l2e"), + (Time("2009-09-30").decimalyear, Time("2009-10-11").decimalyear, "l2f"), +] def _get_laser_bias(time, campaigns, bias): """ Map time (yr) to campaign to correction. """ - camp = [ca for (t1,t2,ca) in campaigns if t1 <= time < t2] # ['c']|[] + camp = [ca for (t1, t2, ca) in campaigns if t1 <= time < t2] # ['c']|[] laser = camp[0][:2] if camp else None + return bias[laser] -get_laser_bias = np.vectorize(_get_laser_bias, excluded=[1,2], cache=True) + +get_laser_bias = np.vectorize(_get_laser_bias, excluded=[1, 2], cache=True) def main(fname): - with h5py.File(fname, 'a') as f: + with h5py.File(fname, "a") as f: t = f[tvar][:] - h = f[hvar][:] b = get_laser_bias(t, campaigns, bias) - f['laser_bias'] = b # save bias - if apply_: f[hvar][:] -= b # remove bias + f["laser_bias"] = b # save bias + + if apply_: + f[hvar][:] -= b # remove bias f.flush() - os.rename(fname, fname.replace('.h5', '_LCOR.h5')) + os.rename(fname, fname.replace(".h5", "_LCOR.h5")) files = sys.argv[1:] -print('Number of files:', len(files)) +print("Number of files:", len(files)) if njobs == 1: - print('running sequential code ...') + print("running sequential code ...") [main(f) for f in files] else: - print('running parallel code (%d jobs) ...' % njobs) + print("running parallel code (%d jobs) ..." % njobs) from joblib import Parallel, delayed + Parallel(n_jobs=njobs, verbose=5)(delayed(main)(f) for f in files) -print('done.') +print("done.") diff --git a/captoolkit/corrscatt.py b/captoolkit/corrscatt.py index bb163a3..58b56f6 100755 --- a/captoolkit/corrscatt.py +++ b/captoolkit/corrscatt.py @@ -694,7 +694,7 @@ def multi_fit_coef(t_, h_, bs_, lew_, tes_): except IOError: print("MULTIVARIATE FIT FAILED, setting params -> 0") - print("VALID DATA POINTS in h:", sum(~np.isnan(h_))) + print(("VALID DATA POINTS in h:", sum(~np.isnan(h_)))) a, b, c, r2, pval, pvals = 0, 0, 0, 0, 1e3, [1e3, 1e3, 1e3] @@ -1008,27 +1008,27 @@ def func(t, c, m, n, a, p): print("Summary:") print("--------") - print("cor applied: ", (h_bs[~np.isnan(h_bs)] != 0).any()) - print( + print(("cor applied: ", (h_bs[~np.isnan(h_bs)] != 0).any())) + print(( "std change: %.3f m (%.1f %%)" % (round(d_std, 3), round(p_std * 100, 1)) - ) - print( + )) + print(( "trend change: %.3f m/yr (%.1f %%)" % (round(d_trend, 3), round(p_trend * 100, 1)) - ) + )) print("") - print("r-squared: ", round(r2, 3)) - print("p-value: ", round(pval, 3)) - print("p-values: ", [round(p, 3) for p in pvals]) + print(("r-squared: ", round(r2, 3))) + print(("p-value: ", round(pval, 3))) + print(("p-values: ", [round(p, 3) for p in pvals])) print("") - print("r_bs: ", round(r_bc, 3)) - print("r_lew: ", round(r_wc, 3)) - print("r_tes: ", round(r_sc, 3)) + print(("r_bs: ", round(r_bc, 3))) + print(("r_lew: ", round(r_wc, 3))) + print(("r_tes: ", round(r_sc, 3))) print("") - print("s_bs: ", round(s_bc, 3)) - print("s_lew: ", round(s_wc, 3)) - print("s_tes: ", round(s_sc, 3)) + print(("s_bs: ", round(s_bc, 3))) + print(("s_lew: ", round(s_wc, 3))) + print(("s_tes: ", round(s_sc, 3))) plt.show() @@ -1054,7 +1054,7 @@ def main( print("* RUNNING IN TEST MODE (PLOTTING ONLY, NOT SAVING DATA) *") print("*********************************************************") - print("processing file:", ifile, "...") + print(("processing file:", ifile, "...")) # Test if parameter file exists if "_scatgrd" in ifile.lower(): @@ -1161,7 +1161,7 @@ def main( for k in range(N_nodes): if (k % 500) == 0: - print("Calculating correction for node", k, "of", N_nodes, "...") + print(("Calculating correction for node", k, "of", N_nodes, "...")) x0, y0 = x_nodes[k], y_nodes[k] @@ -1584,7 +1584,7 @@ def main( bbox = args.bbox[:] print("parameters:") - for arg in vars(args).items(): + for arg in list(vars(args).items()): print(arg) if njobs == 1: @@ -1596,7 +1596,7 @@ def main( for ifile in ifiles ] else: - print("running parallel code (%d jobs) ..." % njobs) + print(("running parallel code (%d jobs) ..." % njobs)) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( diff --git a/captoolkit/corrscatt.py.bak b/captoolkit/corrscatt.py.bak index ff3619c..bb163a3 100755 --- a/captoolkit/corrscatt.py.bak +++ b/captoolkit/corrscatt.py.bak @@ -1,23 +1,37 @@ #!/usr/bin/env python -# -*- coding: utf-8 -*- """ Corrects radar altimetry height to correlation with waveform parameters. -Example: - scattcor.py -v lon lat h_res t_year -w bs lew tes -d 1 -r 4 -q 2 - -p dif -f /path/to/*files.h5 +Estimates the apparent height change originating from radar surface +scattering variations (which changes the waveform shape). Notes: - The (back)scattering correction is applied as: + - This is a complex code; and it has been tailored to a specific use case. + + - The (back)scattering correction is applied as: + h_corr = h - h_bs + + - Edit parameters below. + +Example: + python scattcor.py -v lon lat h_res t_year -w bs lew tes -d 1 -r 4 \ + q 2 -p dif -f /path/to/*files.h5 + + python scattcor.py -h + +Credits: + captoolkit - JPL Cryosphere Altimetry Processing Toolkit - hc_cor = h - h_bs + Fernando Paolo (paolofer@jpl.nasa.gov) + Johan Nilsson (johan.nilsson@jpl.nasa.gov) + Alex Gardner (alex.s.gardner@jpl.nasa.gov) + + Jet Propulsion Laboratory, California Institute of Technology """ import os -import sys import h5py import pyproj -import timeit import warnings import argparse import numpy as np @@ -27,12 +41,12 @@ from scipy.stats import mode from scipy.spatial import cKDTree from scipy.signal import savgol_filter -from sklearn.metrics import mean_squared_error from math import sqrt +from sklearn.metrics import mean_squared_error -#--- Edit ------------------------------------------------------------ +# --- Edit ------------------------------------------------------------ -##NOTE: This uses random cells, plot results, and do not save data +# OTE: This uses random cells, plot results, and do not save data TEST_MODE = False # Use random locations @@ -44,22 +58,24 @@ SEED = 222 USE_NODES = False # Specific locations for testing: Ross, Getz, PIG -NODES = [(-158.71, -78.7584), - #(-124.427, -74.4377), # Getz - #(-100.97, -75.1478), # PIG - (-158.40, -78.80), - (-178.40, -78.80), - (-188.00, -77.95), - (-160.00, -80.40), - (-178.40, -80.60), - (-190.40, -80.60),] +NODES = [ + (-158.71, -78.7584), + # (-124.427, -74.4377), # Getz + # (-100.97, -75.1478), # PIG + (-158.40, -78.80), + (-178.40, -78.80), + (-188.00, -77.95), + (-160.00, -80.40), + (-178.40, -80.60), + (-190.40, -80.60), +] # Suffix for output file -SUFFIX1 = '_SCAT' -SUFFIX2 = '_SCATGRD' +SUFFIX1 = "_SCAT" +SUFFIX2 = "_SCATGRD" # Nave of variable to save bs-correction -H_BS = 'h_bs' +H_BS = "h_bs" # Apply 3-month running median to processed time series BIN_SERIES = True @@ -70,7 +86,7 @@ R_MIN = 0.1 # Minimum points per cell to compute solution MIN_PTS = 50 -# Minimum number of months to compute solution +# Minimum number of months to compute solution MIN_MONTHS = 24 # Default time range @@ -81,74 +97,142 @@ WINDOW = 15 ORDER = 1 DERIV = 1 -#---------------------------------------------------------------- +# ---------------------------------------------------------------- # Supress anoying warnings -warnings.filterwarnings('ignore') +warnings.filterwarnings("ignore") def get_args(): """ Get command-line arguments. """ - msg = 'Correct height data for surface scattering variation.' + msg = "Correct height data for surface scattering variation." parser = argparse.ArgumentParser(description=msg) parser.add_argument( - '-f', metavar='file', dest='files', type=str, nargs='+', - help='HDF5 file(s) to process', - required=True) + "-f", + metavar="file", + dest="files", + type=str, + nargs="+", + help="HDF5 file(s) to process", + required=True, + ) parser.add_argument( - '-d', metavar=('length'), dest='dxy', type=float, nargs=1, - help=('block size of grid cell (km)'), - default=[None], required=True) + "-d", + metavar=("length"), + dest="dxy", + type=float, + nargs=1, + help=("block size of grid cell (km)"), + default=[None], + required=True, + ) parser.add_argument( - '-r', metavar=('radius'), dest='radius', type=float, nargs=1, - help=('search radius (w/relocation) (km)'), - default=[0],) + "-r", + metavar=("radius"), + dest="radius", + type=float, + nargs=1, + help=("search radius (w/relocation) (km)"), + default=[0], + ) parser.add_argument( - '-q', metavar=('n_reloc'), dest='nreloc', type=int, nargs=1, - help=('number of relocations for search radius'), - default=[0],) + "-q", + metavar=("n_reloc"), + dest="nreloc", + type=int, + nargs=1, + help=("number of relocations for search radius"), + default=[0], + ) parser.add_argument( - '-p', metavar=None, dest='proc', type=str, nargs=1, - help=('pre-process series for sensitivity estimation'), - choices=('det', 'dif'), default=[None],) + "-p", + metavar=None, + dest="proc", + type=str, + nargs=1, + help=("pre-process series for sensitivity estimation"), + choices=("det", "dif"), + default=[None], + ) parser.add_argument( - '-v', metavar=('lon','lat', 'h', 't'), dest='vnames', - type=str, nargs=4, - help=('name of x/y/z/t variables in the HDF5'), - default=[None], required=True) + "-v", + metavar=("lon", "lat", "h", "t"), + dest="vnames", + type=str, + nargs=4, + help=("name of x/y/z/t variables in the HDF5"), + default=[None], + required=True, + ) parser.add_argument( - '-w', metavar=('bs', 'lew', 'tes'), dest='wnames', - type=str, nargs=3, - help=('name of sig0/LeW/TeS parameters in HDF5'), - default=[None], required=True) + "-w", + metavar=("bs", "lew", "tes"), + dest="wnames", + type=str, + nargs=3, + help=("name of sig0/LeW/TeS parameters in HDF5"), + default=[None], + required=True, + ) parser.add_argument( - '-j', metavar=('epsg_num'), dest='proj', type=str, nargs=1, - help=('EPSG proj number (AnIS=3031, GrIS=3413)'), - default=['3031'],) + "-j", + metavar=("epsg_num"), + dest="proj", + type=str, + nargs=1, + help=("EPSG proj number (AnIS=3031, GrIS=3413)"), + default=["3031"], + ) parser.add_argument( - '-n', metavar=('n_jobs'), dest='njobs', type=int, nargs=1, - help="number of jobs for parallel processing", - default=[1],) + "-n", + metavar=("n_jobs"), + dest="njobs", + type=int, + nargs=1, + help="number of jobs for parallel processing", + default=[1], + ) parser.add_argument( - '-t', metavar=('tmin','tmax'), dest='tlim', type=float, nargs=2, - help="time interval to compute corrections (dec years)", - default=[TMIN,TMAX],) + "-t", + metavar=("tmin", "tmax"), + dest="tlim", + type=float, + nargs=2, + help="time interval to compute corrections (dec years)", + default=[TMIN, TMAX], + ) parser.add_argument( - '-b', metavar=('e','w','s','n'), dest='bbox', type=float, nargs=4, - help="full bbox in case of processing tiles for consistency", - default=[None],) + "-b", + metavar=("e", "w", "s", "n"), + dest="bbox", + type=float, + nargs=4, + help="full bbox in case of processing tiles for consistency", + default=[None], + ) parser.add_argument( - '-a', dest='apply', action='store_true', - help=('apply correction to height in addition to saving'), - default=False) + "-a", + dest="apply", + action="store_true", + help=("apply correction to height in addition to saving"), + default=False, + ) return parser.parse_args() """ Generic functions """ -def binning(x, y, xmin=None, xmax=None, dx=1/12., window=3/12., - interp=False, median=False): +def binning( + x, + y, + xmin=None, + xmax=None, + dx=1 / 12.0, + window=3 / 12.0, + interp=False, + median=False, +): """ Time-series binning (w/overlapping windows). @@ -159,11 +243,13 @@ def binning(x, y, xmin=None, xmax=None, dx=1/12., window=3/12., window: size of binning window. interp: interpolate binned values to original x points. """ - if xmin is None: xmin = np.nanmin(x) - if xmax is None: xmax = np.nanmax(x) + if xmin is None: + xmin = np.nanmin(x) + if xmax is None: + xmax = np.nanmax(x) - steps = np.arange(xmin, xmax+dx, dx) - bins = [(ti, ti+window) for ti in steps] + steps = np.arange(xmin, xmax + dx, dx) + bins = [(ti, ti + window) for ti in steps] N = len(bins) @@ -173,22 +259,23 @@ def binning(x, y, xmin=None, xmax=None, dx=1/12., window=3/12., nb = np.full(N, np.nan) sb = np.full(N, np.nan) - for i in xrange(N): + for i in range(N): t1, t2 = bins[i] - idx, = np.where((x >= t1) & (x <= t2)) + (idx,) = np.where((x >= t1) & (x <= t2)) - if len(idx) == 0: continue + if len(idx) == 0: + continue ybv = y[idx] - xbv = x[idx] + # xbv = x[idx] if median: yb[i] = np.nanmedian(ybv) else: yb[i] = np.nanmean(ybv) - xb[i] = 0.5 * (t1+t2) + xb[i] = 0.5 * (t1 + t2) eb[i] = mad_std(ybv) nb[i] = np.sum(~np.isnan(ybv)) sb[i] = np.sum(ybv) @@ -205,8 +292,8 @@ def binning(x, y, xmin=None, xmax=None, dx=1/12., window=3/12., def transform_coord(proj1, proj2, x, y): """ Transform coordinates from proj1 to proj2 (EPSG num). """ # Set full EPSG projection strings - proj1 = pyproj.Proj("+init=EPSG:"+str(proj1)) - proj2 = pyproj.Proj("+init=EPSG:"+str(proj2)) + proj1 = pyproj.Proj("+init=EPSG:" + str(proj1)) + proj2 = pyproj.Proj("+init=EPSG:" + str(proj2)) # Convert coordinates return pyproj.transform(proj1, proj2, x, y) @@ -217,8 +304,7 @@ def mad_std(x, axis=None): def mode_filter(x, min_count=10, maxiter=3): - """ - Iterative mode filter. + """Iterative mode filter. Remove values repeated 'min_count' times. """ @@ -228,7 +314,7 @@ def mode_filter(x, min_count=10, maxiter=3): if len(count) < 1: break if count[0] > min_count: - x[x==mod[0]] = np.nan + x[x == mod[0]] = np.nan n_iter += 1 else: n_iter = maxiter @@ -236,53 +322,57 @@ def mode_filter(x, min_count=10, maxiter=3): def median_filter(x, n_median=3): - """ Remove values greater than n * MAD-Std. """ + """ Remove values greater than n * MAD-Std. """ x[np.abs(x) > n_median * mad_std(x)] = np.nan return x -def detrend_binned(x, y, order=1, window=3/12.): - """ Bin data (Med), compute trend (OLS) on binned, detrend original data. """ +def detrend_binned(x, y, order=1, window=3 / 12.0): + """Bin data (Med), compute trend (OLS) on binned, detrend original data.""" x_b, y_b = binning(x, y, median=True, window=window, interp=False)[:2] i_valid = ~np.isnan(y_b) x_b, y_b = x_b[i_valid], y_b[i_valid] try: coef = np.polyfit(x_b, y_b, order) - y_trend = np.polyval(coef, x) ##NOTE: Eval on full time - except: + y_trend = np.polyval(coef, x) # NOTE: Eval on full time + except IOError: y_trend = np.zeros_like(y) - return y-y_trend, y_trend + return y - y_trend, y_trend -def sigma_filter(x, y, order=1, window=3/12., n_iter=3, n_sigma=3): +def sigma_filter(x, y, order=1, window=3 / 12.0, n_iter=3, n_sigma=3): """ Detrend (binned) data, remove 3 sigma from residual, repeat. """ y_filt, y_res = y.copy(), y.copy() for _ in range(n_iter): y_res, y_trend = detrend_binned(x, y_res, order=order, window=window) - idx, = np.where(np.abs(y_res) > mad_std(y_res)*n_sigma) - if len(idx) == 0: break # if no data to filter, stop iterating + (idx,) = np.where(np.abs(y_res) > mad_std(y_res) * n_sigma) + if len(idx) == 0: + break # if no data to filter, stop iterating y_res[idx] = np.nan - if np.sum(~np.isnan(y_res)) < 10: break ##NOTE: Arbitrary min obs - y_filt[np.isnan(y_res)] = np.nan + if np.sum(~np.isnan(y_res)) < 10: + break # NOTE: Arbitrary min obs + y_filt[np.isnan(y_res)] = np.nan return y_filt -def detrend_binned2(x, y, window=3/12.): +def detrend_binned2(x, y, window=3 / 12.0): """ Bin data (Med), detrend original data with binned data. """ x_b, y_b = binning(x, y, median=True, window=window, interp=True)[:2] - return y-y_b, y_b + return y - y_b, y_b -def sigma_filter2(x, y, window=3/12., n_iter=3, n_sigma=3): +def sigma_filter2(x, y, window=3 / 12.0, n_iter=3, n_sigma=3): """ Bin data, remove binned, remove 3 sigma from residual, repeat. """ y_filt, y_res = y.copy(), y.copy() for _ in range(n_iter): y_res, y_trend = detrend_binned2(x, y_res, window=window) - idx, = np.where(np.abs(y_res) > mad_std(y_res)*n_sigma) - if len(idx) == 0: break # if no data to filter, stop iterating + (idx,) = np.where(np.abs(y_res) > mad_std(y_res) * n_sigma) + if len(idx) == 0: + break # if no data to filter, stop iterating y_res[idx] = np.nan - if sum(~np.isnan(y_res)) < 10: break ##NOTE: Arbitrary min obs - y_filt[np.isnan(y_res)] = np.nan + if sum(~np.isnan(y_res)) < 10: + break # NOTE: Arbitrary min obs + y_filt[np.isnan(y_res)] = np.nan return y_filt @@ -293,24 +383,29 @@ def center(*arrs): def normalize(*arrs): """ Normalize array(s) by std. """ - #return [a / np.nanstd(a, ddof=1) for a in arrs] + # return [a / np.nanstd(a, ddof=1) for a in arrs] return [a / mad_std(a) for a in arrs] def corr_coef(h, bs, lew, tes): - """ Get corr coef between h and w/f params. """ - idx, = np.where(~np.isnan(h) & ~np.isnan(bs) & ~np.isnan(lew) & ~np.isnan(tes)) + """ Get corr coef between h and w/f params. """ + (idx,) = np.where( + ~np.isnan(h) & ~np.isnan(bs) & ~np.isnan(lew) & ~np.isnan(tes) + ) h_, bs_, lew_, tes_ = h[idx], bs[idx], lew[idx], tes[idx] - r_bs = np.corrcoef(bs_, h_)[0,1] - r_lew = np.corrcoef(lew_, h_)[0,1] - r_tes = np.corrcoef(tes_, h_)[0,1] + r_bs = np.corrcoef(bs_, h_)[0, 1] + r_lew = np.corrcoef(lew_, h_)[0, 1] + r_tes = np.corrcoef(tes_, h_)[0, 1] return r_bs, r_lew, r_tes def corr_grad(h, bs, lew, tes, normalize=False, robust=False): - """ Get corr gradient (slope) between h and w/f params. """ - idx, = np.where(~np.isnan(h) & ~np.isnan(bs) & ~np.isnan(lew) & ~np.isnan(tes)) - if len(idx) < 3: return np.nan, np.nan, np.nan + """ Get corr gradient (slope) between h and w/f params. """ + (idx,) = np.where( + ~np.isnan(h) & ~np.isnan(bs) & ~np.isnan(lew) & ~np.isnan(tes) + ) + if len(idx) < 3: + return np.nan, np.nan, np.nan h_, bs_, lew_, tes_ = h[idx], bs[idx], lew[idx], tes[idx] try: @@ -329,11 +424,11 @@ def corr_grad(h, bs, lew, tes, normalize=False, robust=False): s_bs /= mad_std(bs_) s_lew /= mad_std(lew_) s_tes /= mad_std(tes_) - except: + except IOError: return np.nan, np.nan, np.nan return s_bs, s_lew, s_tes - + def linefit(x, y, return_coef=False): """ @@ -344,15 +439,16 @@ def linefit(x, y, return_coef=False): assert sum(~np.isnan(y)) > 1 X = sm.add_constant(x, prepend=False) - y_fit = sm.RLM(y, X, M=sm.robust.norms.HuberT(), missing="drop") \ - .fit(maxiter=1, tol=0.001) - + y_fit = sm.RLM(y, X, M=sm.robust.norms.HuberT(), missing="drop").fit( + maxiter=1, tol=0.001 + ) + x_fit = x[~np.isnan(y)] if return_coef: - if len(y_fit.params) < 2: - return y_fit.params[0], 0. - else: + if len(y_fit.params) < 2: + return y_fit.params[0], 0.0 + else: return y_fit.params[:] else: return x_fit, y_fit.fittedvalues @@ -361,15 +457,16 @@ def linefit(x, y, return_coef=False): def is_empty(ifile): """If file is empty/corruted, return True.""" try: - with h5py.File(ifile, 'r') as f: return not bool(f.keys()) - except: + with h5py.File(ifile, "r") as f: + return not bool(list(f.keys())) + except IOError: return True """ Helper functions """ -def filter_data(t, h, bs, lew, tes, n_sigma=3, window=3/12.): +def filter_data(t, h, bs, lew, tes, n_sigma=3, window=3 / 12.0): """ Use various filters to remove outliers. @@ -408,35 +505,35 @@ def interp_params(t, h, bs, lew, tes): npts = [sum(~np.isnan(x)) for x in params] # Determine all the entries that should have valid data - isvalid = ~np.isnan(h) + isvalid = ~np.isnan(h) n_valid = sum(isvalid) # Do nothing if params are empty or have the same valid entries if np.all(npts == n_valid): - return params + return params # Sort indices for interpolation i_sort = np.argsort(t) for k, (n_p, p) in enumerate(zip(npts, [bs, lew, tes])): - + if n_p == n_valid: continue # Get the points that should be interpolated (if any) - i_interp, = np.where(np.isnan(p) & isvalid) + (i_interp,) = np.where(np.isnan(p) & isvalid) - ''' + """ print 'full length: ', len(p) print 'max valid entries:', n_valid print 'valid entries: ', sum(~np.isnan(p)) print 'interp entries: ', len(i_interp) print t[i_interp] - ''' + """ # Get sorted and clean (w/o NaNs) series to interpolate tt, pp = t[i_sort], p[i_sort] - i_valid, = np.where(~np.isnan(pp)) + (i_valid,) = np.where(~np.isnan(pp)) tt, pp = tt[i_valid], pp[i_valid] p_interp = np.interp(t[i_interp], tt, pp) @@ -448,13 +545,13 @@ def interp_params(t, h, bs, lew, tes): return params -def get_bboxs(lon, lat, dxy, proj='3031'): +def get_bboxs(lon, lat, dxy, proj="3031"): """ Define cells (bbox) for estimating corrections. """ # Convert into sterographic coordinates - x, y = transform_coord('4326', proj, lon, lat) + x, y = transform_coord("4326", proj, lon, lat) - # Number of tile edges on each dimension + # Number of tile edges on each dimension Nns = int(np.abs(np.nanmax(y) - np.nanmin(y)) / dxy) + 1 New = int(np.abs(np.nanmax(x) - np.nanmin(x)) / dxy) + 1 @@ -463,11 +560,14 @@ def get_bboxs(lon, lat, dxy, proj='3031'): yg = np.linspace(y.min(), y.max(), Nns) # Vector of bbox for each cell - bboxs = [(w,e,s,n) for w,e in zip(xg[:-1], xg[1:]) - for s,n in zip(yg[:-1], yg[1:])] + bboxs = [ + (w, e, s, n) + for w, e in zip(xg[:-1], xg[1:]) + for s, n in zip(yg[:-1], yg[1:]) + ] del xg, yg - #print 'total grid cells:', len(bboxs) + # print 'total grid cells:', len(bboxs) return bboxs @@ -485,50 +585,51 @@ def get_cell_idx(lon, lat, bbox, proj=3031): # Bounding box of grid cell xmin, xmax, ymin, ymax = bbox - + # Convert lon/lat to sterographic coordinates x, y = transform_coord(4326, proj, lon, lat) # Get the sub-tile (grid-cell) indices - i_cell, = np.where( (x >= xmin) & (x <= xmax) & - (y >= ymin) & (y <= ymax) ) + (i_cell,) = np.where((x >= xmin) & (x <= xmax) & (y >= ymin) & (y <= ymax)) return i_cell -def get_radius_idx(x, y, x0, y0, r, Tree, n_reloc=0): #NOTE: Add min reloc dist????? +def get_radius_idx( + x, y, x0, y0, r, Tree, n_reloc=0 +): # NOTE: Add min reloc dist????? """ Get indices of all data points inside radius. """ # Query the Tree from the node idx = Tree.query_ball_point((x0, y0), r) - ###print 'query #: 1 ( first search )' + # print 'query #: 1 ( first search )' # Either no relocation or not enough points to do relocation if n_reloc < 1 or len(idx) < 2: return idx - # Relocate center of search radius and query again + # Relocate center of search radius and query again for k in range(n_reloc): # Compute new search location => relocate initial center x0_new, y0_new = np.median(x[idx]), np.median(y[idx]) # Compute relocation distance - reloc_dist = np.hypot(x0_new-x0, y0_new-y0) + reloc_dist = np.hypot(x0_new - x0, y0_new - y0) # Do not allow total relocation to be larger than the search radius if reloc_dist > r: break - ###print 'query #:', k+2, '( reloc #:', k+1, ')' - ###print 'relocation dist:', reloc_dist + # print 'query #:', k+2, '( reloc #:', k+1, ')' + # print 'relocation dist:', reloc_dist # Query from the new location idx = Tree.query_ball_point((x0_new, y0_new), r) # If max number of relocations reached, exit - if n_reloc == k+1: + if n_reloc == k + 1: break return idx @@ -559,7 +660,7 @@ def multi_fit_coef(t_, h_, bs_, lew_, tes_): lew: leading-edge width tes: trailing-edge slope - """ + """ # Ensure zero mean of processed series h_, bs_, lew_, tes_ = center(h_, bs_, lew_, tes_) @@ -568,12 +669,12 @@ def multi_fit_coef(t_, h_, bs_, lew_, tes_): # Check for division by zero try: - + # Fit robust linear model on differenced series (w/o NaNs) - model = sm.OLS(h_, A_, missing="drop").fit(method='qr') + model = sm.OLS(h_, A_, missing="drop").fit(method="qr") + + # print model.summary() - #print model.summary() - # Get multivar coefficients for Bs, LeW, TeS a, b, c = model.params[:3] @@ -590,13 +691,13 @@ def multi_fit_coef(t_, h_, bs_, lew_, tes_): pvals = model.pvalues # Set all params to zero if exception detected - except: + except IOError: + + print("MULTIVARIATE FIT FAILED, setting params -> 0") + print("VALID DATA POINTS in h:", sum(~np.isnan(h_))) - print 'MULTIVARIATE FIT FAILED, setting params -> 0' - print 'VALID DATA POINTS in h:', sum(~np.isnan(h_)) - a, b, c, r2, pval, pvals = 0, 0, 0, 0, 1e3, [1e3, 1e3, 1e3] - + return [a, b, c, r2, pval, pvals] @@ -610,31 +711,33 @@ def rmse(t, x1, x2, order=1): def std_change(t, x1, x2, order=1): - """ Variance change from (detrended) x1 to x2 (magnitude and percentage). """ + """Variance change from (detrended) x1 to x2 (magnitude and percentage)""" idx = ~np.isnan(x1) & ~np.isnan(x2) - if sum(idx) < 3: return np.nan, np.nan + if sum(idx) < 3: + return np.nan, np.nan x1_res = detrend_binned(t, x1, order=order)[0] # use OLS poly fit x2_res = detrend_binned(t, x2, order=order)[0] s1 = mad_std(x1_res) s2 = mad_std(x2_res) delta_s = s2 - s1 - return delta_s, delta_s/s1 + return delta_s, delta_s / s1 def trend_change(t, x1, x2): """ Linear-trend change from x1 to x2 (magnitude and percentage). """ idx = ~np.isnan(x1) & ~np.isnan(x2) - if sum(idx) < 3: return np.nan, np.nan + if sum(idx) < 3: + return np.nan, np.nan t_, x1_, x2_ = t[idx], x1[idx], x2[idx] x1_ -= x1_.mean() x2_ -= x2_.mean() a1 = np.polyfit(t_, x1_, 1)[0] # use OLS poly fit a2 = np.polyfit(t_, x2_, 1)[0] delta_a = a2 - a1 - return delta_a, delta_a/np.abs(a1) + return delta_a, delta_a / np.abs(a1) -def sgolay1d(h, window=3, order=1, deriv=0, dt=1.0, mode='nearest', time=None): +def sgolay1d(h, window=3, order=1, deriv=0, dt=1.0, mode="nearest", time=None): """Savitztky-Golay filter with support for NaNs If time is given, interpolate NaNs otherwise pad w/zeros. @@ -642,8 +745,8 @@ def sgolay1d(h, window=3, order=1, deriv=0, dt=1.0, mode='nearest', time=None): dt is spacing between samples. """ h2 = h.copy() - ii, = np.where(np.isnan(h2)) - jj, = np.where(np.isfinite(h2)) + (ii,) = np.where(np.isnan(h2)) + (jj,) = np.where(np.isfinite(h2)) if len(ii) > 0 and time is not None: h2[ii] = np.interp(time[ii], time[jj], h2[jj]) elif len(ii) > 0: @@ -651,7 +754,7 @@ def sgolay1d(h, window=3, order=1, deriv=0, dt=1.0, mode='nearest', time=None): else: pass h2 = savgol_filter(h2, window, order, deriv, delta=dt, mode=mode) - return h2 + return h2 def overlap(x1, x2, y1, y2): @@ -661,24 +764,52 @@ def overlap(x1, x2, y1, y2): def intersect(x1, x2, y1, y2, a1, a2, b1, b2): """ Return True if two (x1,x2,y1,y2) rectangles intersect. """ - return (overlap(x1, x2, a1, a2) & overlap(y1, y2, b1, b2)) - - -def plot(x, y, xc, yc, tc, hc, bc, wc, sc, - hc_, bc_, wc_, sc_, hc_cor, h_bs, - r_bc, r_wc, r_sc, s_bc, s_wc, s_sc, - d_std, p_std, d_trend, p_trend, r2, pval, pvals): + return overlap(x1, x2, a1, a2) & overlap(y1, y2, b1, b2) + + +def plot( + x, + y, + xc, + yc, + tc, + hc, + bc, + wc, + sc, + hc_, + bc_, + wc_, + sc_, + hc_cor, + h_bs, + r_bc, + r_wc, + r_sc, + s_bc, + s_wc, + s_sc, + d_std, + p_std, + d_trend, + p_trend, + r2, + pval, + pvals, +): tc_ = tc.copy() # Bin variables - hc_b = binning(tc, hc_cor, median=True, window=3/12.)[1] ##NOTE: If binning before, this is binning over a binning - bc_b = binning(tc, bc, median=True, window=3/12.)[1] - wc_b = binning(tc, wc, median=True, window=3/12.)[1] - tc_b, sc_b = binning(tc, sc, median=True, window=3/12.)[:2] + hc_b = binning(tc, hc_cor, median=True, window=3 / 12.0)[ + 1 + ] # NOTE: If binning before, this is binning over a binning + bc_b = binning(tc, bc, median=True, window=3 / 12.0)[1] + wc_b = binning(tc, wc, median=True, window=3 / 12.0)[1] + tc_b, sc_b = binning(tc, sc, median=True, window=3 / 12.0)[:2] # Compute trends for plot - ii, = np.where(~np.isnan(tc) & ~np.isnan(hc) & ~np.isnan(hc_cor)) + (ii,) = np.where(~np.isnan(tc) & ~np.isnan(hc) & ~np.isnan(hc_cor)) t_, h1_, h2_ = tc[ii], hc[ii], hc_cor[ii] coefs1 = np.polyfit(t_, h1_, 1) coefs2 = np.polyfit(t_, h2_, 1) @@ -687,41 +818,45 @@ def plot(x, y, xc, yc, tc, hc, bc, wc, sc, # Mask NaNs for plotting mask = np.isfinite(hc_b) - + # Default color cycle cmap = plt.get_cmap("tab10") - plt.figure(figsize=(6,8)) - - plt.subplot(4,1,1) - plt.plot(tc, hc, '.') - plt.plot(tc, hc_cor, '.') - plt.plot(tc_b[mask], hc_b[mask], '-', color=cmap(3), linewidth=2) - plt.plot(tc, trend1, '-', color=cmap(0), linewidth=1.5) - plt.plot(tc, trend2, '-', color=cmap(3), linewidth=1.5) - plt.ylabel('Height (m)') - plt.title('Original time series') - - plt.subplot(4,1,2) - plt.plot(tc, bc, '.') - plt.plot(tc_b[mask], bc_b[mask], '-', linewidth=2) - plt.ylabel('Bs (s.d.)') - - plt.subplot(4,1,3) - plt.plot(tc, wc, '.') - plt.plot(tc_b[mask], wc_b[mask], '-', linewidth=2) - plt.ylabel('LeW (s.d.)') - - plt.subplot(4,1,4) - plt.plot(tc, sc, '.') - plt.plot(tc_b[mask], sc_b[mask], '-', linewidth=2) - plt.ylabel('TeS (s.d.)') + plt.figure(figsize=(6, 8)) + + plt.subplot(4, 1, 1) + plt.plot(tc, hc, ".") + plt.plot(tc, hc_cor, ".") + plt.plot(tc_b[mask], hc_b[mask], "-", color=cmap(3), linewidth=2) + plt.plot(tc, trend1, "-", color=cmap(0), linewidth=1.5) + plt.plot(tc, trend2, "-", color=cmap(3), linewidth=1.5) + plt.ylabel("Height (m)") + plt.title("Original time series") + + plt.subplot(4, 1, 2) + plt.plot(tc, bc, ".") + plt.plot(tc_b[mask], bc_b[mask], "-", linewidth=2) + plt.ylabel("Bs (s.d.)") + + plt.subplot(4, 1, 3) + plt.plot(tc, wc, ".") + plt.plot(tc_b[mask], wc_b[mask], "-", linewidth=2) + plt.ylabel("LeW (s.d.)") + + plt.subplot(4, 1, 4) + plt.plot(tc, sc, ".") + plt.plot(tc_b[mask], sc_b[mask], "-", linewidth=2) + plt.ylabel("TeS (s.d.)") # Bin variables - hc_b_ = binning(tc_, hc_, median=False, window=3/12.)[1] ##NOTE: If original TS were binned, - bc_b_ = binning(tc_, bc_, median=False, window=3/12.)[1] ## this is binning over a binning - wc_b_ = binning(tc_, wc_, median=False, window=3/12.)[1] - tc_b_, sc_b_ = binning(tc_, sc_, median=True, window=3/12.)[:2] + hc_b_ = binning(tc_, hc_, median=False, window=3 / 12.0)[ + 1 + ] # NOTE: If original TS were binned, + bc_b_ = binning(tc_, bc_, median=False, window=3 / 12.0)[ + 1 + ] # this is binning over a binning + wc_b_ = binning(tc_, wc_, median=False, window=3 / 12.0)[1] + tc_b_, sc_b_ = binning(tc_, sc_, median=True, window=3 / 12.0)[:2] # Fit seasonality plot_season = True @@ -732,51 +867,54 @@ def plot(x, y, xc, yc, tc, hc, bc, wc, sc, # Trend + Seasonal model def func(t, c, m, n, a, p): """ Seasonality with amplitude and phase. """ - #return a * np.sin(2*np.pi * t + p) - return c + m * t + n * t**2 + a * np.sin(2*np.pi * t + p) + # return a * np.sin(2*np.pi * t + p) + return c + m * t + n * t ** 2 + a * np.sin(2 * np.pi * t + p) # Fit seasonality ii = ~np.isnan(tc_) & ~np.isnan(hc_) t_season, h_season = tc_[ii], hc_[ii] - params, params_covariance = optimize.curve_fit(func, t_season, h_season, p0=[0., 0., 0., 0.1, 0.1]) + params, params_covariance = optimize.curve_fit( + func, t_season, h_season, p0=[0.0, 0.0, 0.0, 0.1, 0.1] + ) a0, a1, a2, amp, pha = params t_season = np.linspace(t_season.min(), t_season.max(), 100) h_season = func(t_season, a0, a1, a2, amp, pha) - # mask NaNs for plotting mask = np.isfinite(hc_b_) - plt.figure(figsize=(6,8)) - - plt.subplot(4,1,1) - plt.plot(tc_, hc_, '.') - plt.plot(tc_b_[mask], hc_b_[mask], '-') - if plot_season: plt.plot(t_season, h_season, '-r') - plt.ylabel('Height (m)') - plt.title('Processed time series') - if plot_season: plt.title('Amplitude = %.2f, Phase = %.2f' % (amp, pha)) - - plt.subplot(4,1,2) - plt.plot(tc_, bc_, '.') - plt.plot(tc_b_[mask], bc_b_[mask], '-') - plt.ylabel('Bs (s.d.)') - - plt.subplot(4,1,3) - plt.plot(tc_, wc_, '.') - plt.plot(tc_b_[mask], wc_b_[mask], '-') - plt.ylabel('LeW (s.d.)') - - plt.subplot(4,1,4) - plt.plot(tc_, sc_, '.') - plt.plot(tc_b_[mask], sc_b_[mask], '-') - plt.ylabel('TeS (s.d.)') + plt.figure(figsize=(6, 8)) + + plt.subplot(4, 1, 1) + plt.plot(tc_, hc_, ".") + plt.plot(tc_b_[mask], hc_b_[mask], "-") + if plot_season: + plt.plot(t_season, h_season, "-r") + plt.ylabel("Height (m)") + plt.title("Processed time series") + if plot_season: + plt.title("Amplitude = %.2f, Phase = %.2f" % (amp, pha)) + + plt.subplot(4, 1, 2) + plt.plot(tc_, bc_, ".") + plt.plot(tc_b_[mask], bc_b_[mask], "-") + plt.ylabel("Bs (s.d.)") + + plt.subplot(4, 1, 3) + plt.plot(tc_, wc_, ".") + plt.plot(tc_b_[mask], wc_b_[mask], "-") + plt.ylabel("LeW (s.d.)") + + plt.subplot(4, 1, 4) + plt.plot(tc_, sc_, ".") + plt.plot(tc_b_[mask], sc_b_[mask], "-") + plt.ylabel("TeS (s.d.)") plt.figure() - plt.plot(x, y, '.', color='0.6', zorder=1) + plt.plot(x, y, ".", color="0.6", zorder=1) plt.scatter(xc, yc, c=hc_cor, s=5, vmin=-1, vmax=1, zorder=2) - plt.plot(np.nanmedian(xc), np.nanmedian(yc), 'o', color='red', zorder=3) - plt.title('Tracks') + plt.plot(np.nanmedian(xc), np.nanmedian(yc), "o", color="red", zorder=3) + plt.title("Tracks") # Plot Spectrum if 1: @@ -785,13 +923,13 @@ def plot(x, y, xc, yc, tc, hc, bc, wc, sc, from astropy.stats import LombScargle - periods = np.arange(3/12., 1.5, 0.01) - freq = 1/periods + periods = np.arange(3 / 12.0, 1.5, 0.01) + freq = 1 / periods - hc[np.isnan(hc)] = 0. - bc[np.isnan(bc)] = 0. - wc[np.isnan(wc)] = 0. - sc[np.isnan(sc)] = 0. + hc[np.isnan(hc)] = 0.0 + bc[np.isnan(bc)] = 0.0 + wc[np.isnan(wc)] = 0.0 + sc[np.isnan(sc)] = 0.0 ls = LombScargle(tc, hc, nterms=1) power_h = ls.power(freq) @@ -802,116 +940,131 @@ def plot(x, y, xc, yc, tc, hc, bc, wc, sc, ls = LombScargle(tc, sc, nterms=1) power_s = ls.power(freq) - plt.figure(figsize=(6,8)) + plt.figure(figsize=(6, 8)) - plt.subplot(4,1,1) - plt.plot(freq, power_h, linewidth=2) - plt.xlabel('Frequency (cycles/year)') - plt.ylabel('Power (1/RMSE)') - plt.title('Spectrum of processed time series') + plt.subplot(4, 1, 1) + plt.plot(freq, power_h, linewidth=2) + plt.xlabel("Frequency (cycles/year)") + plt.ylabel("Power (1/RMSE)") + plt.title("Spectrum of processed time series") - plt.subplot(4,1,2) - plt.plot(freq, power_b, linewidth=2) - plt.xlabel('Frequency (cycles/year)') - plt.ylabel('Power (1/RMSE)') + plt.subplot(4, 1, 2) + plt.plot(freq, power_b, linewidth=2) + plt.xlabel("Frequency (cycles/year)") + plt.ylabel("Power (1/RMSE)") - plt.subplot(4,1,3) - plt.plot(freq, power_w, linewidth=2) - plt.xlabel('Frequency (cycles/year)') - plt.ylabel('Power (1/RMSE)') + plt.subplot(4, 1, 3) + plt.plot(freq, power_w, linewidth=2) + plt.xlabel("Frequency (cycles/year)") + plt.ylabel("Power (1/RMSE)") - plt.subplot(4,1,4) - plt.plot(freq, power_s, linewidth=2) - plt.xlabel('Frequency (cycles/year)') - plt.ylabel('Power (1/RMSE)') + plt.subplot(4, 1, 4) + plt.plot(freq, power_s, linewidth=2) + plt.xlabel("Frequency (cycles/year)") + plt.ylabel("Power (1/RMSE)") # Plot Crosscorrelation if 1: tc, hc, bc, wc, sc = tc_, hc_, bc_, wc_, sc_ # processed TS - plt.figure(figsize=(6,8)) + plt.figure(figsize=(6, 8)) plt.subplot(311) plt.xcorr(hc, bc) - plt.ylabel('h x bs') - plt.title('Crosscorrelation') + plt.ylabel("h x bs") + plt.title("Crosscorrelation") plt.subplot(312) plt.xcorr(hc, wc) - plt.ylabel('h x LeW') + plt.ylabel("h x LeW") plt.subplot(313) plt.xcorr(hc, sc) - plt.ylabel('h x TeS') + plt.ylabel("h x TeS") tc, hc, bc, wc, sc = tc_, hc_, bc_, wc_, sc_ # processed TS - plt.figure(figsize=(3,9)) + plt.figure(figsize=(3, 9)) plt.subplot(311) - plt.plot(bc, hc, '.') - #plt.title('Correlation Bs x h (%s)' % str(proc)) - plt.xlabel('Bs (s.d.)') - plt.ylabel('h (m)') - plt.title('Correlation of processed time series') + plt.plot(bc, hc, ".") + # plt.title('Correlation Bs x h (%s)' % str(proc)) + plt.xlabel("Bs (s.d.)") + plt.ylabel("h (m)") + plt.title("Correlation of processed time series") plt.subplot(312) - plt.plot(wc, hc, '.') - #plt.title('Correlation Bs x h (%s)' % str(proc)) - plt.xlabel('LeW (s.d.)') - plt.ylabel('h (m)') + plt.plot(wc, hc, ".") + # plt.title('Correlation Bs x h (%s)' % str(proc)) + plt.xlabel("LeW (s.d.)") + plt.ylabel("h (m)") plt.subplot(313) - plt.plot(sc, hc, '.') - #plt.title('Correlation Bs x h (%s)' % str(proc)) - plt.xlabel('TeS (s.d.)') - plt.ylabel('h (m)') - - - print 'Summary:' - print '--------' - print 'cor applied: ', (h_bs[~np.isnan(h_bs)] != 0).any() - print 'std change: %.3f m (%.1f %%)' % (round(d_std, 3), round(p_std*100, 1)) - print 'trend change: %.3f m/yr (%.1f %%)' % (round(d_trend, 3), round(p_trend*100, 1)) - print '' - print 'r-squared: ', round(r2, 3) - print 'p-value: ', round(pval, 3) - print 'p-values: ', [round(p, 3) for p in pvals] - print '' - print 'r_bs: ', round(r_bc, 3) - print 'r_lew: ', round(r_wc, 3) - print 'r_tes: ', round(r_sc, 3) - print '' - print 's_bs: ', round(s_bc, 3) - print 's_lew: ', round(s_wc, 3) - print 's_tes: ', round(s_sc, 3) + plt.plot(sc, hc, ".") + # plt.title('Correlation Bs x h (%s)' % str(proc)) + plt.xlabel("TeS (s.d.)") + plt.ylabel("h (m)") + + print("Summary:") + print("--------") + print("cor applied: ", (h_bs[~np.isnan(h_bs)] != 0).any()) + print( + "std change: %.3f m (%.1f %%)" + % (round(d_std, 3), round(p_std * 100, 1)) + ) + print( + "trend change: %.3f m/yr (%.1f %%)" + % (round(d_trend, 3), round(p_trend * 100, 1)) + ) + print("") + print("r-squared: ", round(r2, 3)) + print("p-value: ", round(pval, 3)) + print("p-values: ", [round(p, 3) for p in pvals]) + print("") + print("r_bs: ", round(r_bc, 3)) + print("r_lew: ", round(r_wc, 3)) + print("r_tes: ", round(r_sc, 3)) + print("") + print("s_bs: ", round(s_bc, 3)) + print("s_lew: ", round(s_wc, 3)) + print("s_tes: ", round(s_sc, 3)) plt.show() -def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply_=False): +def main( + ifile, + vnames, + wnames, + dxy, + proj, + radius=0, + n_reloc=0, + proc=None, + apply_=False, +): if is_empty(ifile): - print 'empty file... skipping!!!' + print("empty file... skipping!!!") return - + if TEST_MODE: - print '*********************************************************' - print '* RUNNING IN TEST MODE (PLOTTING ONLY, NOT SAVING DATA) *' - print '*********************************************************' + print("*********************************************************") + print("* RUNNING IN TEST MODE (PLOTTING ONLY, NOT SAVING DATA) *") + print("*********************************************************") + + print("processing file:", ifile, "...") - print 'processing file:', ifile, '...' - # Test if parameter file exists - if '_scatgrd' in ifile.lower(): + if "_scatgrd" in ifile.lower(): return xvar, yvar, zvar, tvar = vnames bpar, wpar, spar = wnames # Load full data into memory (only once) - with h5py.File(ifile, 'r') as fi: + with h5py.File(ifile, "r") as fi: t = fi[tvar][:] h = fi[zvar][:] @@ -922,7 +1075,7 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply tes = fi[spar][:] # Convert into sterographic coordinates - x, y = transform_coord('4326', proj, lon, lat) + x, y = transform_coord("4326", proj, lon, lat) # Get bbox from data xmin_d, xmax_d, ymin_d, ymax_d = x.min(), x.max(), y.min(), y.max() @@ -945,77 +1098,92 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply N_nodes = len(x_nodes) # Values for each data point - r2fit = np.full(N_data, 0.0) # r2 of the multivar fit - pval = np.full(N_data, np.nan) # r2 of the multivar fit - dstd = np.full(N_data, np.nan) # magnitude std change after cor - dtrend = np.full(N_data, np.nan) # magnitude trend change after cor - pstd = np.full(N_data, np.inf) # perc std change after cor ##NOTE: Init w/inf - ptrend = np.full(N_data, np.nan) # perc trend change after cor - hbs = np.full(N_data, np.nan) # scatt cor from multivar fit - rbs = np.full(N_data, np.nan) # corr coef h x Bs - rlew = np.full(N_data, np.nan) # corr coef h x LeW - rtes = np.full(N_data, np.nan) # corr coef h x TeS - sbs = np.full(N_data, np.nan) # sensit h x Bs - slew = np.full(N_data, np.nan) # sensit h x LeW - stes = np.full(N_data, np.nan) # sensit h x TeS - bbs = np.full(N_data, np.nan) # multivar fit coef a.Bs - blew = np.full(N_data, np.nan) # multivar fit coef b.LeW - btes = np.full(N_data, np.nan) # multivar fit coef c.TeS + r2fit = np.full(N_data, 0.0) # r2 of the multivar fit + pval = np.full(N_data, np.nan) # r2 of the multivar fit + dstd = np.full(N_data, np.nan) # magnitude std change after cor + dtrend = np.full(N_data, np.nan) # magnitude trend change after cor + pstd = np.full( + N_data, np.inf + ) # perc std change after cor ##NOTE: Init w/inf + ptrend = np.full(N_data, np.nan) # perc trend change after cor + hbs = np.full(N_data, np.nan) # scatt cor from multivar fit + rbs = np.full(N_data, np.nan) # corr coef h x Bs + rlew = np.full(N_data, np.nan) # corr coef h x LeW + rtes = np.full(N_data, np.nan) # corr coef h x TeS + sbs = np.full(N_data, np.nan) # sensit h x Bs + slew = np.full(N_data, np.nan) # sensit h x LeW + stes = np.full(N_data, np.nan) # sensit h x TeS + bbs = np.full(N_data, np.nan) # multivar fit coef a.Bs + blew = np.full(N_data, np.nan) # multivar fit coef b.LeW + btes = np.full(N_data, np.nan) # multivar fit coef c.TeS # Values for each node r2fitc = np.full(N_nodes, 0.0) - pvalc = np.full(N_nodes, np.nan) dstdc = np.full(N_nodes, np.nan) dtrendc = np.full(N_nodes, np.nan) - pstdc = np.full(N_nodes, np.inf) ##NOTE: Init w/inf + pstdc = np.full(N_nodes, np.inf) # NOTE: Init w/inf ptrendc = np.full(N_nodes, np.nan) rbsc = np.full(N_nodes, np.nan) - rlewc = np.full(N_nodes, np.nan) - rtesc = np.full(N_nodes, np.nan) - sbsc = np.full(N_nodes, np.nan) - slewc = np.full(N_nodes, np.nan) - stesc = np.full(N_nodes, np.nan) - bbsc = np.full(N_nodes, np.nan) - blewc = np.full(N_nodes, np.nan) - btesc = np.full(N_nodes, np.nan) - lonc = np.full(N_nodes, np.nan) - latc = np.full(N_nodes, np.nan) + rlewc = np.full(N_nodes, np.nan) + rtesc = np.full(N_nodes, np.nan) + sbsc = np.full(N_nodes, np.nan) + slewc = np.full(N_nodes, np.nan) + stesc = np.full(N_nodes, np.nan) + bbsc = np.full(N_nodes, np.nan) + blewc = np.full(N_nodes, np.nan) + btesc = np.full(N_nodes, np.nan) + lonc = np.full(N_nodes, np.nan) + latc = np.full(N_nodes, np.nan) # Select cells at random (for testing) if TEST_MODE: if USE_NODES: # Convert into sterographic coordinates - x_nodes = [transform_coord('4326', '3031', xp, yp)[0] for xp, yp in NODES] - y_nodes = [transform_coord('4326', '3031', xp, yp)[1] for xp, yp in NODES] + x_nodes = [ + transform_coord("4326", "3031", xp, yp)[0] for xp, yp in NODES + ] + y_nodes = [ + transform_coord("4326", "3031", xp, yp)[1] for xp, yp in NODES + ] else: - if USE_SEED: np.random.seed(SEED) # not so random! + if USE_SEED: + np.random.seed(SEED) # not so random! # Select a few random nodes ii = np.random.randint(0, N_nodes, N_CELLS) - x_nodes, y_nodes = x_nodes[ii], y_nodes[ii] + x_nodes, y_nodes = x_nodes[ii], y_nodes[ii] N_nodes = len(x_nodes) # Build KD-Tree with polar stereo coords x, y = transform_coord(4326, proj, lon, lat) - Tree = cKDTree(zip(x, y)) + Tree = cKDTree(list(zip(x, y))) # Loop through nodes - for k in xrange(N_nodes): + for k in range(N_nodes): - if (k%500) == 0: - print 'Calculating correction for node', k, 'of', N_nodes, '...' + if (k % 500) == 0: + print("Calculating correction for node", k, "of", N_nodes, "...") x0, y0 = x_nodes[k], y_nodes[k] # If search radius doesn't contain data, skip - if not intersect(x0-radius, x0+radius, y0-radius, y0+radius, - xmin_d, xmax_d, ymin_d, ymax_d): + if not intersect( + x0 - radius, + x0 + radius, + y0 - radius, + y0 + radius, + xmin_d, + xmax_d, + ymin_d, + ymax_d, + ): continue # Get indices of data within search radius i_cell = get_radius_idx(x, y, x0, y0, radius, Tree, n_reloc=n_reloc) # If cell empty or not enough data go to next node - if len(i_cell) < MIN_PTS: continue + if len(i_cell) < MIN_PTS: + continue # Get all data within the grid search radius tc = t[i_cell] @@ -1027,73 +1195,101 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply sc = tes[i_cell] # Keep original (unfiltered) data - tc_orig, hc_orig, bc_orig, wc_orig, sc_orig = tc.copy(), hc.copy(), bc.copy(), wc.copy(), sc.copy() + tc_orig, hc_orig = tc.copy(), hc.copy() PLOT_INTERP = False - if PLOT_INTERP: bc0, wc0, sc0 = bc.copy(), wc.copy(), sc.copy() #NOTE: for plotting + if PLOT_INTERP: + bc0, wc0, sc0 = ( + bc.copy(), + wc.copy(), + sc.copy(), + ) # NOTE: for plotting # Filter invalid points - tc, hc, bc, wc, sc = filter_data(tc, hc, bc, wc, sc, n_sigma=3, window=3/12.) + tc, hc, bc, wc, sc = filter_data( + tc, hc, bc, wc, sc, n_sigma=3, window=3 / 12.0 + ) - if PLOT_INTERP: bc1, wc1, sc1 = bc.copy(), wc.copy(), sc.copy() #NOTE: for plotting + if PLOT_INTERP: + bc1, wc1, sc1 = ( + bc.copy(), + wc.copy(), + sc.copy(), + ) # NOTE: for plotting # Test minimum number of obs in all params nobs = min([len(v[~np.isnan(v)]) for v in [hc, bc, wc, sc]]) # Test for enough points - if (nobs < MIN_PTS): continue + if nobs < MIN_PTS: + continue # Bin at monthly intervals to check temporal sampling - h_bin = binning(tc, hc, dx=1/12., window=3/12., interp=False)[1] - if sum(~np.isnan(h_bin)) < MIN_MONTHS: continue + h_bin = binning(tc, hc, dx=1 / 12.0, window=3 / 12.0, interp=False)[1] + if sum(~np.isnan(h_bin)) < MIN_MONTHS: + continue # Interpolate missing w/f params based on h series - #bc, wc, sc = interp_params(tc, hc, bc, wc, sc) #NOTE: Leave this out for now (only use valid entries in ALL series) <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Is this a good idea? + # bc, wc, sc = interp_params(tc, hc, bc, wc, sc) + # NOTE: Leave this out for now (only use valid entries in ALL series) - if PLOT_INTERP: bc2, wc2, sc2 = bc.copy(), wc.copy(), sc.copy() #NOTE: for plotting + if PLOT_INTERP: + bc2, wc2, sc2 = ( + bc.copy(), + wc.copy(), + sc.copy(), + ) # NOTE: for plotting # Plot interpolated time-series points if PLOT_INTERP: - bc2[bc2==bc1] = np.nan - wc2[wc2==wc1] = np.nan - sc2[sc2==sc1] = np.nan - - ax0 = plt.subplot(411) - plt.plot(tc, hc_orig, '.') - plt.plot(tc, hc, '.') - ax1 = plt.subplot(412) - plt.plot(tc, bc0, '.') - plt.plot(tc, bc1, '.') - plt.plot(tc, bc2, 'x') - ax2 = plt.subplot(413) - plt.plot(tc, wc0, '.') - plt.plot(tc, wc1, '.') - plt.plot(tc, wc2, 'x') - ax3 = plt.subplot(414) - plt.plot(tc, sc0, '.') - plt.plot(tc, sc1, '.') - plt.plot(tc, sc2, 'x') + bc2[bc2 == bc1] = np.nan + wc2[wc2 == wc1] = np.nan + sc2[sc2 == sc1] = np.nan + + plt.subplot(411) + plt.plot(tc, hc_orig, ".") + plt.plot(tc, hc, ".") + plt.subplot(412) + plt.plot(tc, bc0, ".") + plt.plot(tc, bc1, ".") + plt.plot(tc, bc2, "x") + plt.subplot(413) + plt.plot(tc, wc0, ".") + plt.plot(tc, wc1, ".") + plt.plot(tc, wc2, "x") + plt.subplot(414) + plt.plot(tc, sc0, ".") + plt.plot(tc, sc1, ".") + plt.plot(tc, sc2, "x") plt.figure() - plt.plot(x, y, '.', color='.5', rasterized=True) - plt.plot(xc, yc, '.') - plt.plot(xc[~np.isnan(hc)], yc[~np.isnan(hc)], '.') + plt.plot(x, y, ".", color=".5", rasterized=True) + plt.plot(xc, yc, ".") + plt.plot(xc[~np.isnan(hc)], yc[~np.isnan(hc)], ".") plt.show() continue if BIN_SERIES: - hc_bin = binning(tc, hc, median=True, window=3/12., interp=True)[1] - bc_bin = binning(tc, bc, median=True, window=3/12., interp=True)[1] - wc_bin = binning(tc, wc, median=True, window=3/12., interp=True)[1] - sc_bin = binning(tc, sc, median=True, window=3/12., interp=True)[1] + hc_bin = binning( + tc, hc, median=True, window=3 / 12.0, interp=True + )[1] + bc_bin = binning( + tc, bc, median=True, window=3 / 12.0, interp=True + )[1] + wc_bin = binning( + tc, wc, median=True, window=3 / 12.0, interp=True + )[1] + sc_bin = binning( + tc, sc, median=True, window=3 / 12.0, interp=True + )[1] else: hc_bin = hc bc_bin = bc wc_bin = wc sc_bin = sc - + # Ensure zero mean on all variables hc, bc, wc, sc = center(hc, bc, wc, sc) hc_bin, bc_bin, wc_bin, sc_bin = center(hc_bin, bc_bin, wc_bin, sc_bin) @@ -1102,7 +1298,7 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply bc, wc, sc = normalize(bc, wc, sc) bc_bin, wc_bin, sc_bin = normalize(bc_bin, wc_bin, sc_bin) - if proc == 'det': + if proc == "det": # Detrend time series hc_res = detrend_binned(tc, hc_bin, order=2)[0] bc_res = detrend_binned(tc, bc_bin, order=2)[0] @@ -1110,24 +1306,31 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply sc_res = detrend_binned(tc, sc_bin, order=2)[0] else: # Savitzky-Golay numerical diff - hc_res = sgolay1d(hc_bin, WINDOW, ORDER, DERIV) #FIXME: Think what dt to use here (uneven scattered points)!!!!!!! + hc_res = sgolay1d( + hc_bin, WINDOW, ORDER, DERIV + ) # FIXME: Think what dt to use here (uneven scattered points) bc_res = sgolay1d(bc_bin, WINDOW, ORDER, DERIV) wc_res = sgolay1d(wc_bin, WINDOW, ORDER, DERIV) sc_res = sgolay1d(sc_bin, WINDOW, ORDER, DERIV) # Get coefs from multivariate fit - b_bc, b_wc, b_sc, r2, pval, pvals = multi_fit_coef(tc, hc_res, bc_res, wc_res, sc_res) + b_bc, b_wc, b_sc, r2, pval, pvals = multi_fit_coef( + tc, hc_res, bc_res, wc_res, sc_res + ) - if sum([b_bc, b_wc, b_sc]) == 0: continue + if sum([b_bc, b_wc, b_sc]) == 0: + continue - # Get linear combination of original FILTERED series => h_bs = a Bs + b LeW + c TeS + # Get linear combination of original FILTERED series + # h_bs = a Bs + b LeW + c TeS hc_bs = np.dot(np.vstack((bc, wc, sc)).T, [b_bc, b_wc, b_sc]) - #NOTE 1: The correction is generated using the original w/f series - #NOTE 2: Use np.dot so NaNs from original (filtered) series are preserved + # NOTE 1: The correction is generated using the original w/f series + # NOTE 2: Use np.dot to preseved NaNs from original (filtered) series + + if np.isnan(hc_bs).all(): + continue - if np.isnan(hc_bs).all(): continue - # Apply correction to height hc_cor = hc - hc_bs @@ -1135,7 +1338,9 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply r_bc, r_wc, r_sc = corr_coef(hc_res, bc_res, wc_res, sc_res) # Calculate sensitivity values (corr grad) - s_bc, s_wc, s_sc = corr_grad(hc_res, bc_res, wc_res, sc_res, normalize=False) + s_bc, s_wc, s_sc = corr_grad( + hc_res, bc_res, wc_res, sc_res, normalize=False + ) # Calculate variance change (magnitude and perc) d_std, p_std = std_change(tc, hc, hc_cor, order=1) @@ -1143,75 +1348,107 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply # Calculate trend change (magnitude and perc) d_trend, p_trend = trend_change(tc, hc, hc_cor) - if np.isnan([d_std, p_std, d_trend, p_std]).any(): continue + if np.isnan([d_std, p_std, d_trend, p_std]).any(): + continue # Calculate RMSE between detrended series - #d_rmse = rmse(tc, hc, hc_cor, order=1) - + # d_rmse = rmse(tc, hc, hc_cor, order=1) + # Test if at least one correlation is significant - #r_cond = (np.abs(r_bc) < R_MIN and np.abs(r_wc) < R_MIN and np.abs(r_sc) < R_MIN) + # r_cond = (np.abs(r_bc) < R_MIN and np.abs(r_wc) \ + # < R_MIN and np.abs(r_sc) < R_MIN) - ##NOTE: We are ignoring the 'pval' (significance of fit). - ##NOTE: All we care about is the reduction in variance. + # NOTE: We are ignoring the 'pval' (significance of fit). + # NOTE: All we care about is the reduction in variance. # Do not apply correction if: std increases by more than 5% if p_std > 0.05: # Cor is set to zero (keeping NaNs) - hc_cor = hc.copy() # back to original - hc_bs[~np.isnan(hc_bs)] = 0. # hc_bs keeps NaNs from filtered out values + hc_cor = hc.copy() # back to original + hc_bs[ + ~np.isnan(hc_bs) + ] = 0.0 # hc_bs keeps NaNs from filtered out values - ''' + """ # All params are set to zero/one b_bc, b_wc, b_sc = 0., 0., 0. r_bc, r_wc, r_sc = 0., 0., 0. s_bc, s_wc, s_sc = 0., 0., 0. r2, pval, pvals = 0., 1e3, (1e3, 1e3, 1e3) d_std, p_std, d_trend, p_trend = 0., 0., 0., 0. - ''' - + """ + # Plot individual grid cells for testing if TEST_MODE: - plt.figure(figsize=(6,2)) - plt.plot(tc_orig, hc_orig, '.', color='0.3') - plot(x, y, xc, yc, tc, hc, bc, wc, sc, - hc_res, bc_res, wc_res, sc_res, hc_cor, hc_bs, - r_bc, r_wc, r_sc, s_bc, s_wc, s_sc, - d_std, p_std, d_trend, p_trend, - r2, pval, pvals) + plt.figure(figsize=(6, 2)) + plt.plot(tc_orig, hc_orig, ".", color="0.3") + plot( + x, + y, + xc, + yc, + tc, + hc, + bc, + wc, + sc, + hc_res, + bc_res, + wc_res, + sc_res, + hc_cor, + hc_bs, + r_bc, + r_wc, + r_sc, + s_bc, + s_wc, + s_sc, + d_std, + p_std, + d_trend, + p_trend, + r2, + pval, + pvals, + ) """ Store results (while checking previously stored estimates) """ - # Check where previously stored values need update - #i_update, = np.where(r2fit[i_cell] <= r2) # r2_prev <= r2_new (new r2 is larger) - i_update, = np.where(pstd[i_cell] > p_std) # std_prev > std_new (new std is lower) + # Check where previously stored values need update + # r2_prev <= r2_new (new r2 is larger) + # i_update, = np.where(r2fit[i_cell] <= r2) + (i_update,) = np.where( + pstd[i_cell] > p_std + ) # std_prev > std_new (new std is lower) - # Remove indices and values that will not be updated + # Remove indices and values that will not be updated i_cell_new = [i_cell[i] for i in i_update] hc_bs_new = hc_bs[i_update] # Store correction for cell: only the improved values - hbs[i_cell_new] = hc_bs_new # set of values - r2fit[i_cell_new] = r2 # one value (same for all) - dstd[i_cell_new] = d_std - pstd[i_cell_new] = p_std - dtrend[i_cell_new] = d_trend - ptrend[i_cell_new] = p_trend - rbs[i_cell_new] = r_bc # corr coef - rlew[i_cell_new] = r_wc + hbs[i_cell_new] = hc_bs_new # set of values + r2fit[i_cell_new] = r2 # one value (same for all) + dstd[i_cell_new] = d_std + pstd[i_cell_new] = p_std + dtrend[i_cell_new] = d_trend + ptrend[i_cell_new] = p_trend + rbs[i_cell_new] = r_bc # corr coef + rlew[i_cell_new] = r_wc rtes[i_cell_new] = r_sc - sbs[i_cell_new] = s_bc # sensitivity + sbs[i_cell_new] = s_bc # sensitivity slew[i_cell_new] = s_wc stes[i_cell_new] = s_sc - bbs[i_cell_new] = b_bc # multivar fit coef + bbs[i_cell_new] = b_bc # multivar fit coef blew[i_cell_new] = b_wc btes[i_cell_new] = b_sc - # Compute centroid of cell + # Compute centroid of cell xc_ = np.nanmedian(xc) yc_ = np.nanmedian(yc) - # Convert x/y -> lon/lat + # Convert x/y -> lon/lat lon_c, lat_c = transform_coord(proj, 4326, xc_, yc_) # Store one s and r value per cell @@ -1234,40 +1471,43 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply """ Correct h (full dataset) with best values """ - if apply_: h[~np.isnan(hbs)] -= hbs[~np.isnan(hbs)] + if apply_: + h[~np.isnan(hbs)] -= hbs[~np.isnan(hbs)] """ Save data """ if not TEST_MODE: - print 'saving data ...' + print("saving data ...") - with h5py.File(ifile, 'a') as fi: + with h5py.File(ifile, "a") as fi: # Update h in the file and save correction (all cells at once) - if apply_: fi[zvar][:] = h - + if apply_: + fi[zvar][:] = h + # Try to create varibales try: # Save params for each point fi[H_BS] = hbs - fi['r2'] = r2fit - fi['d_std'] = dstd - fi['p_std'] = pstd - fi['d_trend'] = dtrend - fi['p_trend'] = ptrend - fi['r_bs'] = rbs - fi['r_lew'] = rlew - fi['r_tes'] = rtes - fi['s_bs'] = sbs - fi['s_lew'] = slew - fi['s_tes'] = stes - fi['b_bs'] = bbs - fi['b_lew'] = blew - fi['b_tes'] = btes - - #FIXME: Check if this is a good idea. Content of input file is being deleted!!! Removing for now!!! - # Update variabels instead - except: + fi["r2"] = r2fit + fi["d_std"] = dstd + fi["p_std"] = pstd + fi["d_trend"] = dtrend + fi["p_trend"] = ptrend + fi["r_bs"] = rbs + fi["r_lew"] = rlew + fi["r_tes"] = rtes + fi["s_bs"] = sbs + fi["s_lew"] = slew + fi["s_tes"] = stes + fi["b_bs"] = bbs + fi["b_lew"] = blew + fi["b_tes"] = btes + + # FIXME: Check if this is a good idea. + # Content of input file is being deleted! Removing for now! + # Update variables instead + except IOError: """ # Save params for each point fi[H_BS][:] = hbs @@ -1286,76 +1526,84 @@ def main(ifile, vnames, wnames, dxy, proj, radius=0, n_reloc=0, proc=None, apply fi['b_lew'][:] = blew fi['b_tes'][:] = btes """ - print 'SOME PROBLEM WITH THE FILE' - print 'PARAMETERS NOT SAVED!' - print ifile + print("SOME PROBLEM WITH THE FILE") + print("PARAMETERS NOT SAVED!") + print(ifile) return # Only rename file if _SCAT has not been added - if ifile.find(SUFFIX1+'.h5') < 0: - os.rename(ifile, ifile.replace('.h5', SUFFIX1+'.h5')) - - # Save bs params as external file - with h5py.File(ifile.replace('.h5', SUFFIX2+'.h5'), 'w') as fo: - + if ifile.find(SUFFIX1 + ".h5") < 0: + os.rename(ifile, ifile.replace(".h5", SUFFIX1 + ".h5")) + + # Save bs params as external file + with h5py.File(ifile.replace(".h5", SUFFIX2 + ".h5"), "w") as fo: + # Try to svave variables try: # Save varibales - fo['lon'] = lonc - fo['lat'] = latc - fo['r2'] = r2fitc - fo['d_std'] = dstdc - fo['p_std'] = pstdc - fo['d_trend'] = dtrendc - fo['p_trend'] = ptrendc - fo['r_bs'] = rbsc - fo['r_lew'] = rlewc - fo['r_tes'] = rtesc - fo['s_bs'] = sbsc - fo['s_lew'] = slewc - fo['s_tes'] = stesc - fo['b_bs'] = bbsc - fo['b_lew'] = blewc - fo['b_tes'] = btesc - - # Catch any exceptions - except: + fo["lon"] = lonc + fo["lat"] = latc + fo["r2"] = r2fitc + fo["d_std"] = dstdc + fo["p_std"] = pstdc + fo["d_trend"] = dtrendc + fo["p_trend"] = ptrendc + fo["r_bs"] = rbsc + fo["r_lew"] = rlewc + fo["r_tes"] = rtesc + fo["s_bs"] = sbsc + fo["s_lew"] = slewc + fo["s_tes"] = stesc + fo["b_bs"] = bbsc + fo["b_lew"] = blewc + fo["b_tes"] = btesc + + # Catch any exceptions + except IOError: # Exit program - print 'COUND NOT SAVE PARAMETERS FOR EACH CELL (SCATGRD)' + print("COUND NOT SAVE PARAMETERS FOR EACH CELL (SCATGRD)") return -if __name__ == '__main__': +if __name__ == "__main__": - # Pass arguments + # Pass arguments args = get_args() - ifiles = args.files[:] # input files - vnames = args.vnames[:] # lon/lat/h/time variable names - wnames = args.wnames[:] # variables to use for correction - dxy = args.dxy[0] * 1e3 # grid-cell length (km -> m) - radius = args.radius[0] * 1e3 # search radius (km -> m) - nreloc = args.nreloc[0] # number of relocations - proc = args.proc[0] # det, dif, bin or None series - proj = args.proj[0] # EPSG proj number - apply_ = args.apply # Apply cor in addition to saving - njobs = args.njobs[0] # parallel writing - tmin = args.tlim[0] # min time in decimal years - tmax = args.tlim[1] # max time in decimal years - bbox = args.bbox[:] - - print 'parameters:' - for arg in vars(args).iteritems(): - print arg + ifiles = args.files[:] # input files + vnames = args.vnames[:] # lon/lat/h/time variable names + wnames = args.wnames[:] # variables to use for correction + dxy = args.dxy[0] * 1e3 # grid-cell length (km -> m) + radius = args.radius[0] * 1e3 # search radius (km -> m) + nreloc = args.nreloc[0] # number of relocations + proc = args.proc[0] # det, dif, bin or None series + proj = args.proj[0] # EPSG proj number + apply_ = args.apply # Apply cor in addition to saving + njobs = args.njobs[0] # parallel writing + tmin = args.tlim[0] # min time in decimal years + tmax = args.tlim[1] # max time in decimal years + bbox = args.bbox[:] + + print("parameters:") + for arg in vars(args).items(): + print(arg) if njobs == 1: - print 'running sequential code ...' - [main(ifile, vnames, wnames, dxy, proj, radius, nreloc, proc, apply_) \ - for ifile in ifiles] + print("running sequential code ...") + [ + main( + ifile, vnames, wnames, dxy, proj, radius, nreloc, proc, apply_ + ) + for ifile in ifiles + ] else: - print 'running parallel code (%d jobs) ...' % njobs + print("running parallel code (%d jobs) ..." % njobs) from joblib import Parallel, delayed + Parallel(n_jobs=njobs, verbose=5)( - delayed(main)(ifile, vnames, wnames, dxy, proj, radius, nreloc, proc, apply_) \ - for ifile in ifiles) + delayed(main)( + ifile, vnames, wnames, dxy, proj, radius, nreloc, proc, apply_ + ) + for ifile in ifiles + ) - print 'done!' + print("done!") diff --git a/captoolkit/filtst.py b/captoolkit/filtst.py index 62b7159..fa16e61 100755 --- a/captoolkit/filtst.py +++ b/captoolkit/filtst.py @@ -158,7 +158,7 @@ def get_bbox(fname, key="bbox"): fname = fname.split("_") # fname -> list i = fname.index(key) - return map(float, fname[i + 1 : i + 5]) # m + return list(map(float, fname[i + 1 : i + 5])) # m ##NOTE: For stfilter we don't want to pass a subgrid because we @@ -286,7 +286,7 @@ def binning( nb = np.full(N, np.nan) sb = np.full(N, np.nan) - for i in xrange(N): + for i in range(N): t1, t2 = bins[i] (idx,) = np.where((x >= t1) & (x <= t2)) @@ -352,8 +352,7 @@ def get_residuals(tc, hc, order=1, dx=3 / 12.0, window=5 / 12.0): def stfilter( - data, - (xi, yi), + data, xxx_todo_changeme, radius=None, min_obs=25, n_std=3, @@ -361,6 +360,7 @@ def stfilter( window=1 / 12.0, ): + (xi, yi) = xxx_todo_changeme t, x, y, z = data # full file/tile data xi_uniq, yi_uniq = np.unique(xi), np.unique(yi) @@ -377,15 +377,15 @@ def stfilter( np.nanmax(y), ) - print "building KDTree ..." + print("building KDTree ...") Tree = cKDTree(np.column_stack((x, y))) - print "entering spatial loop ..." + print("entering spatial loop ...") - for i_node in xrange(xi.shape[0]): + for i_node in range(xi.shape[0]): if i_node % 500 == 0: - print "node:", i_node + print("node:", i_node) x0, y0 = xi[i_node], yi[i_node] # prediction pt (grid node) @@ -443,8 +443,7 @@ def stfilter( def stfilter2( - data, - (xi, yi), + data, xxx_todo_changeme1, radius=None, min_obs=25, n_std=3, @@ -452,7 +451,7 @@ def stfilter2( window=1 / 12.0, ): """Does the same as above but globally.""" - + (xi, yi) = xxx_todo_changeme1 t, x, y, z = data # full file/tile data # Create output container @@ -486,7 +485,7 @@ def absfilter(t, h, max_abs=50, order=1, step=1 / 12.0, window=1 / 12.0): def main(ifile, args): - print ifile + print(ifile) # ifile = args.ifile[0] bbox = args.bbox[:] @@ -500,7 +499,7 @@ def main(ifile, args): tvar, xvar, yvar, zvar = vnames - print "loading data ..." + print("loading data ...") time, lon, lat, obs = load_data(ifile, tvar, xvar, yvar, zvar, step=1) if len(obs) < MINOBS: @@ -579,18 +578,18 @@ def main(ifile, args): files = args.ifile[:] njobs = args.njobs[0] -print "parameters:" +print("parameters:") -for p in vars(args).iteritems(): - print p +for p in vars(args).items(): + print(p) if njobs == 1: - print "running serial code ..." + print("running serial code ...") [main(f, args) for f in files] else: - print "running parallel code (%d jobs) ..." % njobs + print("running parallel code (%d jobs) ..." % njobs) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=1)(delayed(main)(f, args) for f in files) diff --git a/captoolkit/interpgaus.py b/captoolkit/interpgaus.py index 1b02109..099f4d8 100755 --- a/captoolkit/interpgaus.py +++ b/captoolkit/interpgaus.py @@ -190,7 +190,7 @@ def spatial_filter(x, y, z, dx, dy, sigma=5.0): # Print parameters to screen print('parameters:') -for p in vars(args).items(): print(p) +for p in list(vars(args).items()): print(p) print("reading data ...") diff --git a/captoolkit/interpmed.py b/captoolkit/interpmed.py index 5e80446..ac69fa2 100755 --- a/captoolkit/interpmed.py +++ b/captoolkit/interpmed.py @@ -185,7 +185,7 @@ def spatial_filter(x, y, z, dx, dy, sigma=5.0): # Print parameters to screen print('parameters:') -for p in vars(args).items(): print(p) +for p in list(vars(args).items()): print(p) print("reading data ...") diff --git a/captoolkit/joingrd.py b/captoolkit/joingrd.py index d9a9d4f..0f93f6d 100755 --- a/captoolkit/joingrd.py +++ b/captoolkit/joingrd.py @@ -75,7 +75,7 @@ def get_args(): def print_args(args): print('Input arguments:') - for arg in vars(args).items(): + for arg in list(vars(args).items()): print(arg) @@ -195,7 +195,7 @@ def join(ifiles, suffix=''): # Iterate over tiles for ifile in ifiles: - print('tile:', ifile) + print(('tile:', ifile)) tile_bbox = get_tile_bbox(ifile) i1,i2,j1,j2 = get_tile_position(x_grid, y_grid, tile_bbox) @@ -205,8 +205,8 @@ def join(ifiles, suffix=''): if flipy: flipud(ofile_, vnames) - print('joined tiles:', len(ifiles)) - print('out ->', ofile_) + print(('joined tiles:', len(ifiles))) + print(('out ->', ofile_)) @@ -236,7 +236,7 @@ def join(ifiles, suffix=''): [join(ifiles) for ifiles in allfiles] else: - print('Running parallel code (%d jobs) ...' % njobs) + print(('Running parallel code (%d jobs) ...' % njobs)) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( delayed(join)(ifiles, time_key) for ifiles in allfiles) diff --git a/captoolkit/joingrd.py.bak b/captoolkit/joingrd.py.bak index 1cc592d..d9a9d4f 100755 --- a/captoolkit/joingrd.py.bak +++ b/captoolkit/joingrd.py.bak @@ -74,16 +74,16 @@ def get_args(): def print_args(args): - print 'Input arguments:' - for arg in vars(args).iteritems(): - print arg + print('Input arguments:') + for arg in vars(args).items(): + print(arg) def get_tile_bbox(fname, key='bbox'): """ Extract bbox info from file name. """ fname = fname.split('_') # fname -> list i = fname.index(key) - return map(float, fname[i+1:i+5]) # m + return list(map(float, fname[i+1:i+5])) # m def get_tile_proj(fname, key='epsg'): @@ -123,7 +123,7 @@ def get_grid_shape(tile_shape, num_tiles): def get_grid_names(fname): """ Return all 2d '/variable' names in the HDF5. """ with h5py.File(fname, 'r') as f: - vnames = [k for k in f.keys() if f[k].ndim == 2] + vnames = [k for k in list(f.keys()) if f[k].ndim == 2] return vnames @@ -166,7 +166,7 @@ def group_by_key(files, key='bin'): def flipud(fname, vnames): with h5py.File(fname, 'a') as f: for v in vnames: f[v][:] = np.flipud(f[v][:]) - print 'final grids flipped upside-down' + print('final grids flipped upside-down') def join(ifiles, suffix=''): @@ -175,7 +175,7 @@ def join(ifiles, suffix=''): # Sort input files on keyword number if provided if key: - print 'sorting input files ...' + print('sorting input files ...') natkey = lambda s: int(re.findall(key+'_\d+', s)[0].split('_')[-1]) ifiles.sort(key=natkey) @@ -195,7 +195,7 @@ def join(ifiles, suffix=''): # Iterate over tiles for ifile in ifiles: - print 'tile:', ifile + print('tile:', ifile) tile_bbox = get_tile_bbox(ifile) i1,i2,j1,j2 = get_tile_position(x_grid, y_grid, tile_bbox) @@ -205,8 +205,8 @@ def join(ifiles, suffix=''): if flipy: flipud(ofile_, vnames) - print 'joined tiles:', len(ifiles) - print 'out ->', ofile_ + print('joined tiles:', len(ifiles)) + print('out ->', ofile_) @@ -232,11 +232,11 @@ except: multiple_grids = len(allfiles) > 1 if njobs == 1 or not multiple_grids: - print 'Running sequential code ...' + print('Running sequential code ...') [join(ifiles) for ifiles in allfiles] else: - print 'Running parallel code (%d jobs) ...' % njobs + print('Running parallel code (%d jobs) ...' % njobs) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( delayed(join)(ifiles, time_key) for ifiles in allfiles) diff --git a/captoolkit/merge.py b/captoolkit/merge.py index 06d47b6..064f48e 100755 --- a/captoolkit/merge.py +++ b/captoolkit/merge.py @@ -107,7 +107,7 @@ def merge(ifiles, ofile, vnames, comp): # Iterate over the input files k1 = 0 for ifile in ifiles: - print('reading', ifile) + print(('reading', ifile)) # Write next chunk (the input file) with h5py.File(ifile) as f2: @@ -119,8 +119,8 @@ def merge(ifiles, ofile, vnames, comp): k1 = k2 - print('merged', len(ifiles), 'files') - print('output ->', ofile) + print(('merged', len(ifiles), 'files')) + print(('output ->', ofile)) # Sort input files by key @@ -161,7 +161,7 @@ def sort_files(ifiles, key=None): ifile, ofile = [ifile], [ofile] if njobs > 1 and nfiles > 1: - print('Running parallel code (%d jobs) ...' % njobs) + print(('Running parallel code (%d jobs) ...' % njobs)) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( delayed(merge)(fi, fo, vnames, comp) \ diff --git a/captoolkit/merge.py.bak b/captoolkit/merge.py.bak index 9785cde..06d47b6 100755 --- a/captoolkit/merge.py.bak +++ b/captoolkit/merge.py.bak @@ -64,14 +64,14 @@ def get_total_len(ifiles): N = 0 for fn in ifiles: with h5py.File(fn) as f: - N += f.values()[0].shape[0] + N += list(f.values())[0].shape[0] return N def get_var_names(ifile): """ Return all '/variable' names in the HDF5. """ with h5py.File(ifile, 'r') as f: - vnames = f.keys() + vnames = list(f.keys()) return vnames @@ -81,7 +81,7 @@ def get_multi_io(ifiles, ofile, nfiles): ifiles = [list(arr) for arr in np.array_split(ifiles, nfiles)] # List of output file names fname = os.path.splitext(ofile)[0] + '_%02d.h5' - ofiles = [(fname % k) for k in xrange(len(ifiles))] + ofiles = [(fname % k) for k in range(len(ifiles))] return ifiles, ofiles @@ -95,7 +95,7 @@ def merge(ifiles, ofile, vnames, comp): vnames (list): name of vars to merge. """ # Get length of output containers (from all input files) - print 'Calculating lenght of output from all input files ...' + print('Calculating lenght of output from all input files ...') N = get_total_len(ifiles) with h5py.File(ofile, 'w') as f: @@ -107,11 +107,11 @@ def merge(ifiles, ofile, vnames, comp): # Iterate over the input files k1 = 0 for ifile in ifiles: - print 'reading', ifile + print('reading', ifile) # Write next chunk (the input file) with h5py.File(ifile) as f2: - k2 = k1 + f2.values()[0].shape[0] # first var/first dim + k2 = k1 + list(f2.values())[0].shape[0] # first var/first dim # Iterate over all variables for key in vnames: @@ -119,8 +119,8 @@ def merge(ifiles, ofile, vnames, comp): k1 = k2 - print 'merged', len(ifiles), 'files' - print 'output ->', ofile + print('merged', len(ifiles), 'files') + print('output ->', ofile) # Sort input files by key @@ -128,7 +128,7 @@ def sort_files(ifiles, key=None): """ Sort files by numbers *after* the key in the file name. """ if key: import re - print 'sorting input files ...' + print('sorting input files ...') natkey = lambda s: int(re.findall(key+'_\d+', s)[0].split('_')[-1]) ifiles.sort(key=natkey) @@ -161,12 +161,12 @@ if __name__ == '__main__': ifile, ofile = [ifile], [ofile] if njobs > 1 and nfiles > 1: - print 'Running parallel code (%d jobs) ...' % njobs + print('Running parallel code (%d jobs) ...' % njobs) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( delayed(merge)(fi, fo, vnames, comp) \ for fi,fo in zip(ifile, ofile)) else: - print 'Running sequential code ...' + print('Running sequential code ...') [merge(fi, fo, vnames, comp) for fi,fo in zip(ifile, ofile)] diff --git a/captoolkit/tile.py b/captoolkit/tile.py index 392b0bb..00e9f23 100755 --- a/captoolkit/tile.py +++ b/captoolkit/tile.py @@ -77,7 +77,7 @@ def get_args(): def print_args(args): print('Input arguments:') - for arg in vars(args).items(): + for arg in list(vars(args).items()): print(arg) @@ -181,7 +181,7 @@ def get_tile_data(ifile, x, y, bbox, buff=1, proj='3031', tile_num=0): fo.flush() - if npts != 0: print('tile %03d: #points' % tile_num, npts, '...') + if npts != 0: print(('tile %03d: #points' % tile_num, npts, '...')) try: fo.close() @@ -225,9 +225,9 @@ def count_files(ifiles, key='*tile*'): fxys = [(f,x,y) for f,(x,y) in zip(ifiles, xys)] # [(f1,x1,y1), (f2,x2,y2)..] fxybs = [(f,x,y,b,n+1) for (f,x,y) in fxys for n,b in enumerate(bboxs)] # [(f1,x1,y1,b1,1), (f1,x1,y1,b2,2)..] -print('number of files:', len(ifiles)) -print('number of tiles:', len(bboxs)) -print('number of tasks:', len(fxybs)) +print(('number of files:', len(ifiles))) +print(('number of tiles:', len(bboxs))) +print(('number of tasks:', len(fxybs))) # NOTE: For each bbox scan full file @@ -236,9 +236,9 @@ def count_files(ifiles, key='*tile*'): [get_tile_data(f, x, y, b, dr, proj, n) for f,x,y,b,n in fxybs] else: - print('running parallel code (%d jobs) ...' % njobs) + print(('running parallel code (%d jobs) ...' % njobs)) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( delayed(get_tile_data)(f, x, y, b, dr, proj, n) for f,x,y,b,n in fxybs) -print('number of tiles with data:', count_files(ifiles)) +print(('number of tiles with data:', count_files(ifiles))) diff --git a/captoolkit/tile.py.bak b/captoolkit/tile.py.bak index 23db633..392b0bb 100755 --- a/captoolkit/tile.py.bak +++ b/captoolkit/tile.py.bak @@ -14,7 +14,16 @@ complete the tilling. Notes: Bedmap boundaries: -b -3333000 3333000 -3333000 3333000 Ross boundaries: -b -600000 410000 -1400000 -400000 - ITS_LIVE boundaries: -b -2678407.5 2816632.5 -2154232.5 2259847.5 + ITS_LIVE boundaries: -b -2678407.5 2816632.5 -2154232.5 2259847.5 + +Credits: + captoolkit - JPL Cryosphere Altimetry Processing Toolkit + + Fernando Paolo (paolofer@jpl.nasa.gov) + Johan Nilsson (johan.nilsson@jpl.nasa.gov) + Alex Gardner (alex.s.gardner@jpl.nasa.gov) + + Jet Propulsion Laboratory, California Institute of Technology """ import os @@ -67,9 +76,9 @@ def get_args(): def print_args(args): - print 'Input arguments:' - for arg in vars(args).iteritems(): - print arg + print('Input arguments:') + for arg in vars(args).items(): + print(arg) def transform_coord(proj1, proj2, x, y): @@ -172,7 +181,7 @@ def get_tile_data(ifile, x, y, bbox, buff=1, proj='3031', tile_num=0): fo.flush() - if npts != 0: print 'tile %03d: #points' % tile_num, npts, '...' + if npts != 0: print('tile %03d: #points' % tile_num, npts, '...') try: fo.close() @@ -205,7 +214,7 @@ print_args(args) if len(ifiles) == 1: ifiles = glob(ifiles[0]) # pass str if "Argument list too long" # Generate list of items (only once!) for parallel proc -print 'generating list of tasks (files x tiles) ...' +print('generating list of tasks (files x tiles) ...') # NOTE: To avoid reprojecting coords (lon/lat -> x/y) several # times for each file (one conversion per tile), coordinates are @@ -216,20 +225,20 @@ xys = [get_xy(f, vnames, proj) for f in ifiles] # [(x1, fxys = [(f,x,y) for f,(x,y) in zip(ifiles, xys)] # [(f1,x1,y1), (f2,x2,y2)..] fxybs = [(f,x,y,b,n+1) for (f,x,y) in fxys for n,b in enumerate(bboxs)] # [(f1,x1,y1,b1,1), (f1,x1,y1,b2,2)..] -print 'number of files:', len(ifiles) -print 'number of tiles:', len(bboxs) -print 'number of tasks:', len(fxybs) +print('number of files:', len(ifiles)) +print('number of tiles:', len(bboxs)) +print('number of tasks:', len(fxybs)) # NOTE: For each bbox scan full file if njobs == 1: - print 'running sequential code ...' + print('running sequential code ...') [get_tile_data(f, x, y, b, dr, proj, n) for f,x,y,b,n in fxybs] else: - print 'running parallel code (%d jobs) ...' % njobs + print('running parallel code (%d jobs) ...' % njobs) from joblib import Parallel, delayed Parallel(n_jobs=njobs, verbose=5)( delayed(get_tile_data)(f, x, y, b, dr, proj, n) for f,x,y,b,n in fxybs) -print 'number of tiles with data:', count_files(ifiles) +print('number of tiles with data:', count_files(ifiles)) diff --git a/captoolkit/utils.py b/captoolkit/utils.py index 9724da3..bf5116a 100755 --- a/captoolkit/utils.py +++ b/captoolkit/utils.py @@ -16,7 +16,7 @@ def print_args(args): """Print arguments passed to argparse.""" print("Input arguments:") - for arg in vars(args).items(): + for arg in list(vars(args).items()): print(arg) diff --git a/captoolkit/utils.py.bak b/captoolkit/utils.py.bak index c3ccbcb..9724da3 100755 --- a/captoolkit/utils.py.bak +++ b/captoolkit/utils.py.bak @@ -16,7 +16,7 @@ from scipy import signal def print_args(args): """Print arguments passed to argparse.""" print("Input arguments:") - for arg in vars(args).iteritems(): + for arg in vars(args).items(): print(arg) @@ -37,7 +37,7 @@ def save_h5(fname, vardict, mode="a"): vardict : {'name1': var1, 'name2': va2, 'name3': var3} """ with h5py.File(fname, mode) as f: - for k, v in vardict.items(): + for k, v in list(vardict.items()): if k in f: f[k][:] = np.squeeze(v) else: @@ -48,7 +48,7 @@ def is_empty(ifile): """Test if file is corruted or empty""" try: with h5py.File(ifile, "r") as f: - if bool(f.keys()): + if bool(list(f.keys())): return False else: return True