<jupyter_start><jupyter_text>Our Tile Directory Path<jupyter_code>tile_name = 'NR08'

tile_dir = Path(f'out/{tile_name}')
tile_dir_4326 = Path(f'out/{tile_name}_4326')

list(tile_dir_4326.glob('*'))<jupyter_output><empty_output><jupyter_text>Our GRWL Directory Path<jupyter_code>grwl_data = Path('grwl_data')
# Ignore hidden files
list(grwl_data.glob('[!.]*'))<jupyter_output><empty_output><jupyter_text>Open Data Rasters<jupyter_code>with rasterio.open(tile_dir_4326/f'segments_{tile_name}.tif') as ds:
    segments = ds.read(1)
with rasterio.open(tile_dir_4326/f'ocean_mask_{tile_name}.tif') as ds:
    ocean_mask = ds.read(1)
    
with rasterio.open(tile_dir_4326/f'segments_{tile_name}.tif') as ds:
    segments = ds.read(1)
    
with rasterio.open(tile_dir_4326/f'width_differences.tif') as ds:
    width_differences = ds.read(1)
    
with rasterio.open(tile_dir_4326/f'water_mask_{tile_name}.tif') as ds:
    water_mask = ds.read(1)
    profile= ds.profile

ocean_layer = ocean_mask.astype(dtype=float)
ocean_layer[~ocean_mask.astype(bool)] = np.nan

plt.imshow(water_mask, interpolation='none')<jupyter_output><empty_output><jupyter_text>Network Data<jupyter_code>%%time

edges_grwl = gpd.read_file(tile_dir_4326/f'{tile_name}.geojson')
widths_df = gpd.read_file(tile_dir_4326/'width_geometries.geojson')
edges_orinoco = gpd.read_file(tile_dir_4326/f'{tile_name}_edges.geojson')
edges_orinoco.plot()<jupyter_output>CPU times: user 5.16 s, sys: 102 ms, total: 5.27 s
Wall time: 5.53 s<jupyter_text>Inspect a subsetWe will use the below to select subsets of our ROI.<jupyter_code>sy = np.s_[500:1_250]
sx = np.s_[2300:2_950]
plt.imshow(water_mask[sy, sx])<jupyter_output><empty_output><jupyter_text>Automating the plotsWe selected two slices using the `imshow` above. We have two slices indices encoded via `if ... then...` below. We create all the relevant plots over the subsets in what follows. We have encoded two subsets via `0` and `1`.<jupyter_code># Options are 0 and 1 for different areas over the Mackenzie.
slice_index = 1
# Mackenzie
figsize=(12, 10)

if (tile_name is 'NR08') & (slice_index == 0):
    sy = np.s_[2_000:3_000]
    sx = np.s_[4_000:5_000]
    loc='upper left'

elif (tile_name is 'NR08') & (slice_index == 1):
    sy = np.s_[5_600:6_600]
    sx = np.s_[7_750:8_750]
    loc='lower right'<jupyter_output><empty_output><jupyter_text>Making a Plot DirectoryWe make a separate directory for each subset (i.e. slice index) we have selected.<jupyter_code>plot_dir = Path(f'plot/{tile_name}_{slice_index}')
plot_dir.mkdir(exist_ok=True, parents=True)<jupyter_output><empty_output><jupyter_text>Getting the Geographic Metadata<jupyter_code>profile_cropped = get_cropped_profile(profile, sx, sy)
transform_c = profile_cropped['transform']
x_min, y_min, x_max, y_max = array_bounds(profile_cropped['height'], 
                                          profile_cropped['width'], 
                                          profile_cropped['transform'])<jupyter_output><empty_output><jupyter_text>Corpping our Proudcts to the Area<jupyter_code>edges_grwl_c = edges_grwl.cx[x_min: x_max, y_min: y_max].copy()
edges_orinoco_c = edges_orinoco.cx[x_min: x_max, y_min: y_max].copy()
widths_df_c = widths_df.cx[x_min: x_max, y_min: y_max].copy()

ocean_layer_c = ocean_layer[sy, sx]
water_mask_c = water_mask[sy, sx]
M = 4
labels = ['Land', 'Inland Water', 'Ocean', 'Orinoco']

cNorm  = colors.Normalize(vmin=0, vmax=(M))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap='viridis')

legend_elements = [Line2D([0], 
                          [0], 
                          marker='s',
                          color='w',
                          label=f'{labels[k]}',
                          markerfacecolor=scalarMap.to_rgba(k+1),
                          markeredgecolor='black',
                          markersize=10) for k in range(M-2)] +\
                  [Line2D([0], 
                          [0], 
                          marker=None,
                          color=scalarMap.to_rgba(M),
                          label=f'{labels[M-1]}',
                          linewidth=5,
                          #markerfacecolor=scalarMap.to_rgba(k),
                          #markeredgecolor='black',
                          #edgecolor='black',
                          markersize=10)]
fig, ax = plt.subplots(figsize=(12, 10))
plot.show(water_mask_c, 
          ax=ax, 
          transform=transform_c, 
          zorder=1,  
          vmin=0, 
          vmax=2, 
          interpolation='none',
          alpha=1)

plot.show(ocean_layer_c, 
          ax=ax, 
          transform=transform_c, 
          zorder=1,  
          vmin=0, 
          vmax=2, 
          interpolation='none',
          alpha=1)

ax.set_xlim((x_min, x_max))
ax.set_ylim((y_min, y_max))

