<jupyter_start><jupyter_text>We use `notebooks/examples/1 - Make Channel Network` removing extraneous discussion and visualizations wherever possible. We also include a node attribute that uses the distance transform in a node's associated segment. This will be used to very roughly compare [RivGraph](https://github.com/jonschwenk/RivGraph/tree/master/rivgraph) and [ChanGeom](https://www.burchfisher.com/data.html)'s approach to measuring width.<jupyter_code>import rasterio
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
import networkx as nx
from orinoco import (get_distance_in_channel,
                     get_distance_segments,
                     get_undirected_channel_network,
                     direct_channel_network,
                     export_edges_to_geodataframe,
                     export_nodes_to_geodataframe,
                     get_map_centroid_from_binary_mask,
                     add_flow_attributes,
                     split_tuple_pairs,
                     get_segment_df,
                     get_geo_width_df,
                     update_graph_with_widths,
                     get_width_features_from_segments,
                     add_width_features_to_graph,
                     get_array_from_features
                     )
from skimage.color import label2rgb
import random
import geopandas as gpd
from shapely.geometry import Point
from tqdm import tqdm<jupyter_output><empty_output><jupyter_text>Initialize the Tile Directory<jupyter_code>tile_name = 'NH15'

tile_dir = Path(f'out/{tile_name}')
tile_dir.exists()<jupyter_output><empty_output><jupyter_text>Read Data from the Tile Directory<jupyter_code>with rasterio.open(tile_dir/f'water_mask_{tile_name}.tif') as ds:
    water_mask = ds.read(1)
with rasterio.open(tile_dir/f'ocean_mask_{tile_name}.tif') as ds:
    ocean_mask = ds.read(1)
    profile = ds.profile
plt.imshow(water_mask, interpolation='none')<jupyter_output><empty_output><jupyter_text>Generate Orinoco Products Obtaining the Distance Function<jupyter_code>transform = profile['transform']
dx, dy = transform.a, -transform.e
dist = get_distance_in_channel(water_mask,
                               ocean_mask,
                               dx=dx,
                               dy=dy,
                               # removes areas with less than percentage of total size
                               min_rel_area=0.0,
                               # Applys a 1 pixel buffer for distance computation
                               # for artifical 8-connectivity.
                               # See: https://github.com/scikit-fmm/scikit-fmm/issues/32
                               # Removes buffer after computation complete.
                               # Only recommended for GRWL
                               apply_mask_buffer=True
                              )
plt.imshow(dist, interpolation='none')<jupyter_output><empty_output><jupyter_text>Save the distance product.<jupyter_code>p = profile.copy()
p['dtype'] = 'float32'
with rasterio.open(tile_dir/f'distance_{tile_name}.tif', 'w', **p) as ds:
    ds.write(dist.astype(np.float32), 1)<jupyter_output><empty_output><jupyter_text>Segmentation of the Channel<jupyter_code>pixel_step = 5
segments, interface_adj_segments = get_distance_segments(dist,
                                                         pixel_step,
                                                         dx=dx,
                                                         dy=dy,
                                                         connectivity=8,
                                                         min_size=4)
segments_rgb = label2rgb(segments, bg_label=0)
plt.imshow(segments_rgb, interpolation='none')<jupyter_output><empty_output><jupyter_text>We save both the `segments` and the `segments_rgb` to disk.<jupyter_code>p = profile.copy()
p['dtype'] = 'int32'
with rasterio.open(tile_dir/f'segments_{tile_name}.tif', 'w', **p) as ds:
    ds.write(segments.astype(np.int32), 1)

p = profile.copy()
p['dtype'] = 'float32'
p['count'] = 3
with rasterio.open(tile_dir/f'segments_rgb_{tile_name}.tif', 'w', **p) as ds:
    ds.write(segments_rgb.transpose([2, 0, 1]).astype(np.float32))<jupyter_output><empty_output><jupyter_text>Obtaining the Undirected Network<jupyter_code>chanG_undirected = get_undirected_channel_network(segments, 
                                                  dist, 
                                                  profile, 
                                                  interface_adj_segments,
                                                  connectivity=8)
node_data = chanG_undirected.nodes(data=True)
[(node, data) for (node, data) in node_data if data['label'] == 268]<jupyter_output><empty_output><jupyter_text>Obtaining the Directed Network<jupyter_code>chanG = direct_channel_network(chanG_undirected, 
                               # The keywords below are how the pruning occurs
                               # Specifies pruning will be done
                               remove_dangling=True, 
                               # Do not prune nodes within 1 km of interface
                               interface_buffer_m=1_000,     
                               # Remove edge groups with an degree 1 endpoint and size <=3
                               group_min_size=3,
                               # How many times to do this pruning
                               dangling_iterations=1
                               )<jupyter_output><empty_output><jupyter_text>Let's do a quick sanity check.<jupyter_code>df_edges = export_edges_to_geodataframe(chanG, profile['crs'])
