Figure 7¶

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

pej0_plot_length = 200   # This determines how many vertices from the terminus should be plotted (200 vertices = 10 km)

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')
            ax1[0, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['b'], color='xkcd:brown')
            ax1[1, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['u'], color='xkcd:light green')
            ax1[2, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['pe_ignore_dslope'], color='xkcd:light red')
            ax1[3, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['j0_ignore_dslope'], color='xkcd:light blue')
            ax1[4, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['udiff'], color='xkcd:gray')
            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')
            # 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=30)
        else:
            ax1[1, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['u'], color='xkcd:dark green')
            ax1[2, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['pe_ignore_dslope'], color='xkcd:dark red')
            ax1[3, 0].plot(results[basin_no][key]['d'], results[basin_no][key]['j0_ignore_dslope'], color='xkcd:dark blue')
            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=3, 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 
    ax1[0, 0].set_title('Basin {}'.format(basin_no))
    ax1[0, 0].set_ylabel('Surface (cyan) and bed (brown) \n elevation (m)')
    ax1[1, 0].set_ylabel('Speed 1996 (m/yr)')
    ax1[2, 0].set_ylabel(r'$\frac{P_e}{\ell}$ (1/m)')
    ax1[3, 0].set_ylabel(r'$J_0$ (m/yr)')
    ax1[4, 0].set_xlabel('Distance from terminus (km)')
    ax1[4, 0].set_ylabel('Speed change 1996-2018 (m/yr)')
    ax1[4, 0].set_xlim([0, results[basin_no]['avg']['d'][-1]])
    axbig.set_xlabel(r'$\frac{P_e}{\ell}$ (1/m)')
    axbig.set_ylabel(r'$J_0$ (m/yr)')
    axbig.set_title('Dot spacing: 50 m; \n Big dot indicates the first non-NaN value (the one closest to the terminus)')
    
    plt.savefig('../data/results/Austfonna_glacier{}.png'.format(basin_no))

colors = np.array([[5,113,176,255],
                   [171,229,255,255],
                   [255,240,189,255],
                   [222,130,38,255],
                   [202,0,32,255],
                   [220,22,177,255],
                   [132,17,177,255]])
colors = colors / 255
nodes = (np.array([-500, -300, 100, 500, 1500, 2500, 3500]) + 500) / 4000
mycmap = LinearSegmentedColormap.from_list("mycmap", list(zip(nodes, colors)))
mycmap