Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

feat: add test GZ plots for GLA12 #44

Merged
merged 1 commit into from
Jun 19, 2024
Merged

feat: add test GZ plots for GLA12 #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions GZ/calculate_GZ_ICESat2_ATL06.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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