Figure 7
Contents
Figure 7¶
This notebook produces Figure 7 and other auxiliary figures.
Get data¶
(1) From the Ice Thickness Models Intercomparison eXperiment (ITMIX) data set (Farinotti et al., 2017), available at https://doi.org/10.5905/ethz-1007-92. Redistributed under the CC BY-NC-SA 4.0 License (see the previous link for details).
Reference: Farinotti, D., Brinkerhoff, D. J., Clarke, G. K. C., Fürst, J. J., Frey, H., Gantayat, P., Gillet-Chaulet, F., Girard, C., Huss, M., Leclercq, P. W., Linsbauer, A., Machguth, H., Martin, C., Maussion, F., Morlighem, M., Mosbeux, C., Pandit, A., Portmann, A., Rabatel, A., …Andreassen, L. M. (2017). How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment. Cryosphere, 11(2), 949–970. https://doi.org/10.5194/tc-11-949-2017
%%bash
if [ ! -d ../data/ITMIX/ ]; then
wget -P ../data/ITMIX/ --content-disposition --no-verbose 'https://github.com/whyjz/pejzero/releases/download/v0.1/ITMIX_dataset.zip'
unzip ../data/ITMIX/ITMIX_dataset.zip -d ../data/ITMIX/
unzip ../data/ITMIX/01_ITMIX_input_data.zip -d ../data/ITMIX/01_ITMIX_input_data/
unzip ../data/ITMIX/02_ITMIX_results.zip -d ../data/ITMIX/02_ITMIX_results/
unzip ../data/ITMIX/03_ITMIX_measured_thickness.zip -d ../data/ITMIX/03_ITMIX_measured_thickness/
rm -rf ../data/ITMIX/*.zip
fi
(2) From ITS_LIVE: A NASA MEaSUREs project to provide automated, low latency, global glacier flow and elevation change datasets. Original data available at https://its-live.jpl.nasa.gov/. We use the 2018 speed mosaics.
References:
Gardner, A. S., M. A. Fahnestock, and T. A. Scambos (2019) [Accessed October 21, 2021]. ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities. Data archived at National Snow and Ice Data Center. https://doi.org/10.5067/6II6VW8LLWJ7
Gardner, A. S., G. Moholdt, T. Scambos, M. Fahnstock, S. Ligtenberg, M. van den Broeke, and J. Nilsson (2018). Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years, Cryosphere, 12(2), 521–547. https://doi.org/10.5194/tc-12-521-2018
%%bash
if [ ! -f ../data/SRA_G0240_2018_v_EPSG32633.tif ]; then
wget https://github.com/whyjz/pejzero/releases/download/v0.1/SRA_G0240_2018_v_EPSG32633.tif -P ../data/ --no-verbose
fi
(3) The flowline shapefile available at ../data/flowline_statpts4_clean_EPSG32633.shp
for flowline vertices, which is included in the data folder.
Analysis¶
import geopandas as gpd
import rasterio as rio
import numpy as np
import pejzero
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import os
Here are all file locations:
flowline_file = '../data/flowline_statpts4_clean_EPSG32633.shp'
bed_file = '../data/ITMIX/02_ITMIX_results/Farinotti_Austfonna_bedrock.asc'
surface_file = '../data/ITMIX/01_ITMIX_input_data/Austfonna/02_surface_Austfonna_2007_UTM33.asc'
speed_file = '../data/ITMIX/01_ITMIX_input_data/Austfonna/06_speed_Austfonna_1995-1996_UTM33.asc'
speed2018_file = '../data/SRA_G0240_2018_v_EPSG32633.tif'
(1) Data ingest¶
Sample all the necessary measurements at the vertex locations:
flowlines = gpd.read_file(flowline_file)
d_collection = []
s_collection = []
b_collection = []
u_collection = []
udiff_collection = []
with rio.open(surface_file) as src_surface, rio.open(bed_file) as src_bed, rio.open(speed_file) as src_speed, rio.open(speed2018_file) as src_speed2018:
for idx, row in flowlines.iterrows():
# ==== creating distance marks
d_stop = (len(row.geometry.coords) - 1)* 50
d = np.linspace(0, d_stop, len(row.geometry.coords))
d_collection.append(d)
# ==== sample surface elevation
surface_gen = src_surface.sample(row.geometry.coords)
s = np.array([float(record) for record in surface_gen])
s[s < -9998] = np.nan
s_collection.append(s)
# ==== sample bed elevation
bed_gen = src_bed.sample(row.geometry.coords)
b = np.array([float(record) for record in bed_gen])
b[b < -9998] = np.nan
b_collection.append(b)
# ==== sample glacier speed
speed_gen = src_speed.sample(row.geometry.coords)
u = np.array([float(record) for record in speed_gen])
u[u < 0] = np.nan
u_collection.append(u)
# ==== sample and calculate glacier speed change
speed2018_gen = src_speed2018.sample(row.geometry.coords)
u2018 = np.array([float(record) for record in speed2018_gen])
u2018[u2018 < 0] = np.nan
udiff = u2018 - u
udiff_collection.append(udiff)
flowlines['d'] = d_collection
flowlines['s'] = s_collection
flowlines['b'] = b_collection
flowlines['u'] = u_collection
flowlines['udiff'] = udiff_collection
flowlines[:8] # shows the first 8 rows for a data structure overview
glacier_id | geometry | d | s | b | u | udiff | |
---|---|---|---|---|---|---|---|
0 | 11 | LINESTRING (661744.501 8812864.548, 661772.209... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [124.80580139160156, 125.9896011352539, 126.51... | [6.389999866485596, 6.630000114440918, 6.76000... | [9.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0... | [2.063974380493164, 1.063974380493164, 1.06397... |
1 | 12 | LINESTRING (664882.963 8810146.886, 664950.156... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [105.45079803466797, 107.12020111083984, 107.8... | [-66.0, -75.98999786376953, -78.91999816894531... | [12.0, 13.0, 13.0, 13.0, 13.0, 14.0, 14.0, 14.... | [0.3671684265136719, 0.06629180908203125, 0.06... |
2 | 13 | LINESTRING (669704.623 8810585.218, 669770.209... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [130.61459350585938, 132.28590393066406, 132.4... | [-107.91999816894531, -113.83999633789062, -11... | [18.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.... | [-6.172802925109863, -5.395359039306641, -5.77... |
3 | 14 | LINESTRING (673035.951 8808972.154, 673095.381... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [127.32119750976562, 127.9738998413086, 128.56... | [-41.369998931884766, -45.77000045776367, -47.... | [14.0, 14.0, 14.0, 14.0, 14.0, 14.0, 14.0, 14.... | [-0.8997259140014648, -0.8997259140014648, -0.... |
4 | 15 | LINESTRING (676753.012 8808446.155, 676815.893... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [144.21249389648438, 146.13330078125, 147.1508... | [-46.41999816894531, -50.880001068115234, -53.... | [9.0, 8.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, ... | [5.6007080078125, 6.6007080078125, 7.600708007... |
5 | 16 | LINESTRING (680680.473 8808481.221, 680771.006... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [142.48159790039062, 142.54119873046875, 142.3... | [-84.48999786376953, -88.70999908447266, -89.3... | [8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, ... | [27.496936798095703, 27.496936798095703, 27.49... |
6 | 31 | LINESTRING (709733.163 8827112.114, 709698.322... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... |
7 | 32 | LINESTRING (712538.492 8829608.857, 712479.960... | [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0,... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... | [nan, nan, nan, nan, nan, nan, nan, nan, nan, ... |
(2) Derive and ¶
Fig5.ipynb
, we have saved the results as mega_results_Austfonna.h5
and have commented out the cell below. If you wish to repeat/verify the mega_results_Austfonna.h5
file.
# results = {'1': {}, '3': {}, '4': {}, '5': {}, '6': {}, '7': {}, '10': {}, '17': {}}
# for idx, row in flowlines.iterrows():
# data_group = pejzero.cal_pej0_for_each_flowline_raw(row['d'], row['s'], row['b'], row['u'], savgol_winlength=251)
# if data_group is not None:
# data_group['udiff'] = row['udiff']
# udiff_sm = pejzero.my_savgol_filter(row['udiff'], window_length=151, polyorder=1, deriv=0, delta=50, mode='interp')
# data_group['udiff_sm'] = udiff_sm
# basin_no = str(row['glacier_id'] // 10)
# flowline_no = str(row['glacier_id'] % 10)
# results[basin_no][flowline_no] = data_group
# #### Calculate average
# for key in results:
# data_group = results[key]
# avg = pejzero.cal_avg_for_each_basin(data_group)
# results[key]['avg'] = avg
# pejzero.save_pej0_results(results, "../data/results/mega_results_Austfonna.h5")
Now we read the data from the HDF5 file. If you rerun the cell above and already get the results
variable, you can safely comment this cell out.
results = pejzero.load_pej0_results("../data/results/mega_results_Austfonna.h5")
(3) plot for single glacier basin¶
You can change basin_no_list
to plot more basins. The result PNGs will be saved in ../data/results/single_basins/
.
pej0_plot_length = 200 # This determines how many vertices from the terminus should be plotted (200 vertices = 10 km)
matplotlib.rc('font', size=24)
matplotlib.rc('axes', linewidth=2)
basin_no_list = ['7']
# basin_no_list = ['1', '3', '4', '5', '6', '7', '10', '17']
for basin_no in basin_no_list:
# ==== Create axes (5 in the left and a big one in the right)
fig, ax1 = plt.subplots(5, 2, sharex=True, figsize=(24, 20))
gs = ax1[1, 1].get_gridspec()
for ax in ax1[:, 1]:
ax.remove()
axbig = fig.add_subplot(gs[1:4, 1])
# ==== Plot data
for key in results[basin_no]:
if key != 'avg':
ax1[0, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['s'], color='xkcd:aquamarine', linewidth=2)
ax1[0, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['b'], color='xkcd:brown', linewidth=2)
ax1[1, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['u'], color='xkcd:light green', linewidth=2)
ax1[2, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['pe_ignore_dslope'], color='xkcd:light red', linewidth=2)
ax1[3, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['j0_ignore_dslope'], color='xkcd:light blue', linewidth=2)
ax1[4, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['udiff'], color='xkcd:gray', linewidth=2)
axbig.plot(results[basin_no][key]['pe_ignore_dslope'][:pej0_plot_length], results[basin_no][key]['j0_ignore_dslope'][:pej0_plot_length], '.-', color='xkcd:light purple', linewidth=2)
# plot first non-NaN value (the one closest to the terminus)
axbig.plot(next(x for x in results[basin_no][key]['pe_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)),
next(x for x in results[basin_no][key]['j0_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)), '.', color='xkcd:light purple', markersize=25)
else:
ax1[1, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['u'], color='xkcd:dark green', linewidth=4)
ax1[2, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['pe_ignore_dslope'], color='xkcd:dark red', linewidth=4)
ax1[3, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['j0_ignore_dslope'], color='xkcd:dark blue', linewidth=4)
axbig.plot(results[basin_no][key]['pe_ignore_dslope'][:pej0_plot_length], results[basin_no][key]['j0_ignore_dslope'][:pej0_plot_length], '.-', color='xkcd:dark purple', linewidth=4, markersize=10)
# plot first non-NaN value (the one closest to the terminus)
axbig.plot(next(x for x in results[basin_no][key]['pe_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)),
next(x for x in results[basin_no][key]['j0_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)), '.', color='xkcd:dark purple', markersize=30)
# ==== Add supplemental text and labels
letter_specs = {'fontsize': 30, 'fontweight': 'bold', 'va': 'top', 'ha': 'center'}
ax1[0, 0].set_title('Basin {}'.format(basin_no))
ax1[0, 0].set_ylabel('Elevantion (m): \n Surface (cyan) \n bed (brown)')
ax1[0, 0].text(0.04, 0.95, 'A', transform=ax1[0, 0].transAxes, **letter_specs)
ax1[1, 0].set_ylabel('Speed 1996 \n (m yr$^{-1}$)')
ax1[1, 0].text(0.96, 0.95, 'B', transform=ax1[1, 0].transAxes, **letter_specs)
ax1[2, 0].set_ylabel(r'$\frac{P_e}{\ell}$ (m$^{-1}$)')
ax1[2, 0].text(0.96, 0.95, 'C', transform=ax1[2, 0].transAxes, **letter_specs)
ax1[3, 0].set_ylabel(r'$J_0$ (m yr$^{-1}$)')
ax1[3, 0].text(0.96, 0.95, 'D', transform=ax1[3, 0].transAxes, **letter_specs)
ax1[4, 0].set_xlabel('Distance from terminus (km)')
ax1[4, 0].set_ylabel('Speed change \n 1996–2018 (m yr$^{-1}$)')
ax1[4, 0].set_xlim([0, results[basin_no]['avg']['d'][-1]])
ax1[4, 0].text(0.96, 0.95, 'E', transform=ax1[4, 0].transAxes, **letter_specs)
axbig.set_xlabel(r'$\frac{P_e}{\ell}$ (m$^{-1}$)')
axbig.set_ylabel(r'$J_0$ (m yr$^{-1}$)')
axbig.set_title('Dot spacing: 50 m; \n Big dot indicates the first non-NaN value \n (closest to the terminus)')
axbig.text(0.03, 0.985, 'F', transform=axbig.transAxes, **letter_specs)
pe_labels = ['{:.6f}'.format(x) for x in axbig.get_xticks()]
axbig.set_xticklabels(pe_labels, rotation=45)
outdir = '../data/results/single_basins/'
if not os.path.exists(outdir):
os.makedirs(outdir)
plt.savefig(outdir + 'Austfonna_glacier{}.png'.format(basin_no))
<ipython-input-8-9197152a31aa>:60: UserWarning: FixedFormatter should only be used together with FixedLocator
axbig.set_xticklabels(pe_labels, rotation=45)
Visualization¶
We use the same color map from Figure 5 but stretch it differently.
# Version 2 (revised submission, based on the roma colormap https://github.com/GenericMappingTools/gmt/blob/master/share/cpt/roma.cpt)
colors = np.array([[56,156,198,230],
[100,198,213,210],
[192,234,195,180],
[207,229,168,190],
[200,180,85,230],
[126,23,0,255],
[63,12,0,255],
[32,6,0,255]])
colors = colors / 255
nodes = (np.array([-500, -300, 0, 100, 500, 1500, 2500, 3500]) + 500) / 4000
mycmap = LinearSegmentedColormap.from_list("mycmap", list(zip(nodes, colors)))
mycmap
Finally, we can make Figure 7 by color coding each average glacier flowline using mycmap
on the profile and the
# alias (to keep consistent with Figure 5)
mega_results = results
matplotlib.rc('font', size=20)
matplotlib.rc('axes', linewidth=2)
fig, ax4 = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 3]}, figsize=(9, 12))
for key in mega_results:
z_max = np.nanmax(mega_results[key]['avg']['udiff_sm'])
z_min = np.nanmin(mega_results[key]['avg']['udiff_sm'])
z_value = z_max if abs(z_max) > abs(z_min) else z_min
#### The color map is stretched between 0 and 2000 m/yr! ####
z_value_scaled = z_value / 2000
# z_value_scaled = (z_value + 500) / 4000
#############################################################
rgba = mycmap(z_value_scaled)
# rgba2 = np.array(rgba) / 1.5
ax4[0].plot(mega_results[key]['avg']['d'], mega_results[key]['avg']['udiff_sm'], color=rgba)
ax4[0].set_xlim([3, 20])
ax4[0].set_xlabel('Distance from the terminus (km)')
label_top = ax4[0].set_ylabel('Speed change \n 1996-2018 \n (m yr$^{-1}$)')
ax4[0].tick_params(width=2, length=5)
ax4[1].plot(mega_results[key]['avg']['pe_ignore_dslope'][:100], mega_results[key]['avg']['j0_ignore_dslope'][:100], color=rgba, linewidth=2)
# Not plotting small dots here for better readibility
# ax4[1].plot(mega_results[key]['avg']['pe_ignore_dslope'][:100], mega_results[key]['avg']['j0_ignore_dslope'][:100], '.', color=rgba, markersize=6)
ax4[1].plot(next(x for x in mega_results[key]['avg']['pe_ignore_dslope'] if not np.isnan(x)),
next(x for x in mega_results[key]['avg']['j0_ignore_dslope'] if not np.isnan(x)), '.', color=rgba, markersize=10)
ax4[1].set_xlabel(r'$\frac{P_e}{\ell}$ (m$^{-1}$)')
ax4[1].set_ylabel(r'$J_0$ (m yr$^{-1}$)')
ax4[1].tick_params(width=2, length=5)
# Fine tune label location
if key == '5':
ax4[1].text(next(x for x in mega_results[key]['avg']['pe_ignore_dslope'] if not np.isnan(x)) - 0.00002,
next(x for x in mega_results[key]['avg']['j0_ignore_dslope'] if not np.isnan(x)) - 0.2, str(key), size=20)
else:
ax4[1].text(next(x for x in mega_results[key]['avg']['pe_ignore_dslope'] if not np.isnan(x)) + 0.00001,
next(x for x in mega_results[key]['avg']['j0_ignore_dslope'] if not np.isnan(x)) - 0.2, str(key), size=20)
letter_specs = {'fontsize': 30, 'fontweight': 'bold', 'va': 'top', 'ha': 'center'}
ax4[0].text(0.04, 0.95, 'A', transform=ax4[0].transAxes, **letter_specs)
ax4[1].text(0.04, 0.08, 'B', transform=ax4[1].transAxes, **letter_specs)
# ==== Place a color bar
cbaxes = ax4[1].inset_axes([0.34, 0.15, 0.6, 0.1])
norm = matplotlib.colors.Normalize(vmin=0, vmax=2000)
bounds = [0, 1000, 2000]
fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=mycmap),
cax=cbaxes, orientation='horizontal', label='Speed change (m yr$^{-1}$)', ticks=bounds)
cbaxes.tick_params(width=2, length=5)
fig.savefig('../data/results/Fig7.pdf', bbox_extra_artists=(label_top,), bbox_inches='tight')