From e0632188e6a2de2d17d24cd9ae22e731d62a1ba2 Mon Sep 17 00:00:00 2001 From: alessandrofelder Date: Fri, 21 May 2021 17:32:20 +0100 Subject: [PATCH 1/6] small fixes remove unused import. ensure diff_threshold arg is parsed to float. --- pycascadia/remove_restore.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pycascadia/remove_restore.py b/pycascadia/remove_restore.py index b625d03..5f0ded3 100644 --- a/pycascadia/remove_restore.py +++ b/pycascadia/remove_restore.py @@ -5,7 +5,7 @@ It's just a way for us to learn how to use pyGMT and how to read/write the sample data. """ -from pygmt import blockmedian, surface, grdtrack, grdcut, grdfilter +from pygmt import blockmedian, surface, grdtrack, grdfilter import os import pandas as pd import matplotlib.pyplot as plt @@ -63,7 +63,7 @@ def main(): parser.add_argument('filenames', nargs='+', help='sources to combine with the base grid') parser.add_argument('--base', required=True, help='base grid') parser.add_argument('--spacing', type=float, help='output grid spacing') - parser.add_argument('--diff_threshold', default=0.0, help='value above which differences will be added to the base grid') + parser.add_argument('--diff_threshold', type=float, default=0.0, help='value above which differences will be added to the base grid') parser.add_argument('--plot', action='store_true', help='plot final output before saving') parser.add_argument('--output', required=True, help='filename of final output') parser.add_argument('--region_of_interest', required=False, nargs=4, type=float, From f3d6aeafedb644cec7691dad0fd4f624dd47c0c6 Mon Sep 17 00:00:00 2001 From: alessandrofelder Date: Fri, 21 May 2021 17:35:40 +0100 Subject: [PATCH 2/6] define function for grid preprocessing --- pycascadia/remove_restore.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/pycascadia/remove_restore.py b/pycascadia/remove_restore.py index 5f0ded3..4f98bcc 100644 --- a/pycascadia/remove_restore.py +++ b/pycascadia/remove_restore.py @@ -56,6 +56,24 @@ def load_base_grid(fname, region=None, spacing=None): return base_grid +def preprocess_base_grid(base_grid, update_grid, final_spacing): + xyz1 = base_grid.as_xyz() + xyz2 = update_grid.as_xyz() + xyz1.append(xyz2, ignore_index=True) + + max_spacing = max(update_grid.spacing, final_spacing) + minimal_region = min_regions(update_grid.region, base_grid.region) + bmd = blockmedian(xyz1, spacing=4*max_spacing, region=minimal_region, Q=True) + combined = surface(bmd.x, bmd.y, bmd.z, spacing=4*max_spacing, region=minimal_region) + + filter_combined = grdfilter(combined, D=0, F=f"c{12*max_spacing}") + fname = "filtered.nc" + filter_combined.to_netcdf(fname) + combined_grid = Grid(fname) + os.remove(fname) + combined_grid.resample(max_spacing) + + return combined_grid def main(): # Handle arguments From 7c4a698b781292993a8522a2ac553ea8d9c606d9 Mon Sep 17 00:00:00 2001 From: alessandrofelder Date: Fri, 21 May 2021 17:38:22 +0100 Subject: [PATCH 3/6] document preprocess_base_grid --- pycascadia/remove_restore.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pycascadia/remove_restore.py b/pycascadia/remove_restore.py index 4f98bcc..bed53c5 100644 --- a/pycascadia/remove_restore.py +++ b/pycascadia/remove_restore.py @@ -57,6 +57,11 @@ def load_base_grid(fname, region=None, spacing=None): return base_grid def preprocess_base_grid(base_grid, update_grid, final_spacing): + """ + Combines and smooths data from the base and update grid, + as described in workflow steps B+C of the GEBCO cookbook. + This is a preprocessing step to remove-restore. + """ xyz1 = base_grid.as_xyz() xyz2 = update_grid.as_xyz() xyz1.append(xyz2, ignore_index=True) From 51063cee632ce109d9fd2f21e78c4eb4767a7383 Mon Sep 17 00:00:00 2001 From: alessandrofelder Date: Fri, 21 May 2021 17:39:14 +0100 Subject: [PATCH 4/6] don't resample base_grid on loading will be resampled as part of preprocessing step instead. --- pycascadia/remove_restore.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pycascadia/remove_restore.py b/pycascadia/remove_restore.py index bed53c5..23b4112 100644 --- a/pycascadia/remove_restore.py +++ b/pycascadia/remove_restore.py @@ -51,9 +51,6 @@ def load_base_grid(fname, region=None, spacing=None): base_grid = Grid(fname, convert_to_xyz=False) if region: base_grid.crop(region) - if spacing: - base_grid.resample(spacing) - return base_grid def preprocess_base_grid(base_grid, update_grid, final_spacing): From 353069ecbddb8f7294ca7792677157e443893ad3 Mon Sep 17 00:00:00 2001 From: alessandrofelder Date: Fri, 21 May 2021 17:40:19 +0100 Subject: [PATCH 5/6] filter any null difference values this is needed because low-pass filter in preprocessing may not have values everywhere. --- pycascadia/remove_restore.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pycascadia/remove_restore.py b/pycascadia/remove_restore.py index 23b4112..746aadd 100644 --- a/pycascadia/remove_restore.py +++ b/pycascadia/remove_restore.py @@ -25,13 +25,14 @@ def calc_diff_grid(base_grid, update_grid, diff_threshold=0.0): print("Find z in base grid") base_pts = grdtrack(bmd, base_grid.grid, 'base_z', interpolation='l') - print ("Create difference grid") + print("Create difference grid") diff = pd.DataFrame() diff['x'] = base_pts['x'] diff['y'] = base_pts['y'] diff['z'] = base_pts['z'] - base_pts['base_z'] diff[diff.z.abs() < diff_threshold]['z'] = 0.0 # Filter out small differences + diff[diff.z.isnull()] = 0.0 # Filter out NaNs. diff_xyz_fname = "diff.xyz" diff_grid_fname = "diff.nc" From 75187734428dfdb0457693751bf04018f9abe0b0 Mon Sep 17 00:00:00 2001 From: alessandrofelder Date: Fri, 21 May 2021 17:41:20 +0100 Subject: [PATCH 6/6] incorporate preprocessing functionality into main() --- pycascadia/remove_restore.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/pycascadia/remove_restore.py b/pycascadia/remove_restore.py index 746aadd..72bd0a4 100644 --- a/pycascadia/remove_restore.py +++ b/pycascadia/remove_restore.py @@ -105,23 +105,27 @@ def main(): print("Loading update grid") update_grid = Grid(fname, convert_to_xyz=True) - diff_grid = calc_diff_grid(base_grid, update_grid, diff_threshold=diff_threshold) + print("Combining grids") + combined_base_grid = preprocess_base_grid(base_grid, update_grid, args.spacing) + + diff_grid = calc_diff_grid(combined_base_grid, update_grid, diff_threshold=diff_threshold) print("Update base grid") - base_grid.grid.values += diff_grid.values + combined_base_grid.grid.values += diff_grid.values - base_grid.save_grid(args.output) + # TODO How should this work with several input?! + combined_base_grid.save_grid(args.output) if args.plot: fig, axes = plt.subplots(2,2) initial_base_grid = load_base_grid(base_fname, region=args.region_of_interest) initial_base_grid.plot(ax=axes[0,0]) axes[0,0].set_title("Initial Grid") - base_grid.plot(ax=axes[0,1]) + combined_base_grid.plot(ax=axes[0,1]) axes[0,1].set_title("Final Grid") - base_grid.grid.differentiate('x').plot(ax=axes[1,0]) + combined_base_grid.grid.differentiate('x').plot(ax=axes[1,0]) axes[1,0].set_title("x Derivative of Final Grid") - base_grid.grid.differentiate('y').plot(ax=axes[1,1]) + combined_base_grid.grid.differentiate('y').plot(ax=axes[1,1]) axes[1,1].set_title("y Derivative of Final Grid") plt.show()