Skip to content

Commit

Permalink
Update todo, convert py2 -> py3
Browse files Browse the repository at this point in the history
  • Loading branch information
fspaolo committed Feb 12, 2020
1 parent 20cb0e0 commit 2931e54
Show file tree
Hide file tree
Showing 16 changed files with 930 additions and 646 deletions.
2 changes: 1 addition & 1 deletion TODO.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 2 additions & 2 deletions captoolkit/corrlaser.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
118 changes: 73 additions & 45 deletions captoolkit/corrlaser.py.bak
Original file line number Diff line number Diff line change
@@ -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 ([email protected])
Johan Nilsson ([email protected])
Alex Gardner ([email protected])
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.")
38 changes: 19 additions & 19 deletions captoolkit/corrscatt.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand Down Expand Up @@ -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()

Expand All @@ -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():
Expand Down Expand Up @@ -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]

Expand Down Expand Up @@ -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:
Expand All @@ -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)(
Expand Down
Loading

0 comments on commit 2931e54

Please sign in to comment.