Merge pull request #520 from matthewhoffman/landice/hydro_stats_plot

matthewhoffman · web-flow · commit 598d26766800 · 2025-11-10T11:48:55.000-07:00
Add script to plot hydro mass balance terms

This PR adds new scripts related to MALI subglacial hydrology:
* script with plots for evaluating subglacial hydrology mass balance from global stats
* script for preprocessing ice geometry input file to remove surface depressions. This greatly helps with lake-related instabilities
diff --git a/landice/mesh_tools_li/fill_supraglacial_lakes.py b/landice/mesh_tools_li/fill_supraglacial_lakes.py
@@ -0,0 +1,272 @@
+#!/usr/bin/env python3
+"""
+fill_supraglacial_lakes.py
+
+fill depressions in the grounded ice-sheet upper surface
+"""
+
+import argparse
+import heapq
+import mosaic
+import numpy as np
+import sys
+import xarray as xr
+
+from matplotlib import pyplot as plt
+import matplotlib.cm as cm
+import matplotlib.colors as mcolors
+from compass.landice.mesh import mpas_flood_fill
+
+rhoi = 910.0
+rhoo = 1028.0
+
+
+class CyclicNormalize(mcolors.Normalize):
+    def __init__(self, cmin=0, cmax=1, vmin=0, vmax=1, clip=False):
+        self.cmin = cmin
+        self.cmax = cmax
+        mcolors.Normalize.__init__(self, vmin, vmax, clip=clip)
+
+    def __call__(self, value, clip=False):
+        x, y = [self.cmin, self.cmax], [0, 1]
+        return np.ma.masked_array(np.interp(value, x, y, period=self.cmax - self.cmin))
+
+
+def parse_arguments():
+    """
+    Parse command-line arguments.
+
+    Returns
+    -------
+    argparse.Namespace
+        Parsed arguments containing input file paths, parameters, and options.
+    """
+    parser = argparse.ArgumentParser(
+        description="Run a priority-flood fill on a 2D unstructured potential surface."
+    )
+
+    parser.add_argument(
+        "-i", "--input",
+        type=str,
+        required=True,
+        help="Path to input file"
+    )
+
+    parser.add_argument(
+        "-o", "--output",
+        type=str,
+        required=True,
+        help="Path to output file"
+    )
+
+    parser.add_argument(
+        "--verbose",
+        action="store_true",
+        help="Enable verbose output for debugging."
+    )
+
+    return parser.parse_args()
+
+def run(args):
+    """
+    Load MALI geometry data and calculates needed derived fields
+
+    Returns
+    -------
+    xarray.dataset
+        
+        
+    """
+    ds = xr.open_dataset(args.input, decode_times=False, decode_cf=False)
+    thk = ds.thickness[0, :]
+    bed = ds.bedTopography[0, :]
+    ds["grd_mask"] = ((thk * rhoi / rhoo +bed) > 0.0) * (thk > 0.0)
+    ds["upperSurfaceGrd"] = ds.grd_mask * (bed + thk)
+    usrf = ds['upperSurfaceGrd'].values[:]
+    nCells = ds.nCells.values
+    cellsOnCell = ds.cellsOnCell.values
+    nEdgesOnCell = ds.nEdgesOnCell.values
+
+    # Create mask for grounded margin, defined here to include transition
+    # between floating and grounded
+    ds["grd_margin_mask"] = xr.zeros_like(ds.upperSurfaceGrd, dtype=int)
+    for iCell in np.argwhere(ds.grd_mask.values == 1):
+        neighbor_indices = cellsOnCell[iCell, :nEdgesOnCell[iCell][0]] - 1
+        neighbor_indices = neighbor_indices[neighbor_indices >= 0]
+        if ds.grd_mask.values[neighbor_indices].sum() != len(neighbor_indices):
+            ds.grd_margin_mask[iCell] = 1
+
+    # create mask for hydrologically connected region to margin
+    #print("Creating mask of hydrologically connected region to margin")
+    #ds['mgn_draining_mask'] = xr.zeros_like(ds.upperSurfaceGrd, dtype=int)
+    #iter = 0
+    #seed_mask = ds['grd_margin_mask'].values[:]
+    #print(f"Size of grd_margin_mask={len(seed_mask)}")
+    #grow_mask = ds["grd_mask"].values[:]
+    #new_keep_mask = seed_mask.copy()
+    #keep_mask = new_keep_mask * 0
+    #n_mask_cells = keep_mask.sum()
+    #grow_iters=sys.maxsize
+    #new_ind = np.nonzero((new_keep_mask - keep_mask) == 1)[0]
+    #for iter in range(grow_iters):
+    #    new_keep_mask = keep_mask.copy()
+    #    print(f'iter={iter}, keep_mask size={keep_mask.sum()}')
+
+    #    for iCell in new_ind:
+    #        neighs = cellsOnCell[iCell, :nEdgesOnCell[iCell]] - 1
+    #        neighs = neighs[neighs >= 0]  # drop garbage cell
+    #        for jCell in neighs:
+    #            if grow_mask[jCell] == 1 and usrf[jCell] >= usrf[iCell]:
+    #                new_keep_mask[jCell] = 1
+
+    #    # only search over new cells added previous iteration
+    #    new_ind = np.nonzero((new_keep_mask - keep_mask) == 1)[0]
+    #    keep_mask = new_keep_mask.copy()
+
+    #    n_mask_cells_new = keep_mask.sum()
+    #    if n_mask_cells_new == n_mask_cells:
+    #        break
+    #    n_mask_cells = n_mask_cells_new
+    #    iter += 1
+    #ds["mgn_draining_mask"][:] = keep_mask
+
+    # ---------------
+    # find depressions
+    grd_mask = ds['grd_mask'].values[:]
+    margin = ds['grd_margin_mask'].values[:]
+    filled = usrf.copy()
+    visited = np.zeros_like(usrf, dtype=bool)
+    pq = []
+
+    # add margin cells to heap to start
+    for iCell in np.where(margin == 1)[0]:
+        heapq.heappush(pq, (usrf[iCell], iCell))
+        visited[iCell] = True
+
+    while pq:
+        ht, iCell = heapq.heappop(pq)
+        neighbor_indices = cellsOnCell[iCell, 0:nEdgesOnCell[iCell]] - 1
+        neighbor_indices = neighbor_indices[neighbor_indices >= 0]
+        for jCell in neighbor_indices:
+            if not visited[jCell] and grd_mask[jCell]:
+                visited[jCell] = True
+                filled[jCell] = max(usrf[jCell], ht)
+                heapq.heappush(pq, (filled[jCell], jCell))
+    ds['usrf_filled'] = xr.zeros_like(ds.upperSurfaceGrd)
+    ds['usrf_filled'][:] = filled
+
+    diff = ds['usrf_filled'] - ds['upperSurfaceGrd']
+    diff = diff.values
+
+    print("Filling statistics:")
+    print(f"#cells adjusted by >  0m: {(diff>0.).sum()}")
+    print(f"#cells adjusted by >  5m: {(diff>5.).sum()}")
+    print(f"#cells adjusted by > 10m: {(diff>10.).sum()}")
+    print(f"#cells adjusted by > 20m: {(diff>20.).sum()}")
+    print(f"#cells adjusted by > 50m: {(diff>50.).sum()}")
+    print(f"#cells adjusted by >100m: {(diff>100.).sum()}")
+
+
+    dsout = xr.open_dataset(args.input, decode_times=False, decode_cf=False)
+    dsout['thickness'] = dsout['thickness'] + (ds['usrf_filled'] - ds['upperSurfaceGrd'])
+    dsout.to_netcdf(args.output)
+
+
+    if args.verbose:
+        # plot result
+        fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(14, 10),
+                                 constrained_layout=True,
+                                 sharex=True, sharey=True)
+        axes = axes.flatten()
+        descriptor = mosaic.Descriptor(ds)
+
+        ax = 0
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               c=(thk > 0.0), aa=False)
+        axes[ax].set_title(f"ice mask")
+
+        ax += 1
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               ds.grd_mask, aa=False)
+        axes[ax].set_title(f"grd mask")
+
+        ax += 1
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               ds.grd_margin_mask, aa=False)
+        axes[ax].set_title(f"grounded margin mask")
+
+        #pc = mosaic.polypcolor(axes[3], descriptor,
+        #                       ds.mgn_draining_mask + ds.grd_mask, aa=False)
+        #axes[3].set_title(f"margin draining mask + grd mask")
+
+        ax += 1
+        period = 200.0
+        norm = mcolors.Normalize(vmin=0, vmax=period)
+        cyclicnorm = CyclicNormalize(cmin=0., cmax=400.0, vmin=0.0, vmax=4500.0)
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               ds['upperSurfaceGrd'], aa=False,
+                               cmap="hsv", norm=cyclicnorm)
+        #pc = mosaic.polypcolor(axes[ax], descriptor,
+        #                       ds.mgn_draining_mask * np.nan,
+        #                       ec=np.where(ds.mgn_draining_mask.values == 1, "white", "black"))
+        ##fig.colorbar(pc, ax=axes[ax], fraction=0.1,
+        ##             label=f"surface elevation (m)")
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               ds['usrf_filled'] * np.nan,
+                               ec=np.where(usrf == filled, "white", "black"),
+                               lw=np.where(usrf == filled, 0.01, 0.5))
+        axes[ax].set_title(f"upper surface")
+
+        ax += 1
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               ds['usrf_filled'], aa=False,
+                               cmap="hsv", norm=cyclicnorm)
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               ds['usrf_filled'] * np.nan,
+                               ec=np.where(usrf == filled, "white", "black"),
+                               lw=np.where(usrf == filled, 0.01, 0.5))
+        axes[ax].set_title(f"upper surface filled")
+
+        ax += 1
+        diff = ds['usrf_filled'] - ds['upperSurfaceGrd']
+        diff = diff.where(diff != 0.0)
+        cmap = cm.get_cmap("rainbow", 10).copy()
+        cmap.set_bad(color="white")
+        pc = mosaic.polypcolor(axes[ax], descriptor,
+                               diff,
+                               vmin=0, vmax=50,
+                               aa=False, cmap=cmap)
+        axes[ax].set_title(f"upper surface adjustment")
+        fig.colorbar(pc, ax=axes[ax], fraction=0.1)
+
+
+        for ax in axes:
+            ax.set_aspect('equal')
+            ax.set_xlabel("X (m)")
+            ax.set_ylabel("Y (m)")
+
+
+
+
+def main():
+    """
+    Main entry point for the script.
+    """
+    args = parse_arguments()
+
+    if args.verbose:
+        print("Arguments parsed successfully:")
+        print(args)
+
+    # perform processing
+    run(args)
+
+    if args.verbose:
+        print("Execution complete.")
+        plt.show()
+
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/landice/output_processing_li/plot_subglacial_water_mass_balance.py b/landice/output_processing_li/plot_subglacial_water_mass_balance.py
@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+'''
+Script to plot subglacial hydrology time-series from a landice globalStats file.
+'''
+
+from __future__ import absolute_import, division, print_function, unicode_literals
+
+import netCDF4 as nc
+import numpy as np
+from matplotlib import pyplot as plt
+import argparse
+
+rhow = 1000.0
+secyr = 3600.0 * 24.0 * 365.0
+
+parser = argparse.ArgumentParser(description=__doc__)
+parser.add_argument("-f", dest="filename", help="input filename", default="globalStats.nc", metavar="FILENAME")
+parser.add_argument("-u", dest="units", help="units for mass: kg, Gt", default="Gt", metavar="FILENAME")
+args = parser.parse_args()
+
+# Scaling assuming variables are in kg
+if args.units == "kg":
+   massUnit = "kg"
+   fluxUnit = "kg yr$^{-1}$"
+   unitScaling = 1.0
+elif args.units == "Gt":
+   massUnit = "Gt"
+   fluxUnit = "Gt yr$^{-1}$"
+   unitScaling = 1.0e-12
+else:
+   sys.exit("Unknown mass unit")
+
+print("Using mass units of: ", massUnit)
+
+
+dataset = nc.Dataset(args.filename)
+
+# Read variables
+# convert everything to kg and years before unit conversion
+totalSubglacialWaterMass = \
+    dataset.variables['totalSubglacialWaterVolume'][:] * rhow * unitScaling
+melt = dataset.variables['totalBasalMeltInput'][:] * unitScaling * secyr
+distFluxMarine = dataset.variables['totalDistWaterFluxMarineMargin'][:] * unitScaling * secyr
+chnlFluxMarine = dataset.variables['totalChnlWaterFluxMarineMargin'][:] * unitScaling * secyr
+distFluxLand = dataset.variables['totalDistWaterFluxTerrestrialMargin'][:] * unitScaling * secyr
+chnlFluxLand = dataset.variables['totalChnlWaterFluxTerrestrialMargin'][:] * unitScaling * secyr
+chnlMelt = dataset.variables['totalChannelMelt'][:] * unitScaling * secyr
+flotFrac = dataset.variables['avgFlotationFraction'][:]
+lakeArea = dataset.variables['totalSubglacialLakeArea'][:] / 1000.0**2  # km^2
+lakeMass = dataset.variables['totalSubglacialLakeVolume'][:] * rhow * unitScaling
+grdArea = dataset.variables['groundedIceArea'][:] / 1000.0**2  # km^2
+
+deltat = dataset.variables['deltat'][:] / secyr  # in years
+yr = dataset.variables['daysSinceStart'][:] / 365.0
+
+subglacialWaterMassRate = np.zeros((len(melt),))
+
+for i in range(len(totalSubglacialWaterMass) - 1):
+    subglacialWaterMassRate[i] = ((totalSubglacialWaterMass[i+1] -
+                                   totalSubglacialWaterMass[i]) / deltat[i])
+
+# plot mass balance time-series
+fig, ax = plt.subplots(1, 1, layout='tight', figsize=(8,6))
+# input
+plt.plot(yr, melt, 'r:', label='basal melt')
+plt.plot(yr, chnlMelt, 'r--', label='channel melt')
+total_melt = melt + chnlMelt
+plt.plot(yr, total_melt, 'r-', label='total melt')
+
+# output
+plt.plot(yr, distFluxMarine, 'b--', label='marine sheet outflux')
+plt.plot(yr, distFluxLand, 'b:', label='land sheet outflux')
+plt.plot(yr, chnlFluxMarine, 'c--', label='marine chnl outflux')
+plt.plot(yr, chnlFluxLand, 'c:', label='land chnl outflux')
+total_outflux = distFluxMarine + distFluxLand + chnlFluxMarine + chnlFluxLand
+plt.plot(yr, total_outflux, 'b-', lw=2, label='total outflux')
+
+plt.plot(yr[1:-1], subglacialWaterMassRate[1:-1], 'g-', label='dV/dt')
+
+plt.plot(yr, total_melt - total_outflux, 'k', label='I-O')
+
+plt.legend(loc='best')
+plt.xlabel('Year')
+plt.ylabel(f'Mass flux ({fluxUnit})')
+
+# Plot other time-series of interest
+fig, axes = plt.subplots(2, 2, sharex=True, layout='tight',
+                         figsize=(10, 7))
+axes = axes.flatten()
+
+ax = 0
+axes[ax].plot(yr, flotFrac)
+axes[ax].set_ylabel('Flotation fraction')
+
+ax += 1
+axes[ax].plot(yr, totalSubglacialWaterMass)
+axes[ax].set_ylabel(f'Water mass ({massUnit})')
+
+ax += 1
+axes[ax].plot(yr, lakeArea)
+axes[ax].set_ylabel('Lake area (km$^2$)')
+# second axis for % area
+ax2 = axes[ax].twinx()
+ax2.plot(yr, lakeArea / grdArea, ':', color="blue")
+ax2.set_ylabel("Lake area percentage", color="blue")
+ax2.tick_params(axis="y", colors="blue")
+
+ax += 1
+axes[ax].plot(yr, lakeMass)
+axes[ax].set_ylabel(f'Lake mass ({massUnit})')
+
+for ax in axes:
+    ax.grid(True)
+    ax.set_xlabel("Year")
+
+
+plt.show()
+
