Skip to content

Commit

Permalink
feat: add test GZ plots for GLA12 (#44)
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsutterley
tsutterley authored Jun 19, 2024
1 parent abae6a5 commit ec71838
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 22 deletions.
1 change: 1 addition & 0 deletions GZ/calculate_GZ_ICESat2_ATL06.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,7 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
# get output directory from input file
if OUTPUT_DIRECTORY is None:
OUTPUT_DIRECTORY = INPUT_FILE.parent
OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
# set the hemisphere flag based on ICESat-2 granule
HEM = set_hemisphere(GRAN)

Expand Down
33 changes: 18 additions & 15 deletions GZ/calculate_GZ_ICESat2_ATL11.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,6 +204,8 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
# get output directory from input file
if OUTPUT_DIRECTORY is None:
OUTPUT_DIRECTORY = INPUT_FILE.parent
# create output directory if it doesn't exist
OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
# file format for associated auxiliary files
file_format = '{0}_{1}_{2}_{3}{4}_{5}{6}_{7}_{8}{9}.h5'
# set the hemisphere flag based on ICESat-2 granule
Expand Down Expand Up @@ -527,19 +529,19 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
# ocean tide height for scaling model
tide_mean = np.mean(tide_ocean['AT'][i,:],axis=1)
# tide_mean = tide_ocean['AT'].data[i,0]
h_tide = np.ma.array(tide_ocean['AT'].data[i,c] - tide_mean,
h_tide = np.ma.array(tide_ocean['AT'].data[i,c],
fill_value=tide_ocean['AT'].fill_value)
h_tide.mask = tide_ocean['AT'].mask[i,c] | tide_mean.mask
# inverse-barometer response
ib_mean = np.mean(IB['AT'][i,:],axis=1)
# ib_mean = IB['AT'].data[i,0]
h_ib = np.ma.array(IB['AT'].data[i,c] - ib_mean,
h_ib = np.ma.array(IB['AT'].data[i,c],
fill_value=IB['AT'].fill_value)
h_ib.mask = IB['AT'].mask[i,c] | ib_mean.mask
# deflection from mean land ice height in grounding zone
dh_gz = h_gz - h_mean
# quasi-freeboard: WGS84 elevation - geoid height
QFB = h_gz - geoid_h[i]
QFB = h_gz - h_tide - h_ib - geoid_h[i]
# ice thickness from quasi-freeboard and densities
w_thick = QFB*rho_w/(rho_w-rho_ice)
# fit with a hard piecewise model to get rough estimate of GZ
Expand All @@ -566,13 +568,13 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
# set tide values for testing
TIDE = []
i0 = 0 if sco else -1
tplus = h_tide[i0] + h_ib[i0]
tplus = (h_tide[i0] - tide_mean) + (h_ib[i0] - ib_mean)
# 1,3: use tide range values from Padman (2002)
# 2,4: use tide range values from model+ib
TIDE.append([1.2,-3.0,3.0])
TIDE.append([tplus,tplus-0.3,tplus+0.3])
TIDE.append([tplus, tplus-0.3, tplus+0.3])
TIDE.append([1.2,-3.0,3.0])
TIDE.append([tplus,tplus-0.3,tplus+0.3])
TIDE.append([tplus, tplus-0.3, tplus+0.3])
# iterate through tests
for grz,tide in zip(GRZ,TIDE):
# fit physical elastic model
Expand Down Expand Up @@ -641,12 +643,13 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
# add to test plot
if PLOT:
# plot height differences
l, = ax1.plot(ref_pt['AT'][i],dh_gz-PdH[0],'.-',ms=1.5,lw=0,
label=f'Cycle {CYCLE}')
l, = ax1.plot(ref_pt['AT'][i], dh_gz-PdH[0], '.-',
ms=1.5, lw=0, label=f'Cycle {CYCLE}')
# plot downstream tide and IB
hocean = tide_ocean['AT'].data[i0,c] - mean_tide
# hocean += IB['AT'].data[i0,c] - mean_ib
ax1.axhline(hocean,color=l.get_color(),lw=3.0,ls='--')
ax1.axhline(hocean, color=l.get_color(), lw=3.0,
ls='--')
# set valid plot flag
valid_plot = True

Expand All @@ -667,14 +670,14 @@ def calculate_GZ_ICESat2(base_dir, INPUT_FILE,
# add to test plot
if PLOT:
# plot elastic deformation model
ax1.plot(ref_pt['AT'][iout],PEMODEL-PdH[0],
color='0.3',lw=2,zorder=9)
ax1.plot(ref_pt['AT'][iout], PEMODEL-PdH[0],
color='0.3', lw=2, zorder=9)
# plot scaled elastic deformation model
ax1.plot(ref_pt['AT'][iout],flexure[iout,c]-mean_tide,
color='0.8',lw=2,zorder=10)
ax1.plot(ref_pt['AT'][iout], flexure[iout,c]-mean_tide,
color='0.8', lw=2, zorder=10)
# plot grounding line location
ax1.axvline(GZrpt,color=l.get_color(),
ls='--',dashes=(8,4))
ax1.axvline(GZrpt, color=l.get_color(),
ls='--', dashes=(8,4))

# make final plot adjustments and save to file
if valid_plot:
Expand Down
63 changes: 56 additions & 7 deletions GZ/calculate_GZ_ICESat_GLA12.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@
# attempt imports
fiona = gz.utilities.import_dependency('fiona')
h5py = gz.utilities.import_dependency('h5py')
plt = gz.utilities.import_dependency('matplotlib.pyplot')
pyproj = gz.utilities.import_dependency('pyproj')
geometry = gz.utilities.import_dependency('shapely.geometry')
timescale = gz.utilities.import_dependency('timescale')
Expand Down Expand Up @@ -187,6 +188,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
REANALYSIS=None,
SEA_LEVEL=False,
GEOID=None,
PLOT=False,
MODE=0o775
):

Expand Down Expand Up @@ -241,8 +243,12 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
# get output directory from input file
if OUTPUT_DIRECTORY is None:
OUTPUT_DIRECTORY = INPUT_FILE.parent
# create output directory if it doesn't exist
OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True)
# full path to output file
OUTPUT_FILE = OUTPUT_DIRECTORY.joinpath(FILENAME)
# plot counter
plot_number = 0

# open the HDF5 file for reading
fid = gz.io.multiprocess_h5py(INPUT_FILE, mode='r')
Expand Down Expand Up @@ -284,13 +290,16 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
# Elevation (height above TOPEX/Poseidon ellipsoid in meters)
d_elev = fid['Data_40HZ']['Elevation_Surfaces']['d_elev'][:]
fv = fid['Data_40HZ']['Elevation_Surfaces']['d_elev'].fillvalue
# retide the elevation data
d_ocElv = fid['Data_40HZ']['Geophysical']['d_ocElv'][:]
d_ocElv[d_ocElv == fv] = 0.0
# mask for reducing to valid values
quality_summary = fid1['Data_40HZ']['Quality']['d_qa_sum'][:]
# get the transform for converting to the latest ITRF
transform = gz.crs.tp_itrf2008_to_wgs84_itrf2020()
# transform the data to WGS84 ellipsoid in ITRF2020
lon, lat, data, tdec = transform.transform(d_lon, d_lat,
d_elev[:] + sat_corr + bias_corr, ts.year)
d_elev[:] + sat_corr + bias_corr + d_ocElv[:], ts.year)
# mask invalid values
elev = np.ma.array(data, fill_value=fv)
elev.mask = (d_elev == fv) | np.isnan(elev.data)
Expand All @@ -306,7 +315,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
a2 = (PRD,RL,VAR,RGTP,ORB,INST,CYCL,TRK,SEG,GRAN,TYPE)
f2 = OUTPUT_DIRECTORY.joinpath(file_format.format(*a2))
elif (FILE_FORMAT == 'generic'):
file2 = f'{INPUT_FILE.stem}_{VAR}_{INPUT_FILE.suffix}'
file2 = f'{INPUT_FILE.stem}_{VAR}{INPUT_FILE.suffix}'
f2 = OUTPUT_DIRECTORY.joinpath(file2)
# extract grounding zone mask
ice_gz = np.zeros((n_40HZ),dtype=bool)
Expand Down Expand Up @@ -335,7 +344,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
f4 = OUTPUT_DIRECTORY.joinpath(file_format.format(*a4))
elif TIDE_MODEL and (FILE_FORMAT == 'generic'):
VAR = f'{TIDE_MODEL}_TIDES'
file4 = f'{INPUT_FILE.stem}_{VAR}_{INPUT_FILE.suffix}'
file4 = f'{INPUT_FILE.stem}_{VAR}{INPUT_FILE.suffix}'
f4 = OUTPUT_DIRECTORY.joinpath(file4)
if TIDE_MODEL:
fid4 = gz.io.multiprocess_h5py(f4, mode='r')
Expand All @@ -358,7 +367,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
f5 = OUTPUT_DIRECTORY.joinpath(file_format.format(*a5))
elif REANALYSIS and (FILE_FORMAT == 'generic'):
VAR = 'DAC' if (REANALYSIS == 'DAC') else f'{REANALYSIS}_IB'
file5 = f'{INPUT_FILE.stem}_{VAR}_{INPUT_FILE.suffix}'
file5 = f'{INPUT_FILE.stem}_{VAR}{INPUT_FILE.suffix}'
f5 = OUTPUT_DIRECTORY.joinpath(file5)
if REANALYSIS:
fid5 = gz.io.multiprocess_h5py(f5, mode='r')
Expand All @@ -377,7 +386,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
f6 = OUTPUT_DIRECTORY.joinpath(file_format.format(*a6))
elif SEA_LEVEL and (FILE_FORMAT == 'generic'):
VAR = 'DAC' if (REANALYSIS == 'DAC') else f'{REANALYSIS}_IB'
file6 = f'{INPUT_FILE.stem}_{VAR}_{INPUT_FILE.suffix}'
file6 = f'{INPUT_FILE.stem}_{VAR}{INPUT_FILE.suffix}'
f6 = OUTPUT_DIRECTORY.joinpath(file6)
if SEA_LEVEL:
fid6 = gz.io.multiprocess_h5py(f6, mode='r')
Expand All @@ -394,7 +403,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
f7 = OUTPUT_DIRECTORY.joinpath(file_format.format(*a7))
elif GEOID and (FILE_FORMAT == 'generic'):
VAR = f'{GEOID}_GEOID'
file7 = f'{INPUT_FILE.stem}_{VAR}_{INPUT_FILE.suffix}'
file7 = f'{INPUT_FILE.stem}_{VAR}{INPUT_FILE.suffix}'
f7 = OUTPUT_DIRECTORY.joinpath(file7)
if GEOID:
fid7 = gz.io.multiprocess_h5py(f7, mode='r')
Expand Down Expand Up @@ -470,7 +479,7 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
# deflection from mean land ice height in grounding zone
dh_gz = h_gz - h_mean
# quasi-freeboard: WGS84 elevation - geoid height
QFB = h_gz - (h_geoid + mdt[i])
QFB = h_gz - h_tide - h_ib - (h_geoid + mdt[i])
# ice thickness from quasi-freeboard and densities
w_thick = QFB*rho_w/(rho_w-rho_ice)
# fit with a hard piecewise model to get rough estimate of GZ
Expand Down Expand Up @@ -549,6 +558,41 @@ def calculate_GZ_ICESat(base_dir, INPUT_FILE,
# Data_GZ['d_H_ice'].append(PT)
# Data_GZ['d_delta_h'].append(PdH)

# add to test plot
if PLOT:
fig1,ax1 = plt.subplots(num=1,figsize=(13,7))
# plot height differences
l, = ax1.plot(dist, dh_gz - PdH[0],
'.-', ms=1.5, lw=0)
# ax1.plot(dist, PEMODEL)
# plot downstream tide and IB
ax1.axhline(tplus, color='0.3',
lw=3.0, ls='--')
# plot grounding line location
ax1.axvline(PGZ[0], color=l.get_color(),
ls='--', dashes=(8,4))
# adjust figure
fig1.subplots_adjust(left=0.05, right=0.97,
bottom=0.04, top=0.96)
# create plot file of flexural zone
VAR = f'{TIDE_MODEL}_GZ_TIDES_{plot_number:d}'
if (FILE_FORMAT == 'standard'):
a8 = (PRD,RL,VAR,RGTP,ORB,INST,CYCL,TRK,SEG,GRAN,TYPE)
f8 = OUTPUT_DIRECTORY.joinpath(file_format.format(*a8))
elif (FILE_FORMAT == 'generic'):
file8 = f'{INPUT_FILE.stem}_{VAR}{INPUT_FILE.suffix}'
f8 = OUTPUT_DIRECTORY.joinpath(file8)
# log output plot file
output_plot_file = f8.with_suffix('.png')
logging.info(str(output_plot_file))
fig1.savefig(output_plot_file, dpi=240, format='png',
metadata={'Title':pathlib.Path(sys.argv[0]).name})
# clear all figure axes
plt.cla()
plt.clf()
# add to plot counter
plot_number += 1

# if no valid grounding zone fit has been found
# skip saving variables and attributes
if not valid_fit:
Expand Down Expand Up @@ -808,6 +852,10 @@ def arguments():
parser.add_argument('--geoid','-G',
metavar='GEOID', type=str,
help='Geoid height model to use in correction')
# create test plots
parser.add_argument('--plot','-P',
default=False, action='store_true',
help='Create plots of flexural zone')
# verbose will output information about each output file
parser.add_argument('--verbose','-V',
default=False, action='store_true',
Expand Down Expand Up @@ -840,6 +888,7 @@ def main():
REANALYSIS=args.reanalysis,
SEA_LEVEL=args.sea_level,
GEOID=args.geoid,
PLOT=args.plot,
MODE=args.mode)

# run main program
Expand Down

0 comments on commit ec71838

Please sign in to comment.