df_edges.plot(column='segment_id', categorical=True)<jupyter_output><empty_output><jupyter_text>Obtaining the Widths Using the distance transformWe are going to use the distance transform within a segment. We have a function that determines the width according to $2\cdot d - 1$ where d is the distance to the land mask according to the distance transform in a small 1-pixel buffer around the segment (if this is 0, we use the full distance transform on the channel mask).<jupyter_code>width_features = get_width_features_from_segments(segments, profile)
widths = get_array_from_features(segments, width_features)
plt.imshow(widths, interpolation='none')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>We now add these attributes to the graph.<jupyter_code>chanG = add_width_features_to_graph(chanG, width_features.ravel())
node_data =dict(chanG.nodes(data=True))
random.choice(list(node_data.items()))<jupyter_output><empty_output><jupyter_text>Using $\nabla \varphi$We compute width using the gradient of the distance function (and possibly the network structure). First we obtain the orientation of along wich to obtain these measured widths (perpendicular to the artificial flow).<jupyter_code>chanG = add_flow_attributes(chanG, dist, profile['transform'])<jupyter_output>flow vector computation using gradient: 100%|██████████| 31377/31377 [00:00<00:00, 63232.66it/s]
flow vector computation using network: 100%|██████████| 31377/31377 [00:02<00:00, 11006.29it/s]<jupyter_text>To perform intersection, we must extract the polygons of the segments. This takes a bit of time to polygonize the segments, so be patient.<jupyter_code>%%time

df_segments = get_segment_df(segments, chanG, profile, connectivity=8)
df_segments.head()<jupyter_output>CPU times: user 30.5 s, sys: 1.31 s, total: 31.8 s
Wall time: 31.7 s<jupyter_text>We will save the geometries to a file for later inspection.<jupyter_code>df_segments_out = split_tuple_pairs(df_segments)
df_segments_out.to_file(tile_dir/'segments.geojson', 
                        driver='GeoJSON')<jupyter_output>/Users/cmarshak/anaconda3/lib/python3.7/site-packages/geopandas/io/file.py:108: FionaDeprecationWarning: Use fiona.Env() instead.
  with fiona.drivers():<jupyter_text>Now we extract the widths are associated geometry. This is the most expensive computation in the notebook, so be very patient. Could take 5 minutes or so.<jupyter_code>%%time

df_geo_widths = get_geo_width_df(df_segments, chanG, radius=1)
df_geo_widths.head()<jupyter_output>CPU times: user 2min 27s, sys: 515 ms, total: 2min 27s
Wall time: 2min 27s<jupyter_text>Let's save the widths again for later inspection.<jupyter_code>df_geo_widths_out = split_tuple_pairs(df_geo_widths)
df_geo_widths_out = df_geo_widths_out[~df_geo_widths_out.geometry.is_empty]
df_geo_widths_out.to_file(tile_dir/'width_geometries.geojson', 
                          driver='GeoJSON')<jupyter_output><empty_output><jupyter_text>We update the node attributes. And we see all of the attributes we have added.<jupyter_code>len(chanG)
chanG = update_graph_with_widths(chanG, df_geo_widths)

node_data =dict(chanG.nodes(data=True))
random.choice(list(node_data.items()))<jupyter_output><empty_output><jupyter_text>Save NetworkWe now save the nodes and edges to `geojson` files and serialize the NetworkX `DiGraph` with `pickle`.<jupyter_code>df_nodes = export_nodes_to_geodataframe(chanG, profile['crs'])
df_nodes = split_tuple_pairs(df_nodes)
df_edges = export_edges_to_geodataframe(chanG, profile['crs'])
df_edges.to_file(tile_dir/f'{tile_name}_edges.geojson', driver='GeoJSON')
df_nodes.drop(columns=['interface_adj']).to_file(tile_dir/f'{tile_name}_nodes.geojson', driver='GeoJSON')
nx.write_gpickle(chanG, tile_dir/f'{tile_name}_network.pkl')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08760977674320178, 'fraction_numerical': 0.014919056184530738, 'mean_word_length': 2.9285120532003326, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 73, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '30002', 'license': 'Apache-2.0', 'n_tokens_mistral': 2979, 'n_tokens_neox': 2812, 'n_words': 627, 'path': 'notebooks/validation_with_usgs/2 - Orinoco Products.ipynb', 'repo_name': 'simard-landscape-lab/orinoco', 'size': '9451'}