Skip to content

Commit cc30d75

Browse files
trhilletrhille
committed
Merge pull request #417 from 'matthewhoffman/landice/regional_stats_plot'
Add multiple new scripts that were used in ISMIP6-2300 data processing.
2 parents 5cc1b97 + 19eb8fb commit cc30d75

File tree

7 files changed

+1298
-17
lines changed

7 files changed

+1298
-17
lines changed

‎landice/mesh_tools_li/adjust_iceshelf_geometry.py‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
#!/usr/bin/env python
2+
'''
3+
Script to fine tune ice-shelf geometry for Antarctica ISMIP6 meshes:
4+
1. Smooth ice thickness across ice shelves. The first row of cells adjacent to GL is left unmodified in this step.
5+
2. Adjust bathymetry to prevent spurious grounding from thickness adjustment by modifying
6+
bathymetry to maintain original water column thickness.
7+
3. Lower bathymetry beneath ice shelves and in open ocean by a uniform amount. This reduces spurious
8+
grounding line advance from transient model dH/dt and geometry inconsistencies. 5-20 meters appears to be a good amount.
9+
Note that adjusting the seaward cells at the grounding line still impacts the velocity solver due to its
10+
subgrid parameterization of friction near the grounding line (the subgrid position evaluation of the grounding
11+
line position depends on the bed elevation at the first ocean node).
12+
13+
See beginning of script for options to adjust.
14+
15+
Input file is copied with a _modifiedBathymetry.nc suffix before being modified.
16+
17+
Matt Hoffman, 9/11/2022
18+
'''
19+
20+
from __future__ import absolute_import, division, print_function, unicode_literals
21+
22+
import sys
23+
import numpy as np
24+
from netCDF4 import Dataset
25+
from optparse import OptionParser
26+
import matplotlib.pyplot as plt
27+
import shutil
28+
29+
30+
# ---- options to adjust (could become command line options ----
31+
32+
nSmoothing = 3 # number of rounds of smoothing to apply. 3 was the best match at 8km to 30 years of thickness evolution on Ross and FRIS. Set to 0 to disable
33+
34+
elevationDrop = 5.0 # height in meters to drop non-grounded marine bed elevation
35+
36+
# ------------------------
37+
38+
39+
print("** Gathering information. (Invoke with --help for more details. All arguments are optional)")
40+
parser = OptionParser(description=__doc__)
41+
parser.add_option("-f", dest="fileName", help="filename to modify", metavar="FILENAME")
42+
options, args = parser.parse_args()
43+
44+
rhoi=910.0
45+
rhosw=1028.0
46+
47+
# Make a copy of file being modified and work on the copy
48+
outFileName = options.fileName.split('.nc')[0] + '_modifiedBathymetry.nc'
49+
shutil.copy(options.fileName, outFileName)
50+
51+
f = Dataset(outFileName, 'r+')
52+
53+
nCells = len(f.dimensions['nCells'])
54+
mu = f.variables['muFriction'][0,:]
55+
thickness = f.variables['thickness'][0,:]
56+
bedTopography = f.variables['bedTopography'][0,:]
57+
cOnC= f.variables['cellsOnCell'][:]
58+
nEOnC = f.variables['nEdgesOnCell'][:]
59+
xCell = f.variables['xCell'][:]
60+
yCell = f.variables['yCell'][:]
61+
lowerSurface = -rhoi/rhosw*thickness # only works for floating cells
62+
WCT = lowerSurface - bedTopography
63+
64+
floatMask = ((thickness*910/1028+bedTopography)<0.0)*(thickness>0.0)
65+
groundMask = ((thickness*910/1028+bedTopography)>0.0)*(thickness>0.0)
66+
67+
# Find row of cells forming the ocean-side of the grounding line (first floating cells)
68+
oGLmask = np.zeros((nCells,), 'i')
69+
for c in range(nCells):
70+
if floatMask[c] == 1:
71+
for n in range(nEOnC[c]):
72+
if groundMask[cOnC[c,n]-1] == 1:
73+
oGLmask[c] = True
74+
break
75+
76+
# Find floating cells at calving front
77+
oMgnMask = np.zeros((nCells,),'i')
78+
for c in range(nCells):
79+
if floatMask[c] == 1:
80+
for n in range(nEOnC[c]):
81+
if thickness[cOnC[c,n]-1] == 0.0:
82+
oMgnMask[c] = True
83+
break
84+
85+
# Update ocean margin mask to also include floating cells one layer inward from calving front
86+
# This is to avoid using the margin cells as input to the smoothing, because they might be
87+
# non-dynamic and have an unrealistic thickness.
88+
oMgnMaskOld = oMgnMask.copy()
89+
for c in range(nCells):
90+
if floatMask[c] == 1:
91+
for n in range(nEOnC[c]):
92+
if oMgnMaskOld[cOnC[c,n]-1] == 1:
93+
oMgnMask[c] = True
94+
95+
96+
97+
mask = np.logical_and(np.logical_not(np.logical_or(oMgnMask, oGLmask)), floatMask) # mask of the floating area excluding the boundaries
98+
99+
#fig, axs = plt.subplots(1,1, figsize=(9,9), num=1)
100+
#plt.scatter(xCell, yCell, s=2, c=mask*1+floatMask*1)
101+
#plt.colorbar()
102+
#plt.show()
103+
104+
# smooth over the mask
105+
106+
def smooth(H):
107+
Hold = H.copy()
108+
for c in range(nCells):
109+
if mask[c] == 1:
110+
Hsum = 0.0
111+
cnt = nEOnC[c]
112+
for n in range(cnt):
113+
Hsum += Hold[cOnC[c,n]-1]
114+
Hsum += Hold[c]
115+
cnt += 1
116+
H[c] = Hsum / float(cnt)
117+
return H
118+
119+
for i in range(nSmoothing):
120+
print(f"start smoothing {i}")
121+
thickness = smooth(thickness)
122+
print("done")
123+
124+
f.variables['thickness'][0,:] = thickness
125+
126+
# Keep WCT
127+
newLowerSurface = -rhoi/rhosw*thickness
128+
for c in range(nCells):
129+
if floatMask[c] == 1:
130+
bedTopography[c] = newLowerSurface[c] - WCT[c]
131+
bedTopography[c] -= elevationDrop
132+
f.variables['bedTopography'][0,:] = bedTopography
133+
134+
f.close()
135+
136+
print("saved")
137+
138+
fig, axs = plt.subplots(1,1, figsize=(9,9), num=1)
139+
plt.scatter(xCell, yCell, s=2, c=thickness, vmin=200, vmax=1000.0, cmap='RdBu')
140+
plt.colorbar()
141+
plt.show()

