Figure 2 script

Figure 2 script

To reproduce this figure, make sure you have downloaded all necessary input files (velocity maps and static terrain geometries) from https://doi.org/10.17605/OSF.IO/HE7YR and have updated the Vx and Vy columns in notebooks/manifest.csv with the downloaded file paths before starting the analysis.

import glaft
import pandas as pd
import matplotlib as mpl
from matplotlib.ticker import FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
# font and linewidth settings
font = {'size'   : 14}
mpl.rc('font', **font)
axes_settings = {'linewidth'   : 2}
mpl.rc('axes', **axes_settings)

# read and select data
df = pd.read_csv('../manifest.csv', dtype=str)
in_shp = '/home/jovyan/Projects/PX_comparison/shapefiles/bedrock_V2.shp'
selected_cases = df.loc[[126, 130, 134, 97, 73, 28]]
selected_cases
Date Duration (days) Template size (px) Template size (m) Pixel spacing (px) Pixel spacing (m) Prefilter Subpixel Software Vx Vy
126 LS8-20180304-20180405 32 64 960 4 60 None pyrUP autoRIFT /home/jovyan/Projects/PX_comparison/PX/autoRIF... /home/jovyan/Projects/PX_comparison/PX/autoRIF...
130 LS8-20180304-20180405 32 64 960 4 60 Gau pyrUP autoRIFT /home/jovyan/Projects/PX_comparison/PX/autoRIF... /home/jovyan/Projects/PX_comparison/PX/autoRIF...
134 LS8-20180304-20180405 32 64 960 4 60 NAOF pyrUP autoRIFT /home/jovyan/Projects/PX_comparison/PX/autoRIF... /home/jovyan/Projects/PX_comparison/PX/autoRIF...
97 LS8-20180304-20180405 32 65 975 1 15 Gau parabolic Vmap /home/jovyan/Projects/PX_comparison/PX/Vmap/pa... /home/jovyan/Projects/PX_comparison/PX/Vmap/pa...
73 LS8-20180304-20180405 32 varying: multi-pass varying: multi-pass 4.009 60.14 NAOF interest point groups GIV /home/jovyan/Projects/PX_comparison/PX/GIV/u_l... /home/jovyan/Projects/PX_comparison/PX/GIV/v_l...
28 LS8-20180304-20180405 32 64 960 1 15 NAOF 16-node oversampling CARST /home/jovyan/Projects/PX_comparison/PX/CARST/2... /home/jovyan/Projects/PX_comparison/PX/CARST/2...

This cell performs the static terrain analysis and calculates the corresponding metrics.

exps = {}

for idx, row in selected_cases.iterrows():
    exp = glaft.Velocity(vxfile=row.Vx, vyfile=row.Vy, static_area=in_shp, 
                         kde_gridsize=60, thres_sigma=2.0)
    exp.static_terrain_analysis()
    exps[idx] = exp

The following functions plot the braces on the axes with annotations, which is necessary for Figure 2. They are modified from guzey’s answer to this StackOverflow thread.

def draw_brace_x(ax, xspan: tuple=(None, None), yy: float=0.0, text: str=''):
    """
    ax: axes to be drawn.
    xspan: x coordinates of the two brace ending points.
    yy: y coordinate of the two brace ending points. 
        (The brace will be placed horizontally)
    text: annotation text.
    """

    xmin, xmax = xspan
    xspan = xmax - xmin
    ax_xmin, ax_xmax = ax.get_xlim()
    xax_span = ax_xmax - ax_xmin

    ymin, ymax = ax.get_ylim()
    yspan = ymax - ymin
    resolution = int(xspan / xax_span * 100) * 2 + 1 # sampling resolution of the sigmoid brace
    beta = 200. / xax_span                           # the higher this is, the sharper the sigmoid

    x = np.linspace(xmin, xmax, resolution)
    x_half = x[:int(resolution / 2) + 1]
    y_half_brace = (1 / (1. + np.exp(-beta * (x_half - x_half[0] )))
                  + 1 / (1. + np.exp(-beta * (x_half - x_half[-1]))))
    y = np.concatenate((y_half_brace, y_half_brace[-2::-1]))
    y = yy - (.05 * y - .01) * yspan                 # adjust vertical stretch and position

    ax.plot(x, y, color='black', lw=1)
    ax.text((xmax + xmin) / 2., yy - .07 * yspan, text, ha='left', va='top')
    
