From 0889d74fe4d4038444e91c0da92664b0c3564c72 Mon Sep 17 00:00:00 2001 From: Courtney Shafer Date: Fri, 24 Jan 2025 11:01:19 -0700 Subject: [PATCH 1/5] Add code to plot ocean TF within transect --- landice/output_processing_li/plot_transect.py | 59 +++++++++++++++++-- 1 file changed, 55 insertions(+), 4 deletions(-) diff --git a/landice/output_processing_li/plot_transect.py b/landice/output_processing_li/plot_transect.py index e4a5c7bac..48ee1b4be 100644 --- a/landice/output_processing_li/plot_transect.py +++ b/landice/output_processing_li/plot_transect.py @@ -2,8 +2,8 @@ # -*- coding: utf-8 -*- """ Plot bed and surface topography, surface speed, and -(optionally) temperature from MPAS netCDF along a transect. -@author: trevorhillebrand +(optionally) temperature and thermal forcing from MPAS netCDF along a transect. +@author: trevorhillebrand, edited by cashafer """ import numpy as np @@ -36,6 +36,9 @@ parser.add_option("--temperature", dest="interp_temp", action="store_true", help="interpolate temperature") +parser.add_option("--thermal_forcing", dest="interp_thermal_forcing", + action="store_true", help="interpolate thermal forcing") + for option in parser.option_list: if option.default != ("NO", "DEFAULT"): option.help += (" " if option.help else "") + "[default: %default]" @@ -66,12 +69,17 @@ if '-1' in times_list: times_list[times_list.index('-1')] = str(dataset.dimensions['Time'].size - 1) -# Cannot plot temperature for more than one time index. +# Cannot plot temperature or thermal forcing for more than one time index. if options.interp_temp and (len(times) > 1): print('Cannot plot temperature for more than one time index.' + ' Skipping temperature interpolation and plotting.') options.interp_temp = False +if options.interp_thermal_forcing and (len(times) > 1): + print('Cannot plot thermal forcing for more than one time index.' + + ' Skipping thermal forcing interpolation and plotting.') + options.interp_thermal_forcing = False + li_mask_ValueDynamicIce = 2 cellMask = dataset.variables['cellMask'][times,:] cellMask_dynamicIce = (cellMask & li_mask_ValueDynamicIce) // li_mask_ValueDynamicIce @@ -91,9 +99,17 @@ ' Skipping velocity plot.') plot_speed = False - if options.interp_temp: temperature = dataset.variables['temperature'][times,:] + +if options.interp_thermal_forcing: + # Ocean thermal forcing + thermal_forcing = dataset.variables['TFocean'][times,:] + thermal_forcing[thermal_forcing == 1.0e36] = np.nan # Large invalid TF values set to NaN + # Number of ocean layers + nISMIP6OceanLayers = dataset.dimensions['nISMIP6OceanLayers'].size + # Depths associated with thermal forcing field + ismip6shelfMelt_zOcean = dataset.variables['ismip6shelfMelt_zOcean'][:] bedTopo = dataset.variables["bedTopography"][0,:] print('Reading bedTopography from the first time level only. If multiple', @@ -185,6 +201,26 @@ temperature[i,:,lev]) temp_transect[:, lev] = temp_interpolant( np.vstack((x_interp, y_interp)).T) + + if options.interp_thermal_forcing: + ocean_layer_interfaces = np.zeros((len(thk_transect), nISMIP6OceanLayers + 1)) # Ocean layer depths + ocean_layer_interfaces[:,0] = 0. # First interface is the topmost z = 0 sea level interface + + # For the entire length of the transect, set the layer interface depths + for ii in range(len(thk_transect)): + ocean_layer_interfaces[ii,1:] = ismip6shelfMelt_zOcean[:] + ismip6shelfMelt_zOcean[0] + + # Create the thermal forcing transect using the TF interpolator at each ocean lev + thermal_forcing_transect = np.zeros((len(x_interp), nISMIP6OceanLayers)) + for ocean_lev in range(nISMIP6OceanLayers): + print(f'Interpolating ocean thermal forcing for ocean level {ocean_lev}') + thermal_forcing_interpolant = LinearNDInterpolator( + np.vstack((xCell, yCell)).T, + thermal_forcing[i,:,ocean_lev]) + thermal_forcing_transect[:, ocean_lev] = thermal_forcing_interpolant( + np.vstack( (x_interp, y_interp)).T) + + thickAx.plot(distance, bed_transect, color='black') @@ -195,6 +231,15 @@ cmap='YlGnBu_r', vmin=240., vmax=273.15) +if options.interp_thermal_forcing: + thermal_forcing_plot = thickAx.pcolormesh( np.tile(distance, (nISMIP6OceanLayers+1,1)).T, + ocean_layer_interfaces[:,:], thermal_forcing_transect[1:,:], + cmap='RdYlBu_r', vmin=-2, vmax=2 ) + thickAx.fill_between(distance, upper_surf_nan, lower_surf_nan, color='xkcd:ice blue') + thickAx.fill_between(distance, bed_transect, thickAx.get_ylim()[0], color='xkcd:greyish brown') + thickAx.grid(False) + + speedAx.set_xlabel('Distance (km)') speedAx.set_ylabel('Surface\nspeed (m/yr)') thickAx.set_ylabel('Elevation\n(m asl)') @@ -204,6 +249,11 @@ temp_cbar.set_label('Temperature (K)') temp_cbar.ax.tick_params(labelsize=12) +if options.interp_thermal_forcing: + tf_cbar = plt.colorbar(thermal_forcing_plot) + tf_cbar.set_label('Thermal Forcing (C)') + tf_cbar.ax.tick_params(labelsize=12) + if (len(times) > 1): time_cbar = plt.colorbar(cm.ScalarMappable(cmap='plasma'), ax=thickAx) time_cbar.ax.tick_params(labelsize=12) @@ -220,3 +270,4 @@ plt.show() dataset.close() + From 7e10582952586618ebd9bd66e59aeabd38fbf99a Mon Sep 17 00:00:00 2001 From: Courtney Shafer Date: Sat, 25 Jan 2025 14:23:32 -0700 Subject: [PATCH 2/5] Combine if statements for time index condition --- landice/output_processing_li/plot_transect.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/landice/output_processing_li/plot_transect.py b/landice/output_processing_li/plot_transect.py index 48ee1b4be..0f7b4fee3 100644 --- a/landice/output_processing_li/plot_transect.py +++ b/landice/output_processing_li/plot_transect.py @@ -70,14 +70,10 @@ times_list[times_list.index('-1')] = str(dataset.dimensions['Time'].size - 1) # Cannot plot temperature or thermal forcing for more than one time index. -if options.interp_temp and (len(times) > 1): - print('Cannot plot temperature for more than one time index.' + - ' Skipping temperature interpolation and plotting.') +if (options.interp_temp or options.interp_thermal_forcing) and (len(times) > 1): + print('Cannot plot temperature or thermal forcing for more than one time index.' + + ' Skipping temperature or TF interpolation and plotting.') options.interp_temp = False - -if options.interp_thermal_forcing and (len(times) > 1): - print('Cannot plot thermal forcing for more than one time index.' + - ' Skipping thermal forcing interpolation and plotting.') options.interp_thermal_forcing = False li_mask_ValueDynamicIce = 2 From 25604d6d8f8500cac875ec3c1b34f79c50816d29 Mon Sep 17 00:00:00 2001 From: Courtney Shafer Date: Sat, 1 Feb 2025 09:40:16 -0700 Subject: [PATCH 3/5] Fix cumulative distance to use multiple coords in transect --- landice/output_processing_li/plot_transect.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/landice/output_processing_li/plot_transect.py b/landice/output_processing_li/plot_transect.py index 0f7b4fee3..47bb5cc04 100644 --- a/landice/output_processing_li/plot_transect.py +++ b/landice/output_processing_li/plot_transect.py @@ -131,8 +131,8 @@ y = np.array([float(i) for i in options.y_coords.split(',')]) # increase sampling to match highest mesh resolution -total_distance, = np.cumsum( np.sqrt( np.diff(x)**2. + np.diff(y)**2. ) ) -n_samples = int(round(total_distance / np.min(dataset.variables["dcEdge"][:]))) +total_distance = np.cumsum( np.sqrt( np.diff(x)**2. + np.diff(y)**2. ) ) +n_samples = int(round(total_distance[-1] / np.min(dataset.variables["dcEdge"][:]))) x_interp = np.interp(np.linspace(0, len(x)-1, n_samples), np.linspace(0, len(x)-1, len(x)), x) y_interp = np.interp(np.linspace(0, len(y)-1, n_samples), @@ -230,11 +230,11 @@ if options.interp_thermal_forcing: thermal_forcing_plot = thickAx.pcolormesh( np.tile(distance, (nISMIP6OceanLayers+1,1)).T, ocean_layer_interfaces[:,:], thermal_forcing_transect[1:,:], - cmap='RdYlBu_r', vmin=-2, vmax=2 ) + cmap='YlOrRd', vmin=0, vmax=5.5) thickAx.fill_between(distance, upper_surf_nan, lower_surf_nan, color='xkcd:ice blue') thickAx.fill_between(distance, bed_transect, thickAx.get_ylim()[0], color='xkcd:greyish brown') thickAx.grid(False) - + thickAx.set_ylim([None, 750]) speedAx.set_xlabel('Distance (km)') speedAx.set_ylabel('Surface\nspeed (m/yr)') From 5cb1cecc5b07a4f291cc0fe4c6a5679a9f171008 Mon Sep 17 00:00:00 2001 From: Courtney Shafer Date: Mon, 3 Feb 2025 12:55:01 -0700 Subject: [PATCH 4/5] Add Delaunay triangulation and bounding box to speed up code --- landice/output_processing_li/plot_transect.py | 106 ++++++++++-------- 1 file changed, 60 insertions(+), 46 deletions(-) diff --git a/landice/output_processing_li/plot_transect.py b/landice/output_processing_li/plot_transect.py index 47bb5cc04..db13a2a1a 100644 --- a/landice/output_processing_li/plot_transect.py +++ b/landice/output_processing_li/plot_transect.py @@ -11,6 +11,7 @@ from netCDF4 import Dataset from optparse import OptionParser from scipy.interpolate import LinearNDInterpolator +from scipy.spatial import Delaunay import matplotlib.pyplot as plt from matplotlib.pyplot import cm @@ -44,15 +45,52 @@ option.help += (" " if option.help else "") + "[default: %default]" options, args = parser.parse_args() + +# Get list of times for plotting times_list = [i for i in options.times.split(',')] # list of string times for plotting times = [int(i) for i in options.times.split(',')] # list of integer time indices +# Load MPAS netCDF file dataset = Dataset(options.data_file, 'r') dataset.set_always_mask(False) + +# Get x and y coordinates +if options.coords_file is not None: + x = [] + y = [] + with open(options.coords_file, newline='') as csvfile: + reader = csv.reader(csvfile, delimiter=',') + for row in reader: + x.append(float(row[0])) + y.append(float(row[1])) + if [options.x_coords, options.y_coords] != [None, None]: + print('-c and -x/-y options were both provided. Reading from ', + f'{options.coords_file} and ignoring -x and -y settings.') + x = np.asarray(x) + y = np.asarray(y) +else: + x = np.array([float(i) for i in options.x_coords.split(',')]) + y = np.array([float(i) for i in options.y_coords.split(',')]) + +# Determine bounding box using max and min of user-provided coordinates +pad = 5000 +x_max = np.max(x) + pad +x_min = np.min(x) - pad +y_max = np.max(y) + pad +y_min = np.min(y) - pad + xCell = dataset.variables["xCell"][:] yCell = dataset.variables["yCell"][:] -nCells = dataset.dimensions['nCells'].size -areaCell = dataset.variables["areaCell"][:] + +indices = np.where( np.logical_and( + np.logical_and(xCell>=x_min, xCell<=x_max), + np.logical_and(yCell>=y_min, yCell<=y_max)) )[0] + +xCell = xCell[indices] +yCell = yCell[indices] + +tri_mesh = Delaunay( np.vstack((xCell,yCell)).T ) + layerThicknessFractions = dataset.variables["layerThicknessFractions"][:] nVertLevels = dataset.dimensions['nVertLevels'].size if "daysSinceStart" in dataset.variables.keys(): @@ -77,18 +115,19 @@ options.interp_thermal_forcing = False li_mask_ValueDynamicIce = 2 -cellMask = dataset.variables['cellMask'][times,:] +cellMask = dataset.variables['cellMask'][times,indices] cellMask_dynamicIce = (cellMask & li_mask_ValueDynamicIce) // li_mask_ValueDynamicIce # only take thickness of dynamic ice -thk = dataset.variables["thickness"][times,:] * cellMask_dynamicIce +thk = dataset.variables["thickness"][times,indices] * cellMask_dynamicIce + plot_speed = True # Include speed on non-dynamic ice to avoid interpolation artifacts. if "surfaceSpeed" in dataset.variables.keys(): - speed = dataset.variables["surfaceSpeed"][times,:] * 3600. * 24. * 365. + speed = dataset.variables["surfaceSpeed"][times,indices] * 3600. * 24. * 365. elif "surfaceSpeed" not in dataset.variables.keys() and \ all([ii in dataset.variables.keys() for ii in ['uReconstructX', 'uReconstructY']]): - speed = np.sqrt(dataset.variables["uReconstructX"][times,:,0]**2. + - dataset.variables["uReconstructY"][times,:,0]**2.) + speed = np.sqrt(dataset.variables["uReconstructX"][times,indices,0]**2. + + dataset.variables["uReconstructY"][times,indices,0]**2.) speed *= 3600. * 24. * 365. # convert from m/s to m/yr else: print('File does not contain surfaceSpeed or uReconstructX/Y.', @@ -96,43 +135,26 @@ plot_speed = False if options.interp_temp: - temperature = dataset.variables['temperature'][times,:] + temperature = dataset.variables['temperature'][times,indices] if options.interp_thermal_forcing: - # Ocean thermal forcing - thermal_forcing = dataset.variables['TFocean'][times,:] + # Extrapolated Ocean thermal forcing + thermal_forcing = dataset.variables['TFocean'][times,indices] thermal_forcing[thermal_forcing == 1.0e36] = np.nan # Large invalid TF values set to NaN # Number of ocean layers nISMIP6OceanLayers = dataset.dimensions['nISMIP6OceanLayers'].size # Depths associated with thermal forcing field ismip6shelfMelt_zOcean = dataset.variables['ismip6shelfMelt_zOcean'][:] + -bedTopo = dataset.variables["bedTopography"][0,:] +bedTopo = dataset.variables["bedTopography"][0,indices] print('Reading bedTopography from the first time level only. If multiple', 'times are needed, plot_transects.py will need to be updated.') -# Use coordinates from CSV file or -x -y options, but not both. -# CSV file takes precedent if both are present. -if options.coords_file is not None: - x = [] - y = [] - with open(options.coords_file, newline='') as csvfile: - reader = csv.reader(csvfile, delimiter=',') - - for row in reader: - x.append(float(row[0])) - y.append(float(row[1])) - if [options.x_coords, options.y_coords] != [None, None]: - print('-c and -x/-y options were both provided. Reading from ', - f'{options.coords_file} and ignoring -x and -y settings.') - x = np.asarray(x) - y = np.asarray(y) -else: - x = np.array([float(i) for i in options.x_coords.split(',')]) - y = np.array([float(i) for i in options.y_coords.split(',')]) + # increase sampling to match highest mesh resolution total_distance = np.cumsum( np.sqrt( np.diff(x)**2. + np.diff(y)**2. ) ) -n_samples = int(round(total_distance[-1] / np.min(dataset.variables["dcEdge"][:]))) +n_samples = int(round(total_distance[-1] / np.min(dataset.variables["dcEdge"][indices]))) x_interp = np.interp(np.linspace(0, len(x)-1, n_samples), np.linspace(0, len(x)-1, len(x)), x) y_interp = np.interp(np.linspace(0, len(y)-1, n_samples), @@ -154,12 +176,11 @@ plt.rcParams.update({'font.size': 16}) -bed_interpolant = LinearNDInterpolator(np.vstack((xCell, yCell)).T, bedTopo) +bed_interpolant = LinearNDInterpolator(tri_mesh, bedTopo) bed_transect = bed_interpolant(np.vstack((x_interp, y_interp)).T) for i, time in enumerate(times): - thk_interpolant = LinearNDInterpolator( - np.vstack((xCell, yCell)).T, thk[i,:]) + thk_interpolant = LinearNDInterpolator(tri_mesh, thk[i,:]) thk_transect = thk_interpolant(np.vstack((x_interp, y_interp)).T) lower_surf = np.maximum( -910. / 1028. * thk_transect, bed_transect) lower_surf_nan = lower_surf.copy() # for plotting @@ -171,8 +192,7 @@ thickAx.plot(distance, upper_surf_nan, color=timeColors[i]) if plot_speed: - speed_interpolant = LinearNDInterpolator( - np.vstack((xCell, yCell)).T, speed[i,:]) + speed_interpolant = LinearNDInterpolator(tri_mesh, speed[i,:]) speed_transect = speed_interpolant(np.vstack((x_interp, y_interp)).T) speed_transect[thk_transect == 0.] = np.nan speedAx.plot(distance, speed_transect, color=timeColors[i]) @@ -192,9 +212,7 @@ temp_transect = np.zeros((len(x_interp), nVertLevels)) for lev in range(nVertLevels): print(f'Interpolating temperature for level {lev}') - temp_interpolant = LinearNDInterpolator( - np.vstack((xCell, yCell)).T, - temperature[i,:,lev]) + temp_interpolant = LinearNDInterpolator(tri_mesh, temperature[i,:,lev]) temp_transect[:, lev] = temp_interpolant( np.vstack((x_interp, y_interp)).T) @@ -208,15 +226,11 @@ # Create the thermal forcing transect using the TF interpolator at each ocean lev thermal_forcing_transect = np.zeros((len(x_interp), nISMIP6OceanLayers)) + for ocean_lev in range(nISMIP6OceanLayers): print(f'Interpolating ocean thermal forcing for ocean level {ocean_lev}') - thermal_forcing_interpolant = LinearNDInterpolator( - np.vstack((xCell, yCell)).T, - thermal_forcing[i,:,ocean_lev]) - thermal_forcing_transect[:, ocean_lev] = thermal_forcing_interpolant( - np.vstack( (x_interp, y_interp)).T) - - + thermal_forcing_interpolant = LinearNDInterpolator(tri_mesh, thermal_forcing[i,:,ocean_lev]) + thermal_forcing_transect[:, ocean_lev] = thermal_forcing_interpolant(np.vstack( (x_interp,y_interp) ).T) thickAx.plot(distance, bed_transect, color='black') From 5a277f8572ffdd745381a0759e338720533e26e5 Mon Sep 17 00:00:00 2001 From: Courtney <36939203+cshafer@users.noreply.github.com> Date: Tue, 4 Feb 2025 10:23:02 -0500 Subject: [PATCH 5/5] Fix spacing and indexing to be PEP8-compliant Co-authored-by: Trevor Hillebrand --- landice/output_processing_li/plot_transect.py | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/landice/output_processing_li/plot_transect.py b/landice/output_processing_li/plot_transect.py index db13a2a1a..d7c69db85 100644 --- a/landice/output_processing_li/plot_transect.py +++ b/landice/output_processing_li/plot_transect.py @@ -89,7 +89,7 @@ xCell = xCell[indices] yCell = yCell[indices] -tri_mesh = Delaunay( np.vstack((xCell,yCell)).T ) +tri_mesh = Delaunay( np.vstack((xCell, yCell)).T ) layerThicknessFractions = dataset.variables["layerThicknessFractions"][:] nVertLevels = dataset.dimensions['nVertLevels'].size @@ -115,19 +115,19 @@ options.interp_thermal_forcing = False li_mask_ValueDynamicIce = 2 -cellMask = dataset.variables['cellMask'][times,indices] +cellMask = dataset.variables['cellMask'][times, indices] cellMask_dynamicIce = (cellMask & li_mask_ValueDynamicIce) // li_mask_ValueDynamicIce # only take thickness of dynamic ice -thk = dataset.variables["thickness"][times,indices] * cellMask_dynamicIce +thk = dataset.variables["thickness"][times, indices] * cellMask_dynamicIce plot_speed = True # Include speed on non-dynamic ice to avoid interpolation artifacts. if "surfaceSpeed" in dataset.variables.keys(): - speed = dataset.variables["surfaceSpeed"][times,indices] * 3600. * 24. * 365. + speed = dataset.variables["surfaceSpeed"][times, indices] * 3600. * 24. * 365. elif "surfaceSpeed" not in dataset.variables.keys() and \ all([ii in dataset.variables.keys() for ii in ['uReconstructX', 'uReconstructY']]): - speed = np.sqrt(dataset.variables["uReconstructX"][times,indices,0]**2. + - dataset.variables["uReconstructY"][times,indices,0]**2.) + speed = np.sqrt(dataset.variables["uReconstructX"][times, indices, 0]**2. + + dataset.variables["uReconstructY"][times, indices, 0]**2.) speed *= 3600. * 24. * 365. # convert from m/s to m/yr else: print('File does not contain surfaceSpeed or uReconstructX/Y.', @@ -135,11 +135,11 @@ plot_speed = False if options.interp_temp: - temperature = dataset.variables['temperature'][times,indices] + temperature = dataset.variables['temperature'][times, indices] if options.interp_thermal_forcing: # Extrapolated Ocean thermal forcing - thermal_forcing = dataset.variables['TFocean'][times,indices] + thermal_forcing = dataset.variables['TFocean'][times, indices] thermal_forcing[thermal_forcing == 1.0e36] = np.nan # Large invalid TF values set to NaN # Number of ocean layers nISMIP6OceanLayers = dataset.dimensions['nISMIP6OceanLayers'].size @@ -147,7 +147,7 @@ ismip6shelfMelt_zOcean = dataset.variables['ismip6shelfMelt_zOcean'][:] -bedTopo = dataset.variables["bedTopography"][0,indices] +bedTopo = dataset.variables["bedTopography"][0, indices] print('Reading bedTopography from the first time level only. If multiple', 'times are needed, plot_transects.py will need to be updated.') @@ -229,8 +229,8 @@ for ocean_lev in range(nISMIP6OceanLayers): print(f'Interpolating ocean thermal forcing for ocean level {ocean_lev}') - thermal_forcing_interpolant = LinearNDInterpolator(tri_mesh, thermal_forcing[i,:,ocean_lev]) - thermal_forcing_transect[:, ocean_lev] = thermal_forcing_interpolant(np.vstack( (x_interp,y_interp) ).T) + thermal_forcing_interpolant = LinearNDInterpolator(tri_mesh, thermal_forcing[i, :, ocean_lev]) + thermal_forcing_transect[:, ocean_lev] = thermal_forcing_interpolant(np.vstack( (x_interp, y_interp) ).T) thickAx.plot(distance, bed_transect, color='black') @@ -242,8 +242,8 @@ vmin=240., vmax=273.15) if options.interp_thermal_forcing: - thermal_forcing_plot = thickAx.pcolormesh( np.tile(distance, (nISMIP6OceanLayers+1,1)).T, - ocean_layer_interfaces[:,:], thermal_forcing_transect[1:,:], + thermal_forcing_plot = thickAx.pcolormesh(np.tile(distance, (nISMIP6OceanLayers+1, 1)).T, + ocean_layer_interfaces[:, :], thermal_forcing_transect[1:, :], cmap='YlOrRd', vmin=0, vmax=5.5) thickAx.fill_between(distance, upper_surf_nan, lower_surf_nan, color='xkcd:ice blue') thickAx.fill_between(distance, bed_transect, thickAx.get_ylim()[0], color='xkcd:greyish brown')