<jupyter_start><jupyter_text>Figures 3-4This notebook produces Figures 3-4 and other auxiliary figures.        Get data(1) from Felikson et al.: Inland limits to diffusion of thinning along Greenland Ice Sheet outlet glaciers, available at https://zenodo.org/record/4284759. Redistributed under the CC-BY-4.0 License (see the previous link for details).Reference: <NAME>., <NAME>., <NAME>., <NAME>., & <NAME>. (2021). Steep Glacier Bed Knickpoints Mitigate Inland Thinning in Greenland. Geophysical Research Letters, 48(2), 1–10. https://doi.org/10.1029/2020GL090112<jupyter_code>%%bash

if [ ! -d ../data/Felikson2021/ ]; then
    wget https://zenodo.org/record/4284759/files/netcdfs.zip -P ../data/Felikson2021/ --no-verbose
    unzip ../data/Felikson2021/netcdfs.zip -d ../data/Felikson2021/
    rm -rf ../data/Felikson2021/netcdfs.zip
fi<jupyter_output><empty_output><jupyter_text>(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 manually made the file for the 2018-1998 speed difference.References: 1. <NAME>., <NAME>, and <NAME> (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/6II6VW8LLWJ72. Gardner, <NAME>., <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, and <NAME> (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<jupyter_code>%%bash

if [ ! -f ../data/GRE_G0240_1998_v.tif ]; then
    wget http://www.geo.cornell.edu/eas/gstudent/wz278/paperdraft2021/GRE_G0240_1998_v.tif -P ../data/ --no-verbose
    wget http://www.geo.cornell.edu/eas/gstudent/wz278/paperdraft2021/GRE_G0240_diff-2018-1998_v.tif -P ../data/ --no-verbose
fi<jupyter_output><empty_output><jupyter_text>Analysis<jupyter_code>import pejzero
import rasterio
from netCDF4 import Dataset
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from pathlib import Path<jupyter_output><empty_output><jupyter_text>File locations. You can replace the number in `/glacier0207.nc` for processing other glacier data. In the paper, **Figure 3** uses `/glacier0207.nc` and **Figure 4** uses `/glacier0277.nc`.<jupyter_code>netcdf_dir = '../data/Felikson2021'
glacier_file = netcdf_dir + '/glacier0207.nc'
speed_file = '../data/GRE_G0240_1998_v.tif'
vdiff_file = '../data/GRE_G0240_diff-2018-1998_v.tif'<jupyter_output><empty_output><jupyter_text>Open the data set and show 6 primary flowlines.<jupyter_code>ds = Dataset(glacier_file, 'r')
flowline_groups, _ = pejzero.get_flowline_groups(ds)
primary_flowlines = [i for i in flowline_groups if 'iter' not in i.path]
primary_flowlines
# Here's how to access variables in the primary_flowlines object:
# primary_flowlines[0]['geometry']['surface']['GIMP']['nominal']['h'][:]<jupyter_output><empty_output><jupyter_text>Calculate $P_e/\ell$ and $J_0$ for each flowline, and calculate the average of all 6 flowlines:<jupyter_code>results = {}

with rasterio.open(speed_file) as speed_data, rasterio.open(vdiff_file) as vdiff_data:
    for flowline_group in primary_flowlines:
        data_group = pejzero.cal_pej0_for_each_flowline(flowline_group, speed_data, vdiff_data)
        if data_group is not None:
            results[flowline_group.name] = data_group
        
results['avg'] = pejzero.cal_avg_for_each_basin(results)
for key in results:
    print(key)<jupyter_output>flowline03
flowline04
flowline05
flowline06
flowline07
flowline08
avg<jupyter_text>Now we can plot **Figure 3**:<jupyter_code>pej0_plot_length = 200   # This determines how many vertices from the terminus should be plotted (200 vertices = 10 km)

matplotlib.rc('font', size=20)

# ==== Create axes (5 in the left and a big one in the right)
fig, ax3 = plt.subplots(5, 2, sharex=True, figsize=(24, 20))
gs = ax3[1, 1].get_gridspec()
for ax in ax3[:, 1]:
    ax.remove()
axbig = fig.add_subplot(gs[1:4, 1])

# ==== Plot data
for key in results:
    
    if key != 'avg':
        ax3[0, 0].plot(results[key]['d'], results[key]['s'], color='xkcd:aquamarine')
        ax3[0, 0].plot(results[key]['d'], results[key]['b'], color='xkcd:brown')
        ax3[1, 0].plot(results[key]['d'], results[key]['u'], color='xkcd:light green')
        ax3[2, 0].plot(results[key]['d'], results[key]['pe_ignore_dslope'], color='xkcd:light red')
        ax3[3, 0].plot(results[key]['d'], results[key]['j0_ignore_dslope'], color='xkcd:light blue')
        ax3[4, 0].plot(results[key]['d'], results[key]['udiff_sm'], color='xkcd:light grey')
        axbig.plot(results[key]['pe_ignore_dslope'][:pej0_plot_length], results[key]['j0_ignore_dslope'][:pej0_plot_length], '.-', color='xkcd:light purple')
        # plot first non-NaN value (the one closest to the terminus)
        axbig.plot(next(x for x in results[key]['pe_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)),
                   next(x for x in results[key]['j0_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)), '.', color='xkcd:light purple', markersize=30)
    else:
        ax3[1, 0].plot(results[key]['d'], results[key]['u'], color='xkcd:dark green', linewidth=3)
        ax3[2, 0].plot(results[key]['d'], results[key]['pe_ignore_dslope'], color='xkcd:dark red', linewidth=3)
        ax3[3, 0].plot(results[key]['d'], results[key]['j0_ignore_dslope'], color='xkcd:dark blue', linewidth=3)
        ax3[4, 0].plot(results[key]['d'], results[key]['udiff_sm'], color='xkcd:dark grey', linewidth=3)
        axbig.plot(results[key]['pe_ignore_dslope'][:pej0_plot_length], results[key]['j0_ignore_dslope'][:pej0_plot_length], '.-', color='xkcd:dark purple', linewidth=3, markersize=10)
        # plot first non-NaN value (the one closest to the terminus)
        axbig.plot(next(x for x in results[key]['pe_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)),
                   next(x for x in results[key]['j0_ignore_dslope'][:pej0_plot_length] if not np.isnan(x)), '.', color='xkcd:dark purple', markersize=30)
        
# ==== Add supplemental text 
ax3[0, 0].set_title(Path(glacier_file).stem)
ax3[0, 0].set_ylabel('Elevantion (m): \n Surface (cyan) \n bed (brown)')
ax3[1, 0].set_ylabel('Speed 1998 (m yr$^{-1}$)')
ax3[2, 0].set_ylabel(r'$\frac{P_e}{\ell}$ (m$^{-1}$)')
ax3[3, 0].set_ylabel(r'$J_0$ (m yr$^{-1}$)')
ax3[4, 0].set_xlabel('Distance from terminus (km)')
ax3[4, 0].set_ylabel('Speed change \n 1998–2018 (m yr$^{-1}$)')
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)')

plt.savefig('../data/results/' + Path(glacier_file).stem + '.pdf')<jupyter_output><empty_output><jupyter_text>To run the analysis in a batch mode for more than one glacier, uncomment the following cell using `Pe-J0-Greenland-SingleBasin.py`. The default setting prints PNG files but can be switched to pdf easily.<jupyter_code># %%bash

# filelist=`ls --color=never ../data/Felikson2021/*`
# for i in $filelist; do 
#     python Pe-J0-Greenland-SingleBasin.py $i
# done<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15932250068851556, 'fraction_numerical': 0.0497108234646103, 'mean_word_length': 4.987633965375103, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '607069', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2915, 'n_tokens_neox': 2656, 'n_words': 613, 'path': 'workflows/Fig3-4.ipynb', 'repo_name': 'whyjz/sensitive-glacier', 'size': '7357'}