def draw_brace_y(ax, yspan: tuple=(None, None), xx: float=0.0, text: str=''):
    """
    ax: axes to be drawn.
    yspan: y coordinates of the two brace ending points.
    xx: x coordinate of the two brace ending points. 
        (The brace will be placed vertically)
    text: annotation text.
    """

    ymin, ymax = yspan
    yspan = ymax - ymin
    ax_ymin, ax_ymax = ax.get_ylim()
    yax_span = ax_ymax - ax_ymin

    xmin, xmax = ax.get_xlim()
    xspan = xmax - xmin
    resolution = int(yspan / yax_span * 100) * 2 + 1 # sampling resolution of the sigmoid brace
    beta = 200. / yax_span                           # the higher this is, the sharper the sigmoid

    y = np.linspace(ymin, ymax, resolution)
    y_half = y[:int(resolution / 2) + 1]
    x_half_brace = (1 / (1. + np.exp(-beta * (y_half - y_half[0] )))
                  + 1 / (1. + np.exp(-beta * (y_half - y_half[-1]))))
    x = np.concatenate((x_half_brace, x_half_brace[-2::-1]))
    x = xx - (.05 * x - .01) * xspan                 # adjust vertical stretch and position

    ax.plot(x, y, color='black', lw=1)
    ax.text(xx - .05 * xspan, (ymax + ymin) / 2., text, rotation=90, ha='right', va='bottom')

Now starting to make the figure:

fig = plt.figure(figsize=(13, 13))
subfigs = fig.subfigures(2, 4, wspace=0, hspace=0, width_ratios=(0.3, 0.3, 0.3, 0.1))    # last column is for colorbar
all_axs = np.empty((2,4), dtype='object')

# create two subplots in each subfigure
for i, row in enumerate(subfigs[:, :3]):
    for j, subfig in enumerate(row):
        all_axs[i, j] = subfig.subplots(2, 1, gridspec_kw = {'hspace':0, 'height_ratios':(0.3962, 0.6038)})

title_labels = np.array([['$\mathbf{a}$ \t Test #126 \n autoRIFT; No pre-filter', 
                          '$\mathbf{b}$ \t Test #130 \n autoRIFT; Gaussian HPF', 
                          '$\mathbf{c}$ \t Test #134 \n autoRIFT; NAOF'], 
                         ['$\mathbf{d}$ \t Test #97 \n Vmap; Gaussian HPF', 
                          '$\mathbf{e}$ \t Test #73 \n GIV; NAOF', 
                          '$\mathbf{f}$ \t Test #28 \n CARST; NAOF']])