‎landice/mesh_tools_li/bulldoze_missing_troughs.py‎

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
#!/usr/bin/env python
2+
'''
3+
Script to bulldoze marine troughs.
4+
5+
At a number of locations around Antarctica, it is clear that the
6+
portion of subglacial troughs seaward of the grounding line are not continuous with
7+
the portion beneath grounded ice. This is expected from the relative dearth of
8+
sub-ice-shelf bathymetric measurements compared with the high density of radar
9+
measurements beneath grounded ice. The coarse representation of bahthymetry on the
10+
seaward side of the grounding line in these locations leads to strong grounding line
11+
advance at these outlet glaciers and associated grounded mass gain.
12+
13+
In this script, a number of trough locations to be deepened are specified by x/y
14+
coordinate pairs for the start and end of the trough. The smoothTrough function
15+
is called for each, which applies a parabolic transverse cross section over a centerline
16+
trough depth that increases linearly between a minimum and maximum value.
17+
Note that as currently written, the algorithm only works for troughs oriented southwest to
18+
northeast in grid coordinates!
19+
20+
The location of the points defining the troughs and the amount to bulldoze is
21+
hardcoded after the smoothTrough function. Adjustments should be made there.
22+
23+
Input file is copied with a _bulldozedTroughs.nc suffix before being modified.
24+
25+
Matt Hoffman, Fall 2022
26+
'''
27+
28+
from __future__ import absolute_import, division, print_function, unicode_literals
29+
30+
import sys
31+
import numpy as np
32+
from netCDF4 import Dataset
33+
from optparse import OptionParser
34+
import matplotlib.pyplot as plt
35+
import shutil
36+
37+
38+
print("** Gathering information. (Invoke with --help for more details. All arguments are optional)")
39+
parser = OptionParser(description=__doc__)
40+
parser.add_option("-f", dest="fileName", help="filename to modify", metavar="FILENAME")
41+
options, args = parser.parse_args()
42+
43+
rhoi=910.0
44+
rhosw=1028.0
45+
46+
# Make a copy of file being modified and work on the copy
47+
outFileName = options.fileName.split('.nc')[0] + '_bulldozedTroughs.nc'
48+
shutil.copy(options.fileName, outFileName)
49+
print("Writing output to:", outFileName)
50+
f = Dataset(outFileName, 'r+')
51+
52+
nCells = len(f.dimensions['nCells'])
53+
mu = f.variables['muFriction'][0,:]
54+
thickness = f.variables['thickness'][0,:]
55+
bedTopography = f.variables['bedTopography'][0,:]
56+
cOnC= f.variables['cellsOnCell'][:]
57+
nEOnC = f.variables['nEdgesOnCell'][:]
58+
xCell = f.variables['xCell'][:]
59+
yCell = f.variables['yCell'][:]
60+
xEdge = f.variables['xEdge'][:]
61+
yEdge = f.variables['yEdge'][:]
62+
dcEdge = f.variables['dcEdge'][:]
63+
lowerSurface = -rhoi/rhosw*thickness # only works for floating cells
64+
WCT = lowerSurface - bedTopography
65+
66+
floatMask = ((thickness*910/1028+bedTopography)<0.0)*(thickness>0.0)
67+
#groundMask = ((thickness*910/1028+bedTopography)>0.0)*(thickness>0.0)
68+
69+
def smoothTrough(WCT, p1, p2, minWCT, maxWCT):
70+
WCTnew = WCT.copy()
71+
72+
# find representative dCell for this area
73+
ind = np.nonzero( (xEdge>p1[0]) * (xEdge<p2[0]) * (yEdge>p1[1]) * (yEdge<p1[2]) )[0]
74+
dCell = dcEdge[ind].mean()
75+
print(f"using dCell={dCell}")
76+
77+
mask = (WCT<maxWCT) * (floatMask==1) * (xCell>(p1[0]-2*dCell)) * (xCell<(p2[0]+2*dCell)) * (yCell>(p1[1]-2*dCell)) * (yCell<(p2[1]+2*dCell))
78+
ind = np.nonzero(mask==1)[0]
79+
length = ( (p2[0]-p1[0])**2 + (p2[1]-p1[1])**2 )**0.5 # length along flowline
80+
m1 = (p2[1]-p1[1]) / (p2[0]-p1[0]) # xy slope of flowline
81+
b1 = p1[1] - m1*p1[0] # y int of flowline
82+
print(f'# cells={len(ind)}, length={length}, m={m1}')
83+
m2 = -1.0/m1 # slope of all transverse lines
84+
for c in ind:
85+
# find width and center at this position
86+
b2 = yCell[c] - m2 * xCell[c] # y int of transverse line
87+
dist2orthogline = np.absolute(m2*xCell + -1*yCell + b2) / (m2**2 + (-1)**2)**0.5
88+
ind2 = np.nonzero((mask==1) * (dist2orthogline<(1.5*dCell)))[0] # get the points within a certain distance of the orthog line
89+
dist2centerline = np.absolute(m1*xCell[ind2] + -1*yCell[ind2] + b1) / (m1**2 + (-1)**2)**0.5
90+
#dist2centerline = np.absolute( (p2[0]-p1[0]) * (p1[1] - yCell[ind2]) - (p1[0]-xCell[ind2])*(p2[1]-p1[1])) / ( (p2[0]-p1[0])**2 + (p2[1]-p1[1])**2)**0.5
91+
dist2centerline *= np.sign((m1*xCell[ind2]+b1) - yCell[ind2]) # make signed distance, assuming NE direction
92+
print(dist2centerline)
93+
width = dist2centerline.max()-dist2centerline.min()
94+
95+
# find frac x along center line
96+
# first need intersection of longit. and ortho lines
97+
xint = (b2-b1)/(m1-(-1.0/m1))
98+
mag = (xint-p1[0]) / (p2[0]-p1[0]) * (maxWCT-minWCT) + minWCT
99+
100+
widthOrgnIdx = ind2[np.argmin(dist2centerline)] #idx of westmost
101+
widthOrgn = (xCell[widthOrgnIdx], yCell[widthOrgnIdx]) # position of westmost
102+
widthEndIdx = ind2[np.argmax(dist2centerline)] #idx of westmost
103+
widthEnd = (xCell[widthEndIdx], yCell[widthEndIdx]) # position of westmost
104+
widthMidpt = ((widthOrgn[0]+widthEnd[0])/2.0, (widthOrgn[1]+widthEnd[1])/2.0)
105+
print(xCell[c], yCell[c], widthOrgn)
106+
#dist2Orgn = ((widthOrgn[0]-xCell[c])**2 + (widthOrgn[1]-yCell[c])**2)**0.5
107+
#fracWidth = (dist2Orgn / (width) - 0.5) * 2.0 # [-1,1] range
108+
dist2Mid = ((widthMidpt[0]-xCell[c])**2 + (widthMidpt[1]-yCell[c])**2)**0.5
109+
fracWidth = dist2Mid / (0.5*(width+2*dCell)) # [-1,1] range
110+
fracWidth = min(fracWidth, 0.95)
111+
print(f'found {len(ind2)} points along this orthogonal; width={width}; fracWidth={fracWidth}, mag={mag}')
112+
z = mag * (-1.0 * fracWidth**2 + 1.0)
113+
WCTnew[c] = max(z, WCT[c]) # don't make WCT thinner than the original
114+
print(f'old z={WCT[c]}, new z={z}')
115+
116+
return WCTnew
117+
118+
119+
# Define maximum and minimum water column thickness along tough centerline
120+
maxWCT = 300.0
121+
minWCT = 50.0
122+
123+
# Specify start and end x/y coordinates for each trough to be bulldozed
124+
WCTnew = WCT.copy()
125+
p1=[-1261073.6, 137133.0]; p2=[-1173862.4, 199332.15] # Rutford
126+
WCTnew = smoothTrough(WCTnew, p1, p2, minWCT, maxWCT)
127+
p1=[-1315459.94, 202346.75]; p2=[-1257332.3, 291441.3] # Carlson
128+
WCTnew = smoothTrough(WCTnew, p1, p2, minWCT, maxWCT)
129+
p1=[-1513787.6, 311996.4]; p2=[-1386639.4, 337107.2] # Evans
130+
WCTnew = smoothTrough(WCTnew, p1, p2, minWCT, maxWCT)
131+
p1=[-34186, 1983479]; p2=[-9070, 2047275] # Fimbul
132+
WCTnew = smoothTrough(WCTnew, p1, p2, minWCT, maxWCT)
133+
p1=[-1064519, 206317]; p2=[-1055596, 265431] # lil guy - note different max/minWCT
134+
WCTnew = smoothTrough(WCTnew, p1, p2, 40.0, 150.0)
135+
136+
f.variables['bedTopography'][0,:] = lowerSurface - WCTnew
137+
138+
f.close()
139+
140+
print("saved")
141+
s=12
142+
ind=(floatMask==1)
143+
fig, axs = plt.subplots(1,2, figsize=(14,9), num=1, sharex=True, sharey=True)
144+
plt.sca(axs[0])
145+
plt.scatter(xCell[ind], yCell[ind], s=s, c=WCT[ind], vmin=0, vmax=100.0, cmap='RdBu')
146+
plt.colorbar()
147+
plt.sca(axs[1])
148+
plt.scatter(xCell[ind], yCell[ind], s=s, c=WCTnew[ind], vmin=0, vmax=100.0, cmap='RdBu')
149+
plt.plot(p1[0], p1[1], 'k*')
150+
plt.plot(p2[0], p2[1], 'ko')
151+
axs[0].set_aspect('equal', 'box')
152+
axs[1].set_aspect('equal', 'box')
153+
plt.colorbar()
154+
plt.show()