edges_orinoco_c.plot(edgecolor=scalarMap.to_rgba(4),
                     linewidth=2,
                     ax=ax,
                     zorder = 3) 
                  

ax.legend(handles=legend_elements, loc=loc, fontsize=22, framealpha=1)
ax.grid('on')

ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

plt.savefig(plot_dir/'orinoco.png')
M = 4
labels = ['Land', 'Inland Water', 'Ocean', 'GRWL']

legend_elements = [Line2D([0], 
                          [0], 
                          marker='s',
                          color='w',
                          label=f'{labels[k]}',
                          markerfacecolor=scalarMap.to_rgba(k+1),
                          markeredgecolor='black',
                          markersize=10) for k in range(M-2)] +\
                  [Line2D([0], 
                          [0], 
                          marker=None,
                          color=scalarMap.to_rgba(M),
                          label=f'{labels[M-1]}',
                          linewidth=5,
                          #markerfacecolor=scalarMap.to_rgba(k),
                          #markeredgecolor='black',
                          #edgecolor='black',
                          markersize=10)]
fig, ax = plt.subplots(figsize=(12, 10))
plot.show(water_mask_c, 
          ax=ax, 
          transform=transform_c, 
          zorder=1,  
          vmin=0, 
          vmax=2, 
          interpolation='none',
          alpha=1)

plot.show(ocean_layer_c, 
          ax=ax, 
          transform=transform_c, 
          zorder=1,  
          vmin=0, 
          vmax=2, 
          interpolation='none',
          alpha=1)

ax.set_xlim((x_min, x_max))
ax.set_ylim((y_min, y_max))

edges_grwl_c.plot(edgecolor=scalarMap.to_rgba(4),
                  linewidth=2,
                  ax=ax,
                  zorder = 3) 
                  
ax.legend(handles=legend_elements, loc=loc, fontsize=22, framealpha=1)
ax.grid('on')

# Axes

ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

plt.savefig(plot_dir/'grwl.png')

M = 1

cNorm  = colors.Normalize(vmin=0, vmax=(M))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap='viridis')

legend_elements = [Line2D([0], 
                          [0], 
                          marker='s',
                          color='w',
                          label=f'No data',
                          markerfacecolor=scalarMap.to_rgba(1),
                          markeredgecolor='black',
                          markersize=10) ] +\
                  [Line2D([0], 
                          [0], 
                          marker=None,
                          color='green',
                          label=f'Width',
                          linewidth=2,
                          markersize=10)]
import scipy.ndimage as nd
mask_crop = water_mask[sy, sx]
outline_c = nd.morphology.binary_dilation(mask_crop, iterations = 3).astype(float)
outline_c[mask_crop.astype(bool)] = 0.
outline_c[outline_c == 0] = np.nan

nodata_diff = (water_mask & np.isnan(width_differences)).astype(float)
nodata_diff[nodata_diff == 0] = np.nan
nodata_diff_c = nodata_diff[sy, sx]
width_differences_c = width_differences[sy, sx]
cbounds = np.max(np.abs(width_differences_c[~np.isnan(width_differences_c)]))

fig, ax = plt.subplots(figsize=(12, 10))
cmap='seismic'
plot.show(width_differences_c, 
          ax=ax, 
          transform=transform_c, 
          zorder=1,  
          cmap=cmap,
          vmin=-cbounds,
          vmax=cbounds,
          interpolation='none')

plot.show(nodata_diff_c, 
          ax=ax, 
          transform=transform_c, 
          zorder=1,  
          vmin=0, 
          vmax=1,
          interpolation='none')

plot.show(outline_c, 
          ax=ax, 
          transform=transform_c,
          interpolation='none')

ax.set_xlim((x_min, x_max))
ax.set_ylim((y_min, y_max))

    
widths_df_c.plot(color='green', ax=ax, alpha=1,zorder=4)

##########
# Colorbar
##########

# Source: https://matplotlib.org/3.1.0/tutorials/colors/colorbar_only.html

divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='3.5%', pad=0.07)

norm = colors.Normalize(vmin=-cbounds, vmax=cbounds)
cb1 = colorbar.ColorbarBase(cax, 
                            cmap=plt.get_cmap(cmap),
                            norm=norm,
                            orientation='vertical')
cax.set_ylabel('Orinoco - GRWL (m)', fontsize=30)
cax.tick_params(axis='both', which='major', labelsize=20)

# Legends

ax.legend(handles=legend_elements, loc=loc, fontsize=22, framealpha=1)
ax.grid('on')


# Axes
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

plt.savefig(plot_dir/'width_differences.png')
fig, ax = plt.subplots(figsize=(12, 10))
X = width_differences[sy, sx]
labels = segments[sy, sx]
data = get_features_from_array(labels, X)
data = data[~np.isnan(data)]
ax.hist(data, bins=50)
ax.set_xlabel('Orinoco - GRWL (m)')
ax.set_ylabel('\# of Segments')
ax.text(.6, .8, f'Total Segments = {len(data):,}', transform=ax.transAxes)
plt.savefig(plot_dir/'diff_histogram.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13099949264332825, 'fraction_numerical': 0.02475900558092339, 'mean_word_length': 2.335363790186125, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '516408', 'license': 'Apache-2.0', 'n_tokens_mistral': 3673, 'n_tokens_neox': 3450, 'n_words': 528, 'path': 'notebooks/validation_with_grwl/4 - Figures - Connectivity.ipynb', 'repo_name': 'simard-landscape-lab/orinoco', 'size': '9855'}