for idx, i, j in [[126, 0, 0], [130, 0, 1], [134, 0, 2],
                  [97,  1, 0], [73,  1, 1], [28,  1, 2]]:
    exp = exps[idx]
    
    # top panel
    ax_sel = all_axs[i, j][0]
    cm_settings = glaft.show_velocomp(exp.vxfile, ax=ax_sel)
    ax_sel.set_aspect('equal', adjustable='datalim')
    ax_sel.set_title(title_labels[i, j])
    
    # bottom panel
    ax_sel = all_axs[i, j][1]
    exp.plot_zoomed_extent(ax=ax_sel)
    ax_sel.set_aspect('equal', adjustable='box')
    ax_sel.set_xlim(-1, 1)
    ax_sel.set_ylim(-1, 1)
    ax_sel.set_title(None)
    
    # bottom panel ticks
    ax_sel.tick_params(direction="in", bottom=True, top=True, left=True, right=True)
    ax_sel.tick_params(axis='x', pad=10)
    ax_sel.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
    ax_sel.xaxis.set_major_formatter(FormatStrFormatter('%.1f'))
    ax_sel.set_yticks([-1.0, -0.5, 0.0, 0.5, 1.0])
    ax_sel.set_xticks([-1.0, -0.5, 0.0, 0.5, 1.0])
    
    # show percentage of incorrent matches
    ax_sel.text(0.95, 0.95, '{:.1f}% incorrect matches'.format(exp.outlier_percent * 100), ha='right', va='top')
    
    # annotations of du and dv
    draw_brace_x(ax_sel, 
                 xspan=(exp.kdepeak_x - exp.metric_static_terrain_x, exp.kdepeak_x), 
                 yy=exp.kdepeak_y - exp.metric_static_terrain_y, 
                 text='$\delta_u$ = {:.2f}'.format(exp.metric_static_terrain_x))
    draw_brace_y(ax_sel, 
                 yspan=(exp.kdepeak_y - exp.metric_static_terrain_y, exp.kdepeak_y), 
                 xx=exp.kdepeak_x - exp.metric_static_terrain_x, 
                 text='$\delta_v$ = {:.2f}'.format(exp.metric_static_terrain_y))
    

# fine-tune positions of the top panels
bbox_top    = all_axs[0, 0][0].get_position()
bbox_bottom = all_axs[0, 0][1].get_position()
bbox_top.x0 = bbox_bottom.x0
bbox_top.x1 = bbox_bottom.x1
bbox_top.y0 = bbox_bottom.y1
new_bbox_top_y1 = bbox_top.y0 + (bbox_top.y1 - bbox_top.y0) * ((bbox_bottom.x1 - bbox_bottom.x0) / (bbox_top.x1 - bbox_top.x0))
bbox_top.y1 = new_bbox_top_y1

for i, row in enumerate(all_axs[:, :3]):
    for j, subfig in enumerate(row):
        all_axs[i, j][0].set_position(bbox_top)

# add colorbars
mappable = glaft.prep_colorbar_mappable(**cm_settings)
cax = subfigs[0, 3].add_axes([0.0, bbox_top.y0, 0.15, bbox_top.y1 - bbox_top.y0])
subfigs[0, 3].colorbar(mappable, cax=cax, orientation='vertical', label='$V_x$ ({})'.format(exp.velocity_unit), ticks=[-2, -1, 0, 1, 2])
cax = subfigs[1, 3].add_axes([0.0, bbox_top.y0, 0.15, bbox_top.y1 - bbox_top.y0])
subfigs[1, 3].colorbar(mappable, cax=cax, orientation='vertical', label='$V_x$ ({})'.format(exp.velocity_unit), ticks=[-2, -1, 0, 1, 2])

# add axis labels
all_axs[0, 0][1].set_yticklabels(['-1.0', '', '', '', '1.0'])                                      # to prevent label bleeding
all_axs[1, 0][1].set_yticklabels(['-1.0', '', '', '', '1.0'])                                      # ditto
all_axs[0, 0][1].set_ylabel('Static area $V_y$ ({})'.format(exp.velocity_unit), labelpad=-30.0)    # ditto
all_axs[1, 0][1].set_ylabel('Static area $V_y$ ({})'.format(exp.velocity_unit), labelpad=-30.0)    # ditto

all_axs[0, 1][1].set_xlabel('Static area $V_x$ ({})'.format(exp.velocity_unit))
all_axs[1, 1][1].set_xlabel('Static area $V_x$ ({})'.format(exp.velocity_unit))

# save figure
fig.patch.set_facecolor('xkcd:white')
fig.savefig('Fig2.png', dpi=200)
../../_images/Fig2_8_0.png