‎landice/mesh_tools_li/set_von_mises_stress_map.py‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
#!/usr/bin/env python
2+
'''
3+
Simple script for creating an input file with region-by-region values for
4+
the von Mises stress threshold. This can be used to assign optimal
5+
regional values identified through a tuning process.
6+
7+
A region mask file is required as input. Values to assign are hardcoded
8+
below. There is not error checking that the number of values match the
9+
number of regions, so use with care.
10+
11+
The script outputs a file called von_mises_calving_parameters.nc with the
12+
assigned regional values.
13+
14+
Matt Hoffman, 9/19/2022
15+
'''
16+
17+
from __future__ import absolute_import, division, print_function, unicode_literals
18+
19+
import sys
20+
import numpy as np
21+
from netCDF4 import Dataset
22+
from optparse import OptionParser
23+
import matplotlib.pyplot as plt
24+
25+
parser = OptionParser(description=__doc__)
26+
parser.add_option("-n", dest="fileRegions", help="region file name.", metavar="FILENAME")
27+
options, args = parser.parse_args()
28+
29+
f = Dataset(options.fileRegions, 'r')
30+
regionCellMasks = f.variables['regionCellMasks'][:]
31+
nRegions = len(f.dimensions['nRegions'])
32+
nCells = len(f.dimensions['nCells'])
33+
34+
fout = Dataset("von_mises_calving_parameters.nc", 'w')
35+
fout.createDimension('nCells', nCells)
36+
fout.createDimension('Time', None)
37+
grdVM = fout.createVariable('groundedVonMisesThresholdStress', 'd', ('Time', 'nCells',))
38+
fltVM = fout.createVariable('floatingVonMisesThresholdStress', 'd', ('Time', 'nCells',))
39+
40+
values=[
41+
125.0,
42+
200.0,
43+
150.0,
44+
300.0, #300
45+
46+
200.0, #225?
47+
350.0, # 400?
48+
400.0,
49+
125.0, # 130-400
50+
51+
300.0, #400?
52+
300.0, #300-400
53+
300.0, # 300
54+
200.0, # 100?
55+
56+
125.0,
57+
125.0,#?
58+
125.0,#120-400
59+
125.0,#125-300
60+
]
61+
62+
grdVM[:]=100.0e3
63+
fltVM[:]=100.0e3
64+
for r in range(nRegions):
65+
mask = np.nonzero(regionCellMasks[:,r] == 1)[0]
66+
grdVM[0, mask] = values[r] * 1000.0
67+
fltVM[0, mask] = values[r] * 1000.0
68+
69+
fout.close()
70+
f.close()

0 commit comments

Comments
 (0)