diff --git a/src/pyclaw/gauges.py b/src/pyclaw/gauges.py index 27983d62..f908ea88 100644 --- a/src/pyclaw/gauges.py +++ b/src/pyclaw/gauges.py @@ -2,9 +2,13 @@ """Module contains class definitions related to dealing with gauge data""" +# TODO: switch to pathlib for all file handling +from pathlib import Path import os -import sys -import numpy +import argparse +import warnings + +import numpy as np # See if pandas is available try: @@ -75,7 +79,7 @@ def __init__(self, gauge_id=None, path=None, use_pandas=False): def read(self, gauge_id, path=None, use_pandas=False): r"""Read the gauge file at path into this object - Read in the gauge with gauge id `gauge_id` located at `path`. + Read in the gauge with gauge id `gauge_id` located at `path`. If `use_pandas` is `True` then `q` will be a Pandas frame. (not yet implemented) @@ -137,10 +141,10 @@ def read(self, gauge_id, path=None, use_pandas=False): #data = line.split() #nvals = 1 + len(data)-6 # should agree with num_eqn - + # Check to see if there is also data in the .txt file, # otherwise perhaps it's in a binary .bin file. - + line = gauge_file.readline() if len(line) > 0: if 'binary32' in line: @@ -172,7 +176,7 @@ def read(self, gauge_id, path=None, use_pandas=False): if file_format == 'ascii': # data follows header in .txt file: - data = numpy.loadtxt(gauge_path, comments="#") + data = np.loadtxt(gauge_path, comments="#") if data.ndim == 1: # only one line in gauge file, expand to 2d array data = data.reshape((1,len(data))) @@ -184,28 +188,27 @@ def read(self, gauge_id, path=None, use_pandas=False): if not os.path.isfile(gauge_path): msg = 'No data in .txt file and did not find ' \ + '\n binary file %s' % gauge_path - import warnings warnings.warn(msg) return if file_format in ['binary','binary64']: - data = numpy.fromfile(gauge_path, dtype=numpy.float64) + data = np.fromfile(gauge_path, dtype=np.float64) elif file_format == 'binary32': - data = numpy.fromfile(gauge_path, dtype=numpy.float32) + data = np.fromfile(gauge_path, dtype=np.float32) # assume rows are: level, t, q[0:num_eqn] nrows = 2 + num_eqn - assert numpy.mod(len(data),nrows) == 0, \ + assert np.mod(len(data), nrows) == 0, \ '*** unexpected number of values in gauge file' \ + '\n*** expected nrows = %i rows' % nrows ncols = int(len(data)/nrows) data = data.reshape((nrows,ncols), order='F').T - self.level = data[:, 0].astype(numpy.int64) + self.level = data[:, 0].astype(np.int64) self.t = data[:, 1] self.q = data[:, 2:].transpose() - + if num_eqn != self.q.shape[0]: raise ValueError("Number of fields in gauge file does not match", "recorded number in header.") @@ -215,11 +218,11 @@ def read(self, gauge_id, path=None, use_pandas=False): if self.gtype == 'lagrangian': # for lagrangian gauges x,y paths are stored in q in place of u,v # note: only implemented in 2d geoclaw for now! - self.particle_path = numpy.vstack((self.t, self.q[1,:], + self.particle_path = np.vstack((self.t, self.q[1,:], self.q[2,:])).T else: # set the lagrangian path to a single fixed location at t=t0: - self.particle_path = numpy.array([[self.t[0], self.location[0], + self.particle_path = np.array([[self.t[0], self.location[0], self.location[1]]]) @@ -333,11 +336,15 @@ def compare_gauges(paths, gauge_id, fields='all'): axes.set_ylabel("q[%s, :]" % n) axes.legend() - axes = fig.add_subplot(len(fields), 2, 2 * i + 2) - axes.plot(gauges[0].t, - numpy.abs(gauges[0].q[n, :] - gauges[1].q[n, :]), 'r') - axes.set_xlabel("t") - axes.set_ylabel("$|q_{old}[%s, :] - q_{new}[%s, :]|$" % (n, n)) + if gauges[0].t.shape == gauges[1].t.shape: + axes = fig.add_subplot(len(fields), 2, 2 * i + 2) + axes.plot(gauges[0].t, + np.abs(gauges[0].q[n, :] - gauges[1].q[n, :]), 'r') + axes.set_xlabel("t") + axes.set_ylabel("$|q_{old}[%s, :] - q_{new}[%s, :]|$" % (n, n)) + else: + warnings.warn("Gauge time series do not match, skipping direct " + + "comparison.") return fig @@ -363,17 +370,17 @@ def convert_gauges(path, output_path=None): output_path = os.getcwd() # Load old data - data = numpy.loadtxt(path) - old_ids = numpy.asarray(data[:, 0], dtype=int) - unique_ids = numpy.asarray(list(set(old_ids))) + data = np.loadtxt(path) + old_ids = np.asarray(data[:, 0], dtype=int) + unique_ids = np.asarray(list(set(old_ids))) # Create new gauges and compare for gauge_id in unique_ids: - gauge_indices = numpy.nonzero(old_ids == gauge_id)[0] + gauge_indices = np.nonzero(old_ids == gauge_id)[0] new_gauge = GaugeSolution() new_gauge.id = gauge_id - new_gauge.location = (numpy.nan, numpy.nan) - new_gauge.level = numpy.asarray(data[gauge_indices, 1], dtype=int) + new_gauge.location = (np.nan, np.nan) + new_gauge.level = np.asarray(data[gauge_indices, 1], dtype=int) new_gauge.t = data[gauge_indices, 2] new_gauge.q = data[gauge_indices, 3:].transpose() new_gauge.write(output_path) @@ -386,7 +393,7 @@ def convert_gauges(path, output_path=None): def compare_old_gauges(old_path, new_path, gauge_id, plot=False, abs_tol=1e-14, - rel_tol=0.0, + rel_tol=0.0, verbose=False): r"""Compare old gauge data at `path` to new gauge data at same path @@ -413,9 +420,9 @@ def compare_old_gauges(old_path, new_path, gauge_id, plot=False, abs_tol=1e-14, """ # Load old gauge data - data = numpy.loadtxt(old_path) - old_ids = numpy.asarray(data[:, 0], dtype=int) - gauge_indices = numpy.nonzero(old_ids == gauge_id)[0] + data = np.loadtxt(old_path) + old_ids = np.asarray(data[:, 0], dtype=int) + gauge_indices = np.nonzero(old_ids == gauge_id)[0] q = data[gauge_indices, 3:] # Load new data @@ -425,9 +432,9 @@ def compare_old_gauges(old_path, new_path, gauge_id, plot=False, abs_tol=1e-14, if verbose: print("Comparison of gauge %s:" % gauge_id) print(r" ||\Delta q||_2 = ", - numpy.linalg.norm(q - gauge.q.transpose(), ord=2)) + np.linalg.norm(q - gauge.q.transpose(), ord=2)) print(r" arg(||\Delta q||_\infty = ", - numpy.argmax(q - gauge.q.transpose())) + np.argmax(q - gauge.q.transpose())) if plot: import matplotlib.pyplot as plt @@ -439,7 +446,7 @@ def compare_old_gauges(old_path, new_path, gauge_id, plot=False, abs_tol=1e-14, axes.set_xlabel("t (s)") axes.set_ylabel("q(%s, :)" % i) - return numpy.allclose(q, gauge.q.transpose(), rtol=rel_tol, atol=abs_tol) + return np.allclose(q, gauge.q.transpose(), rtol=rel_tol, atol=abs_tol) def check_old_gauge_data(path, gauge_id, new_gauge_path="./regression_data"): @@ -451,27 +458,27 @@ def check_old_gauge_data(path, gauge_id, new_gauge_path="./regression_data"): :Input: - *path* (string) - Path to old gauge data file - *gauge_id* (int) - Gauge id to compare - - *new_gauge_path* (path) - Path to directory containing new gauge files, + - *new_gauge_path* (path) - Path to directory containing new gauge files, defaults to './regression_data'. :Output: - - (figure) Returns a matplotlib figure object plotting the differences in + - (figure) Returns a matplotlib figure object plotting the differences in time. """ import matplotlib.pyplot as plt # Load old gauge data - data = numpy.loadtxt(path) - old_ids = numpy.asarray(data[:, 0], dtype=int) - gauge_indices = numpy.nonzero(old_ids == gauge_id)[0] + data = np.loadtxt(path) + old_ids = np.asarray(data[:, 0], dtype=int) + gauge_indices = np.nonzero(old_ids == gauge_id)[0] q = data[gauge_indices, 3:] # Load new data gauge = GaugeSolution(gauge_id, new_gauge_path) - print(numpy.linalg.norm(q - gauge.q.transpose(), ord=2)) - print(numpy.argmax(q - gauge.q.transpose())) + print(np.linalg.norm(q - gauge.q.transpose(), ord=2)) + print(np.argmax(q - gauge.q.transpose())) fig = plt.figure() for i in range(gauge.q.shape[0]): @@ -481,34 +488,29 @@ def check_old_gauge_data(path, gauge_id, new_gauge_path="./regression_data"): return fig - if __name__ == "__main__": - import matplotlib.pyplot as plt + parser = argparse.ArgumentParser() + parser.add_argument('path', type=Path, nargs=2, + help="Paths to gauge data directories") + parser.add_argument('--id', '-n', type=int, nargs='+', default=None, + help="Gauge IDs to plot, " + + "defaults to all found in first path") + parser.add_argument('--fields', '-f', type=int, nargs='+', default=None, + help="Fields to plot, defaults to all") + args = parser.parse_args() + + if args.id: + gauge_ids = args.id + else: + gauge_ids = [int(str(path.name[5:10])) + for path in args.path[0].glob("gauge*.txt")] + if args.fields: + fields = args.fields + else: + fields = 'all' + + for gauge_id in gauge_ids: + fig = compare_gauges(args.path, gauge_id, fields) - help_msg = \ -"""gauges.py path1 path2 [gauge_id] [fields...] - -Plots a comparison between the gauges at path1 and path2 with gauge_id and the -fields specified. Only one gauge_id can be specified at a time but a number of -fields can be specified including 'all'. -""" - - fields = 'all' - gauge_id = 1 - if len(sys.argv) < 3: - print(help_msg) - sys.exit(0) - elif len(sys.argv) >= 3: - paths = [str(sys.argv[1]), str(sys.argv[2])] - if len(sys.argv) > 3: - gauge_id = int(sys.argv[3]) - if len(sys.argv) > 4: - if sys.argv[4].lower() == 'all': - fields = 'all' - else: - fields = [int(field for field in sys.argv[4:])] - - - fig = compare_gauges(paths, gauge_id, fields) plt.show()