veinsoftheearth / rivgraph Goto Github PK

RivGraph identifies parallel links for its directionality finding, but it doesn't consider them as parallel links in the delta metrics (adjacency matrix). In the weighted (by width) adjacency matrix it currently gives the weight as the width of one of the links. This may or may not be a problem. An example on a very simple channel network is here.

A prior solution #10 used to introduce additional nodes, although this may not be totally necessary. It depends on the metrics being computed. To have this functionality, for RG object named delt with the initial network skeletonized (but not pruned) use the following (from @jonschwenk) order of operations.

delt.prune_network(path_shoreline=path_shrln, path_inletnodes=path_inlet)
delt.compute_link_width_and_length()
delt.assign_flow_directions()
delt.links, delt.nodes = lnu.add_artificial_nodes(delt.links, delt.nodes, delt.gdobj)
delt.compute_link_width_and_length() # Need to re-compute because we've split some links into 2 in order to add artificial nodes.
adj_w = delt.adjacency_matrix(weight='wid_adj') # Can also weight by 'len_adj' or provide a vector of your own weights.

Hi,
I went through the documentation of RivGraph and I am having some difficulty to find clear instructions on how to calculate the migration rates in multi threaded ( especially braided and anabraching) rivers from binary channel masks, using RivGraph.
Could someone please point me to the right resource?

Thank you

Hi.

I am not being able to import the RivGraph in Google Colab. I tried different ways of doing it, and I failed in all.

This is my code: https://colab.research.google.com/drive/1FoW42jdMvues7-BL0CyqQrdjQTXHauFd?usp=sharing

  # Install Anaconda
  !wget -c https://repo.continuum.io/archive/Anaconda3-5.1.0-Linux-x86_64.sh
  !chmod +x Anaconda3-5.1.0-Linux-x86_64.sh
  !bash ./Anaconda3-5.1.0-Linux-x86_64.sh -b -f -p /usr/local
  !conda install -y --prefix /usr/local -c <<<your wish>>>>
  
  !git clone https://github.com/jonschwenk/RivGraph
  !conda env create --file /content/RivGraph/environment.yml
  !conda activate rivgraph

Is there a better way of doing it?

Thank you in advance.

This is an issue I've come across and solved, but good to know for debugging/figuring out why networks don't come out. Not sure of a smarter way to handle it automatically.

Sometimes a cyclical channel skeleton is extracted that wraps around the edge of the image. This is primarily caused when using a watermask that has no data values on the edges. This arises typically when projecting rectangles from WGS84 (EPSG 4326) into a local UTM zone, e.g. with the global surface water masks. Here's an example of the resulting skeleton (tealish is no data, blue is water, grey land, yellow is the shoreline, orange is skeleton even after pruning).

Solution
Define the subaerial extent of the delta as a line vector which'll be used as the argument of the shoreline in prune_network. This will clip out the weird skeleton outside the image. When doing this, the inlet node passed to prune network should have its nearest neighbor within the subaerial delta to work. Example here:

There seems to be a gap between what the various geospatial packages accept as valid CRS information.

For example, the below CRS information is read by RivGraph using GDAL, and can works with the GDAL-based functions like to_geotiff(), but does not work with the geopandas based functions like to_geovectors().

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- undefined
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

A suitable fix for this specific case was to manually set the crs information stored in the delta class to delta.crs = "EPSG:4326" as that is recognized and accepted by geopandas.

Need to do some poking around to figure this out as I doubt it is an uncommon problem... other geospatial projects that depend on GDAL and geopandas must have some sort of solution to handle this mismatch.

EDIT: This bug might just be a weird line in RivGraph somewhere, seems like geopandas is designed to work with pyproj CRS information...

Hello!
I have been trying to install with the different methods proposed in the website. However, it is forever processing without installing. Is it necessary a specific Conda version? I could not find any information about it. I also tried with a fresh installation of Anaconda and Miniconda, but it did not work 😭.
Thanks in advance for your help.

On the first page of the documentation, there is a broken link to the FAQ:

Hi! Any ideas why I am getting this error? Also getting a similar error when running the compute.mesh() function.
Seems to be related to the shapely package (I have the 1.7.1 version).

Here is a screenshot of the full sequence of errors.

Thanks!
Tanya

I want to run some demo of rivgraph in my jupyter server,however when I try to launch the code, it crashed like this:

Traceback (most recent call last)
Input In [2], in <cell line: 1>()
----> 1 from rivgraph.classes import river
2 import matplotlib.pyplot as plt
3 import os

ModuleNotFoundError: No module named 'rivgraph'

And before this I have installed rivgraph's conda environment:

(rivgraph) jupyterhub@vm-jupyterhub-server:~/graphs$ conda env list
rivgraph * /home/jupyterhub/.conda/envs/rivgraph
base /opt/conda
db-base /opt/conda/envs/db-base
gis-base /opt/conda/envs/gis-base
lf /opt/conda/envs/lf
ml-base /opt/conda/envs/ml-base
r-base /opt/conda/envs/r-base

So how to solve this problem?

Hello!

I'm trying to use Rivgraph to process the results of a braided river simulation (similar to what you did in the Tejedor (2022) article). I'm following a combination of the steps shown in the three available examples, and so far I've the network shown below (I manually disabled the labels for nodes and links).

The problem I am facing is that I get an error when I use the compute_mesh() function. It seems to be a problem with a distance calculation, but I can't fully understand it. I would really appreciate any feedback or recommendations on how to proceed to get it working.

If it's helpful, here's the error detail:

Traceback (most recent call last):

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/spyder_kernels/py3compat.py", line 356, in compat_exec
    exec(code, globals, locals)

  File "/mnt/data1/GITHUB/Iber4Rivgraph/main.py", line 30, in <module>
    links, nodes = f.getNetwork(RastersPath)

  File "/mnt/data1/GITHUB/Iber4Rivgraph/functions.py", line 173, in getNetwork
    braidedRiver.compute_mesh()

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/rivgraph/classes.py", line 856, in compute_mesh
    self.meshlines, self.meshpolys, self.centerline_smooth = ru.valleyline_mesh(self.centerline, self.avg_chan_width, buf_halfwidth, grid_spacing, smoothing=smoothing)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/rivgraph/rivers/river_utils.py", line 769, in valleyline_mesh
    llines, lmap = iterative_cl_pt_mapping(cl2, bdists, 'left')

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/rivgraph/rivers/river_utils.py", line 608, in iterative_cl_pt_mapping
    distance, path = fastdtw(Ao, An, dist=euclidean)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 53, in fastdtw
    return __fastdtw(x, y, radius, dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 73, in __fastdtw
    __fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 73, in __fastdtw
    __fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 73, in __fastdtw
    __fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 73, in __fastdtw
    __fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 68, in __fastdtw
    return dtw(x, y, dist=dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 130, in dtw
    return __dtw(x, y, None, dist)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/fastdtw/fastdtw.py", line 141, in __dtw
    dt = dist(x[i-1], y[j-1])

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/scipy/spatial/distance.py", line 518, in euclidean
    return minkowski(u, v, p=2, w=w)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/scipy/spatial/distance.py", line 461, in minkowski
    u = _validate_vector(u)

  File "/home/epfl-lhe/mambaforge/envs/Iber4Rivgraph/lib/python3.9/site-packages/scipy/spatial/distance.py", line 301, in _validate_vector
    raise ValueError("Input vector should be 1-D.")

ValueError: Input vector should be 1-D.

Had an issue with the following tidal inlet where remove_all_spurs in prune_network removes the entire network on the first run. Relevant data attached here.
Note that I have a slightly modified version of remove_all_spurs, which will keep a single link channel network (used for some boring river mouths). Only relevant because it ensures the entire network isn't pruned, but keeps the link attached to the inlet. See the modification below

        # Remove all the nodes with only two links attached
        links, nodes = remove_two_link_nodes(links, nodes, dontremove)
 
        # For single channel (no loops & single outlet) networks,
        # prevent overpruning beyond a single network with 2 nodes.
        if (links['n_networks'] == 1) & (len(nodes['conn']) == 2):
            ct = 0
            
        if ct == 0:
            stopflag = 1

    return links, nodes

Hello. I am following the instructions for building from source.

I believe these instructions should mention that the following command is necessary before step 2: pip install -r requirements.txt. If so, then the requirements.txt file should be updated as it currently does not work:

$ pip install -r requirements.txt 
Collecting numpy
  Using cached numpy-1.19.5-cp38-cp38-manylinux2010_x86_64.whl (14.9 MB)
Collecting scikit-image
  Using cached scikit_image-0.18.1-cp38-cp38-manylinux1_x86_64.whl (30.2 MB)
ERROR: Could not find a version that satisfies the requirement opencv (from -r requirements.txt (line 3)) (from versions: none)
ERROR: No matching distribution found for opencv (from -r requirements.txt (line 3))
(rivgraph) leo@syenite: RivGraph [master]$ pip install opencv
ERROR: Could not find a version that satisfies the requirement opencv (from versions: none)
ERROR: No matching distribution found for opencv

Currently, requirements.txt lists opencv. I believe the correct dependency is opencv-python. After changing opencv for opencv-python, the command pip install -r requirements.txt runs smoothly. (Though it should perhaps be made to match environment.yml.)

After that, I tried the following:

$ python setup.py install
Traceback (most recent call last):
  File "setup.py", line 2, in <module>
    from rivgraph import __version__
  File "/home/leo/code/RivGraph/rivgraph/__init__.py", line 3, in <module>
    from rivgraph.classes import delta, river
  File "/home/leo/code/RivGraph/rivgraph/classes.py", line 9, in <module>
    import gdal
ModuleNotFoundError: No module named 'gdal'

I tried running apt install libgdal-dev, which executed successfully. However, the command pip install gdal now fails.

As a result, the command python setup.py install fails with error ModuleNotFoundError: No module named 'gdal'.

As a sidenote, I tried the installation with conda, and it also failed. This may have nothing to do with your package as I have never used conda before. Note I would much rather be able to install from source than having to download and use conda. Requiring conda adds 500MB+ of space required.

For completeness, after downloading, installing, and verifying my installation of conda, I got the following error:

$ conda env create --file environment.yml

# >>>>>>>>>>>>>>>>>>>>>> ERROR REPORT <<<<<<<<<<<<<<<<<<<<<<

    Traceback (most recent call last):
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda/exceptions.py", line 1079, in __call__
        return func(*args, **kwargs)
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda_env/cli/main.py", line 80, in do_call
        exit_code = getattr(module, func_name)(args, parser)
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda_env/cli/main_create.py", line 87, in execute
        spec = specs.detect(name=name, filename=filename, directory=os.getcwd())
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda_env/specs/__init__.py", line 43, in detect
        if spec.can_handle():
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda_env/specs/yaml_file.py", line 18, in can_handle
        self._environment = env.from_file(self.filename)
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda_env/env.py", line 160, in from_file
        return from_yaml(yamlstr, filename=filename)
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda_env/env.py", line 141, in from_yaml
        data = yaml_safe_load(yamlstr)
      File "/home/leo/anaconda3/lib/python3.8/site-packages/conda/common/serialize.py", line 67, in yaml_safe_load
        return yaml.safe_load(string, version="1.2")
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/main.py", line 980, in safe_load
        return load(stream, SafeLoader, version)
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/main.py", line 935, in load
        return loader._constructor.get_single_data()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/constructor.py", line 109, in get_single_data
        node = self.composer.get_single_node()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/composer.py", line 78, in get_single_node
        document = self.compose_document()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/composer.py", line 104, in compose_document
        self.parser.get_event()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/parser.py", line 163, in get_event
        self.current_event = self.state()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/parser.py", line 239, in parse_document_end
        token = self.scanner.peek_token()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/scanner.py", line 182, in peek_token
        self.fetch_more_tokens()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/scanner.py", line 282, in fetch_more_tokens
        return self.fetch_value()
      File "/home/leo/anaconda3/lib/python3.8/site-packages/ruamel_yaml/scanner.py", line 651, in fetch_value
        raise ScannerError(
    ruamel_yaml.scanner.ScannerError: mapping values are not allowed here
      in "<unicode string>", line 90, column 65:
         ... e" content="{&quot;version&quot;: &quot;4&quot;, &quot;rollouts& ... 
                                             ^ (line: 90)

`$ /home/leo/anaconda3/bin/conda-env create --file environment.yml`

  environment variables:
                 CIO_TEST=<not set>
  CONDA_AUTO_UPDATE_CONDA=false
        CONDA_DEFAULT_ENV=base
                CONDA_EXE=/home/leo/anaconda3/bin/conda
             CONDA_PREFIX=/home/leo/anaconda3
    CONDA_PROMPT_MODIFIER=(base)
         CONDA_PYTHON_EXE=/home/leo/anaconda3/bin/python
               CONDA_ROOT=/home/leo/anaconda3
              CONDA_SHLVL=1
           CURL_CA_BUNDLE=<not set>
                     PATH=/home/leo/anaconda3/bin:/home/leo/anaconda3/bin:/home/leo/anaconda3/co
                          ndabin:/home/leo/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin
                          :/sbin:/bin:/usr/games:/usr/local/games:/snap/bin
       REQUESTS_CA_BUNDLE=<not set>
            SSL_CERT_FILE=<not set>
               WINDOWPATH=2

     active environment : base
    active env location : /home/leo/anaconda3
            shell level : 1
       user config file : /home/leo/.condarc
 populated config files : /home/leo/.condarc
          conda version : 4.9.2
    conda-build version : 3.20.5
         python version : 3.8.5.final.0
       virtual packages : __glibc=2.31=0
                          __unix=0=0
                          __archspec=1=x86_64
       base environment : /home/leo/anaconda3  (writable)
           channel URLs : https://repo.anaconda.com/pkgs/main/linux-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/linux-64
                          https://repo.anaconda.com/pkgs/r/noarch
          package cache : /home/leo/anaconda3/pkgs
                          /home/leo/.conda/pkgs
       envs directories : /home/leo/anaconda3/envs
                          /home/leo/.conda/envs
               platform : linux-64
             user-agent : conda/4.9.2 requests/2.24.0 CPython/3.8.5 Linux/5.6.0-1042-oem ubuntu/20.04.1 glibc/2.31
                UID:GID : 1001:1001
             netrc file : None
           offline mode : False


An unexpected error has occurred. Conda has prepared the above report.

If a channel network mask is provided for which the outlet links are not connected via an ocean or other waterbody, the pruning performed by prune_network() either removes too much of the network or doesn't prune anything.

Example:
Binary mask

Network:

Pruned network (and shoreline):

No actual pruning was performed in this case.

This likely arises from the initial development of pruning algorithms which assumed that outlet links would be connected by a waterbody--typically an ocean or lake. This assumption needs to be corrected.

I have noticed some extra nodes remain especially around islands. This leads to an extra link around those islands. Not sure what may be causing this (maybe spurs near those islands as well). I will look into it more later and if I find anything will post here as well. These should have been removed during the remove spurs call. Data (shoreline, inlet node, nodes/links jsons, skel) is below the examples:

Data:
Yukon_GSW_201407.zip

Ran into the following bug today on Delft3D data where the delta metrics failed due to "multiple inlets" but no sources/sinks were found in the directionality step. It was because there was a node disconnected from the rest of the network in the adjacency matrix due to a mask assumption. I'm attaching the data here as an example and will refer explicitly to the node/link ids here.

Here's the the code snippet below showing what happened. The explanation below details what caused this bug and what needs to be done.

path_mask = delta_dir + "/projected/network_" + str(i) + ".tif"

# Instantiate the delta class
delt = delta(name, path_mask, path_results, verbose=True)
delt.compute_network()
delt.to_geovectors('network', ftype = 'json')


path_shrln = os.path.join(delta_dir, "/shoreline/shoreline_" + str(i) + ".shp")
path_inlet = os.path.join(delta_dir, "inlet_nodes.shp")

path_shrln = delta_dir + "/shoreline/shoreline_" + str(i) + ".json"
path_inlet = delta_dir + "/inlet_nodes.shp"

delt.prune_network(path_shoreline=path_shrln, path_inletnodes=path_inlet)

delt.compute_link_width_and_length()
delt.assign_flow_directions()
# delt.links, delt.nodes = lnu.add_artificial_nodes(delt.links, delt.nodes, delt.gdobj)
delt.compute_link_width_and_length() # Need to re-compute because we've split some links into 2 in order to add artificial nodes.
delt.to_geovectors('network', ftype='json')

delt.plot('directions')

delt.to_geotiff('directions')

adj_w = delt.adjacency_matrix(weight='wid_adj') # Can also weight by 'len_adj' or provide a vector of your own weights.
# print(adj_w)

path_adj = path_results + "/w_adj.csv"

np.savetxt(path_adj, adj_w, delimiter=",")
delt.save_network()

nodes = delt.nodes
links = delt.links

delt.compute_topologic_metrics()
Resolving links and nodes...done.
Computing distance transform...done.
Computing link widths and lengths...done.
Using ~\W0.01T0.01_rep_fixlinks.csv to manually set flow directions.
The following cycles (links) were fixed, but should be manually checked: [[352, 328, 329, 339, 342, 349, 351]]
Computing link widths and lengths...done.
Geotiff written to ~\W0.01T0.01_rep_link_directions.tif.
Links and nodes saved to pickle file: ~\W0.01T0.01_rep_network.pkl.
Traceback (most recent call last):

  File "<ipython-input-153-0cda936bec9b>", line 39, in <module>
    delt.compute_topologic_metrics()

  File "C:\Anaconda3\envs\RG_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg\rivgraph\classes.py", line 670, in compute_topologic_metrics
    self.topo_metrics = dm.compute_delta_metrics(self.links, self.nodes)

  File "C:\Anaconda3\envs\RG_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg\rivgraph\deltas\delta_metrics.py", line 40, in compute_delta_metrics
    deltavars = intermediate_vars(G)

  File "C:\Anaconda3\envs\RG_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg\rivgraph\deltas\delta_metrics.py", line 124, in intermediate_vars
    deltavars['apex'], deltavars['outlets'] = find_inlet_outlet_nodes(deltavars['A_w'])

  File "C:\Anaconda3\envs\RG_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg\rivgraph\deltas\delta_metrics.py", line 221, in find_inlet_outlet_nodes
    raise RuntimeError('The graph contains more than one apex.')

RuntimeError: The graph contains more than one apex.

The inlet node id is 257 and is in the 90th row/col of the adjacency matrix (90th from zero-index). Node 176 (row/col 47) however has no incoming values in the adjacency matrix. It does have 3 links connected to it in the node/link representation, but there's one link (235) crossing a small (3 pixel) island that doesn't click into node for whatever reason. Not sure if this is because the graphiphy function (or the underlying functions) assume the directionality pixels are all water.

The directionality pixel was land because in the skeletonize_graph call in m2g, any holes (i.e. islands) that are 4 pixels or less are filled, that is no islands of 4 pixels or less are allowed in RG. This is a fine assumption, but (1) is it totally necessary? and (2) I was checking the docs and didn't see anything about it is posted!

This bug also resulted in the nodes['inlet'] array becoming empty for some reason, not sure why.

The mask issue can be fixed by putting something into the documents, or removing this assumption, what do you folks think is better? I can push something into the docs for the former.

Additionally should some kind of error catcher be put into the dm.graphiphy to report if there's multiple inlet nodes? This would probably let the user know something is wrong. I can push something for this as well.

Hi. I keep getting error in valley_mesh creation like this:

line 806, in valleyline_mesh
if cs[check_pts[0]] - cs[keep_pts[-1]] > bufferdist:
IndexError: index -20744 is out of bounds for axis 0 with size 13160.

I've tried changing smoothing and buffer parameters, but no luck. I'll try to take a deeper look, but in the mean time any help is appreciated.

Thank you in advance!

The call to import io_utils should have rivgraph as a prefix (line 15). Fixed code below:

import rivgraph.io_utils as io

I just replaced the use of OpenCV's contour finding in rabpro's copy of regionprops. It's not well-tested at this point but I think the use cases should all be simple enough--only thing I'm not sure about are holes in the blobs.

Anyway, assuming it holds, we can replace cv2 functions in RG's regionprops function using the same code as the commit referenced here VeinsOfTheEarth/rabpro#50 points to. There are still some convolutions done by OpenCV in RivGraph that would need to be replaced, but I think there should be drop-ins (or very close) for those in skimage.

When running prune_network in the default configuration, find_inlet_nodes from delta_utils produces an error. It arises from:

inlets_epsg = int(inlets_gpd.crs['init'].split(':')[1])  
#this returns TypeError: 'CRS' object is not subscriptable

I use an inlet nodes file with projection EPSG 32603 - the same one I used in a much older version. I replaced the line with this code below which worked, outputted 32603 and the section ran. I was also given a warning that 'init' is being deprecated, may be a result of that.

inlets_epsg = inlets_gpd.crs.to_epsg()

Data found here
Yukon_inlet_nodes.zip

It seems like one (but not the only) issue causing the broken build is changes to networkx. Might try adding a requirement that it be < v2.7.

It would be very useful to have channels which are not really "channels" removed from the channel mask for a more accurate DCN to be extracted. We want narrow islands to be removed but only channels who are too narrow for the channel they are in. This is because we want bifurcations, not flow around some sandbars that are large enough.

One way to do this could be removing any island whose length is less than the width of the channel it is. Could do this simply thru hole filling below a certain threshold. However we need to KNOW how wide the channel to do this.

My initial idea/way to tackle it.
Run thru channel extraction and get channel widths. Then compute the size of islands in every channel. Apply a threshold (L_I > L_C), for any island objects that don't meet this criteria remove them by filling in the holes in the channel mask. L_I could be max radius or could be sqrt(Area). Then have to basically use the cleaned up channel mask and start back from beginning?!?!

Jon you suggested something of working in raster space - maybe a trick using the distance transform to estimate island width and channel width, and then blanking out islands whose max dist is lower than than channel width up or downstream.

In some sample data I have a delta class with EPSG: 32603 (UTM Zone 3N), but the nodes and link that are written are in EPSG 4326 (unprojected). The skeleton is written to 32603. Not sure if this matters in pruning the nodes/links, but it makes visual inspection/overlaying of the skeleton, nodes, and links data in QGIS a multistep process. Found during debugging/probing why prune network produces an error while removing spurs, will follow up in another issue if necessary.

yukon.epsg
Out[40]: 32603

yukon.to_geovectors('network', ftype='json') # ftype can be either 'shp' or 'json'

C:\Anaconda3\envs\geo_py37\lib\site-packages\pyproj\crs\crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))

Data:
Yukon_data.zip

See below:

C:\Anaconda3\envs\RG_0.3\lib\site-packages\rivgraph\im_utils.py:456: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  out[prop] = np.array(allprop)

Jon, a nice (and easy) feature would be if a non-critical warning should pop if the CRS of the channel mask & the inletnodes and/or shoreline don't match. Now it is a simple reprojection without warning that it happened. Accidentally trimmed a Siberian delta with the Mackenzie and was befuddled why it was so bad.

Example in delta_utils.py:

 inlets_gpd = gpd.read_file(inlets_shp)
    # inlets_epsg = int(inlets_gpd.crs['init'].split(':')[1])
    inlets_epsg = inlets_gpd.crs.to_epsg()
    mask_epsg = gu.get_EPSG(gdobj)
    if inlets_epsg - mask_epsg != 0:
        inlets_gpd = inlets_gpd.to_crs(epsg=mask_epsg)

I see this in find_inlet_nodes and clip_by_shoreline, don't see anywhere else.

Geopandas released a new version: https://geopandas.org/docs/changelog.html#version-0-9-0-february-28-2021 RivGraph is currently pinned to 0.7 due to a pyproj dependency. This dependency should not be an issue in later versions of geopandas, and nothing in the release notes suggests potential issues. However, we should ensure that geopandas 0.9.0 causes no issues or warnings when testing RG.

The solidity property looks very useful for my application:

solidity (float)
Ratio of pixels in the region to pixels of the convex hull image.

Considering opening a PR to add it to im_utils.regionprops.

https://scikit-image.org/docs/stable/api/skimage.measure.html?highlight=regionprops#regionprops

I would like to establish a "showroom" of RivGraph applications. Basically a collection of peer-reviewed publications where RG was a major contributor to the work. Off the top of my head, we have

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019WR026867
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019GL086710
https://arxiv.org/abs/2108.04151 (we're working on resubmitting this)
Jay's flux partitioning work
Alex's eBI paper (preprint to be uploaded within weeks)
Tess mentioned something at AGU, not sure if it's published yet.

It would be cool to have a grid with one cell for each application. Title of the paper at the top (with a link to it), a figure/photo from the work, and a one-liner about how RG was used.

@elbeejay Do you think you could put together a template for this? It could live either in the Docs or not. I can populate it once it's in place.

I ran into the following issue. There is a link present which is so small that when I rasterize my links json in R that it literally has no points (some kind of interpolation issue probably). It's visible in the image below. So the link is caused because there are an excess number of links here in part due to these small islands. It has a length of 30sqrt(2), so basically 2 nodes which are diagonally adjacent on a 30-m resolution grid. I found 4 total links like this on the scene.

I see two solutions:

• better CN mask (remove small islands!) which will prevent this from occurring in the first place
• Some kind of link length checker in the prune_nodes? If you can compute length of line features and the length is < 1 or 2 pix or something should be merged with nearest link? Generally the algorithm would go:
  Compute all link lengths
  For all link lengths < cutoff (say 2 pix) {
  remove the link. Combine the two nodes to be one node
  Connections of new node will be downstream/upstream neighbors (should max be 4 neighbors?)

More examples below:

The installation is stuck during solving environment step!

Jon you had this fixed previously, not sure exactly why my shoreline has a column named LEFT_FID and RG is only checking for 'id_left'. This fix you sent to me works (for delta_utils.py, clip_by_shoreline).

    # Remove the links beyond the shoreline
    shore_int = gpd.sjoin(links_gdf, shore_gdf, op='intersects')
    # Find the proper column ID (left)
    leftid = [lid for lid in shore_int.columns if 'id' in lid.lower() and 'left' in lid.lower()][0]
    # Get ids of intersecting links
    cut_link_ids = shore_int[leftid].values

The (old) script I'm using to extract DCN's is only semi-working to get out a super apex. Leaving a note here to improve documentation and also get a working example to have a super apex.

from rivgraph.classes import delta
from matplotlib import pyplot as plt
import numpy as np
import rivgraph.im_utils as im
import rivgraph.mask_utils as mu
from rivgraph.io_utils import write_geotiff
import os
from osgeo import gdal
from datetime import datetime
import sys
import geopandas as gpd
import rivgraph.ln_utils as lnu
import pandas as pd

name = 'Pearl'
tht = 70

fn = name + "/RG_0.3/thresh_" + str(tht)
path_results = os.path.join(delta_dir, fn)
fn = name + "/" + name + "_filled.tif"
path_mask = os.path.join(delta_dir, fn)

delt = delta(name, path_mask, path_results, verbose=True)
delt.compute_network()
delt.to_geovectors('network')

path_shrln = os.path.join(delta_dir, name + "/outlines/OAM_shoreline.shp")
path_inlet = os.path.join(delta_dir, name + "/outlines/" + name + "_inlet_nodes.shp")

delt.prune_network(path_shoreline=path_shrln, path_inletnodes=path_inlet,
                   prune_less=True)
delt.compute_link_width_and_length()
delt.assign_flow_directions()
delt.links, delt.nodes = lnu.add_artificial_nodes(delt.links, delt.nodes, delt.gdobj)
delt.compute_link_width_and_length() # Need to re-compute because we've split some links into 2 in order to add artificial nodes.
delt.to_geovectors('network')
# Plot the links with the directionality marked
delt.plot('directions')

delt.to_geotiff('directions')
adj_w = delt.adjacency_matrix(weight='wid_adj') # Can also weight by 'len_adj' or provide a vector of your own weights.
print(adj_w)

path_adj = os.path.join(delta_dir, name + "/w_adj.csv")

np.savetxt(path_adj, adj_w, delimiter=",")

delt.save_network()

The output json files visualized in QGIS on the Pearl delta, with green circles representing nodes, pink lines representing nodes, and 3 brown circles the inlet nodes. A random node is highlighted (red).

Up to this point everything is fine except that delt doesn't contain a super apex node and has the 3 inlet nodes.

len(delt.nodes['id'])
Out[16]: 238

delt.nodes.keys()
Out[18]: dict_keys(['idx', 'conn', 'id', 'inlets', 'outlets', 'arts'])

delt.nodes['inlets']
Out[19]: [955, 928, 1043]

Running compute_topologic_metrics modifies delt in place and kills all but the widest inlet nodes:

delt.compute_topologic_metrics()

delt.nodes.keys()
Out[22]: dict_keys(['idx', 'conn', 'id', 'inlets', 'outlets', 'arts'])

delt.nodes['inlets']
Out[23]: []

len(delt.nodes['id'])
Out[24]: 212

Question

• What is RG's utility for returning the adjacency matrix modified to have a dummy row for the super apex? Where do I add add_super_apex at any point?
• compute_topologic_metrics modifies in place, but that only matters for the super apex case?

I'm running the delta_example notebook. Everything runs well until Section 4. Pruning the network:

colville.prune_network()
# Note that we can also specify the location of the shoreline and inlet nodes:
# colville.prune_network(path_shoreline='/path/to/shoreline/file', path_inletnodes='/path/to/inletnodes/file')

# Now that we've pruned, we should re-export the network:
colville.to_geovectors()
# Note that this time we didn't specify the arguments; by default 'network' will be exported as type 'json'.

---------------------------------------------------------------------------
CPLE_OpenFailedError                      Traceback (most recent call last)
fiona/_shim.pyx in fiona._shim.gdal_open_vector()

fiona/_err.pyx in fiona._err.exc_wrap_pointer()

CPLE_OpenFailedError: data/Colville_Delta/Results/Colville_inlet_nodes.shp: No such file or directory

During handling of the above exception, another exception occurred:

DriverError                               Traceback (most recent call last)
<ipython-input-8-2181e62ecfce> in <module>
----> 1 colville.prune_network()
      2 # Note that we can also specify the location of the shoreline and inlet nodes:
      3 # colville.prune_network(path_shoreline='/path/to/shoreline/file', path_inletnodes='/path/to/inletnodes/file')
      4 
      5 # Now that we've pruned, we should re-export the network:

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/rivgraph/classes.py in prune_network(self, path_shoreline, path_inletnodes)
    628             raise AttributeError('Could not inlet_nodes shapefile which should be at {}.'.format(self.paths['inlet_nodes']))
    629 
--> 630         self.links, self.nodes = du.prune_delta(self.links, self.nodes, path_shoreline, path_inletnodes, self.gdobj)
    631 
    632 

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/rivgraph/deltas/delta_utils.py in prune_delta(links, nodes, shoreline_shp, inlets_shp, gdobj)
     45     """
     46     # Get inlet nodes
---> 47     nodes = find_inlet_nodes(nodes, inlets_shp, gdobj)
     48 
     49     # Remove spurs from network (this includes valid inlets and outlets)

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/rivgraph/deltas/delta_utils.py in find_inlet_nodes(nodes, inlets_shp, gdobj)
    100 
    101     # Check that CRSs match; reproject inlet points if not
--> 102     inlets_gpd = gpd.read_file(inlets_shp)
    103     mask_crs = CRS(gdobj.GetProjection())
    104     if inlets_gpd.crs != mask_crs:

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/geopandas/io/file.py in read_file(filename, bbox, mask, rows, **kwargs)
     87 
     88     with fiona_env():
---> 89         with reader(path_or_bytes, **kwargs) as features:
     90 
     91             # In a future Fiona release the crs attribute of features will

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/fiona/env.py in wrapper(*args, **kwargs)
    398     def wrapper(*args, **kwargs):
    399         if local._env:
--> 400             return f(*args, **kwargs)
    401         else:
    402             if isinstance(args[0], str):

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/fiona/__init__.py in open(fp, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt, **kwargs)
    254 
    255         if mode in ('a', 'r'):
--> 256             c = Collection(path, mode, driver=driver, encoding=encoding,
    257                            layer=layer, enabled_drivers=enabled_drivers, **kwargs)
    258         elif mode == 'w':

~/miniconda3/envs/rivgraph/lib/python3.8/site-packages/fiona/collection.py in __init__(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, ignore_fields, ignore_geometry, **kwargs)
    160             if self.mode == 'r':
    161                 self.session = Session()
--> 162                 self.session.start(self, **kwargs)
    163             elif self.mode in ('a', 'w'):
    164                 self.session = WritingSession()

fiona/ogrext.pyx in fiona.ogrext.Session.start()

fiona/_shim.pyx in fiona._shim.gdal_open_vector()

DriverError: data/Colville_Delta/Results/Colville_inlet_nodes.shp: No such file or directory

I believe the correct file is DriverError: data/Colville_Delta/Colville_inlet_nodes.shp. I am running the notebook from within the examples/ directory, and I haven't modified any of the code cells.

Did I do something wrong? Should I execute the notebook from a different directory? (If so, perhaps it would be useful to document which directory)

Seems like there can be instances when remove_two_link_nodes() is called (during network pruning), and during the node connectivity updating sometimes the function tries to "pop" an empty set (Line 932) which throws an error.

I think a quick fix could be sticking those few lines of code in a try-except statement, although I will look into the root cause of the bug to properly resolve it.

I want to install RivGraph with conda and environment.yml,how anaconda prompt tells me this:

(base) C:\Users\XXX> conda env create C:\Users\XXX\rivgraph\enviroment.yml
Collecting package metadata (repodata.json): done
Solving environment: failed

ResolvePackageNotFound:
  - esmpy!=8.1.0
  - subversion

Why it can't find esmpy and subversion?

im_utils.regionprops() prints this warning upon each initial call:

...\rivgraph\im_utils.py:657: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  out[prop] = np.array(coords)

We should fix this.

I'm encountering a Python-based issue with how to update the package from source. I've followed the instructions to install RG using from the new installation instructions . The installation appears to be successful, but pytest fails. Installation and failure below. Moreover, when I run RG or inspect the source of recently modified functions, I don't see the latest version of the function (see second code block below installation).

Installation:

HOME>conda activate RG_0.3

(RG_0.3) HOME>cd documents\deltas\RG_2021_2_23\RivGraph

(RG_0.3) HOME\Documents\deltas\RG_2021_2_23\RivGraph>python setup.py install
running install
running bdist_egg
running egg_info
writing rivgraph.egg-info\PKG-INFO
writing dependency_links to rivgraph.egg-info\dependency_links.txt
writing top-level names to rivgraph.egg-info\top_level.txt
reading manifest file 'rivgraph.egg-info\SOURCES.txt'
writing manifest file 'rivgraph.egg-info\SOURCES.txt'
installing library code to build\bdist.win-amd64\egg
running install_lib
running build_py
creating build\bdist.win-amd64\egg
creating build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\classes.py -> build\bdist.win-amd64\egg\rivgraph
creating build\bdist.win-amd64\egg\rivgraph\deltas
copying build\lib\rivgraph\deltas\delta_directionality.py -> build\bdist.win-amd64\egg\rivgraph\deltas
copying build\lib\rivgraph\deltas\delta_metrics.py -> build\bdist.win-amd64\egg\rivgraph\deltas
copying build\lib\rivgraph\deltas\delta_utils.py -> build\bdist.win-amd64\egg\rivgraph\deltas
copying build\lib\rivgraph\deltas\__init__.py -> build\bdist.win-amd64\egg\rivgraph\deltas
copying build\lib\rivgraph\directionality.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\geo_utils.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\im_utils.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\io_utils.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\ln_utils.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\mask_to_graph.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\mask_utils.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\ordered_set.py -> build\bdist.win-amd64\egg\rivgraph
creating build\bdist.win-amd64\egg\rivgraph\rivers
copying build\lib\rivgraph\rivers\centerline_utils.py -> build\bdist.win-amd64\egg\rivgraph\rivers
copying build\lib\rivgraph\rivers\river_directionality.py -> build\bdist.win-amd64\egg\rivgraph\rivers
copying build\lib\rivgraph\rivers\river_utils.py -> build\bdist.win-amd64\egg\rivgraph\rivers
copying build\lib\rivgraph\rivers\__init__.py -> build\bdist.win-amd64\egg\rivgraph\rivers
copying build\lib\rivgraph\walk.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\_plotting_metrics.py -> build\bdist.win-amd64\egg\rivgraph
copying build\lib\rivgraph\__init__.py -> build\bdist.win-amd64\egg\rivgraph
byte-compiling build\bdist.win-amd64\egg\rivgraph\classes.py to classes.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\deltas\delta_directionality.py to delta_directionality.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\deltas\delta_metrics.py to delta_metrics.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\deltas\delta_utils.py to delta_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\deltas\__init__.py to __init__.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\directionality.py to directionality.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\geo_utils.py to geo_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\im_utils.py to im_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\io_utils.py to io_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\ln_utils.py to ln_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\mask_to_graph.py to mask_to_graph.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\mask_utils.py to mask_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\ordered_set.py to ordered_set.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\rivers\centerline_utils.py to centerline_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\rivers\river_directionality.py to river_directionality.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\rivers\river_utils.py to river_utils.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\rivers\__init__.py to __init__.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\walk.py to walk.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\_plotting_metrics.py to _plotting_metrics.cpython-38.pyc
byte-compiling build\bdist.win-amd64\egg\rivgraph\__init__.py to __init__.cpython-38.pyc
creating build\bdist.win-amd64\egg\EGG-INFO
copying rivgraph.egg-info\PKG-INFO -> build\bdist.win-amd64\egg\EGG-INFO
copying rivgraph.egg-info\SOURCES.txt -> build\bdist.win-amd64\egg\EGG-INFO
copying rivgraph.egg-info\dependency_links.txt -> build\bdist.win-amd64\egg\EGG-INFO
copying rivgraph.egg-info\not-zip-safe -> build\bdist.win-amd64\egg\EGG-INFO
copying rivgraph.egg-info\top_level.txt -> build\bdist.win-amd64\egg\EGG-INFO
creating 'dist\rivgraph-0.3-py3.8.egg' and adding 'build\bdist.win-amd64\egg' to it
removing 'build\bdist.win-amd64\egg' (and everything under it)
Processing rivgraph-0.3-py3.8.egg
removing 'c:\anaconda3\envs\rg_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg' (and everything under it)
creating c:\anaconda3\envs\rg_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg
Extracting rivgraph-0.3-py3.8.egg to c:\anaconda3\envs\rg_0.3\lib\site-packages
rivgraph 0.3 is already the active version in easy-install.pth

Installed c:\anaconda3\envs\rg_0.3\lib\site-packages\rivgraph-0.3-py3.8.egg
Processing dependencies for rivgraph==0.3
Finished processing dependencies for rivgraph==0.3

(RG_0.3) HOME\Documents\deltas\RG_2021_2_23\RivGraph>pytest
ImportError while loading conftest 'HOME\Documents\deltas\RG_2021_2_23\RivGraph\tests\conftest.py'.
tests\conftest.py:21: in <module>
    from rivgraph.classes import delta
rivgraph\__init__.py:3: in <module>
    from rivgraph.classes import delta, river
rivgraph\classes.py:9: in <module>
    import gdal
C:\Anaconda3\lib\site-packages\gdal.py:2: in <module>
    from osgeo.gdal import deprecation_warn
C:\Anaconda3\lib\site-packages\osgeo\gdal.py:1707: in <module>
    from . import ogr
C:\Anaconda3\lib\site-packages\osgeo\ogr.py:241: in <module>
    from . import osr
C:\Anaconda3\lib\site-packages\osgeo\osr.py:26: in <module>
    _osr = swig_import_helper()
C:\Anaconda3\lib\site-packages\osgeo\osr.py:22: in swig_import_helper
    _mod = imp.load_module('_osr', fp, pathname, description)
C:\Anaconda3\lib\imp.py:242: in load_module
    return load_dynamic(name, filename, file)
C:\Anaconda3\lib\imp.py:342: in load_dynamic
    return _load(spec)
E   ImportError: DLL load failed: The specified procedure could not be found.

(RG_0.3) HOME\Documents\deltas\RG_2021_2_23\RivGraph>

Functions aren't updated:

import rivgraph.deltas.delta_metrics as dm

import inspect

ps = inspect.getsource(dm.graphiphy)

print(ps)
def graphiphy(links, nodes, weight=None):
    """Create a networkx graph."""
    if weight is not None and weight not in links.keys():
        raise RuntimeError('Provided weight key not in nodes dictionary.')

    if weight is None:
        weights = np.ones((len(links['conn']), 1))
    else:
        weights = links[weight]

    G = nx.DiGraph()
    G.add_nodes_from(nodes['id'])
    for lc, wt in zip(links['conn'], weights):
        G.add_edge(lc[0], lc[1], weight=wt)

    return G

io_utils.py produces an error due to fiona not explicitly recognizing crs. Fixed by including

import fiona.crs

Need to investigate and fix hosted documentation to show the API Reference information. Example empty page here.

I am trying to prune the network on the Yukon data that I have been using and continuously get this error in the remove_two_link_nodes. To get here, I extract the network up to prune_network as in the example (same as all other comments this week). I then run the functions in remove_all_spurs in the global environment (i.e. not the function remove all spurs but the manual parts). Check 2nd code block.

 ct = 0
        # Remove spurs
        for nid, con in zip(nodes['id'], nodes['conn']):
            
            if len(con) == 1 and nid not in dontremove:
                print(nid)
                ct = ct + 1
                links, nodes = delete_link(links, nodes, con[0])
                
        # Remove self-looping links (a link that starts and ends at the same node)
        for nid, con in zip(nodes['id'], nodes['conn']):
            m = mode(con)
            if m.count[0] > 1:
                
                # Get link 
                looplink = m.mode[0]                
                
                # Delete link
                links, nodes = delete_link(links, nodes, looplink)
                ct = ct + 1             
        
            linkkeys = [lk for lk in links.keys() if type(links[lk]) is not int and len(links[lk]) == len(links['id'])]
ct = 0
        for nidx, nid in enumerate(nodes['id']):
            # Get connectivity of current node        
            conn = nodes['conn'][nidx][:]
            # We want to combine links where a node has only two connections
            if len(conn) == 2 and nid not in dontremove:
                print(nid)
                # Delete the node
                nodes = delete_node(nodes, nid, warn=False)
                
                # The first link in conn will be absorbed by the second
                lid_go = conn[0]
                lid_stay = conn[1]
                lgo_idx = links['id'].index(lid_go)
                lstay_idx = links['id'].index(lid_stay)
                
                # Update the connectivity of the node attached to the link being absorbed
                conn_go = links['conn'][lgo_idx]
                conn_stay = links['conn'][lstay_idx]
                
                # Update node connectivty of go link (stay link doesn't need updating)
                node_id_go = (set(conn_go) - set([nid])).pop()
                nodes['conn'][nodes['id'].index(node_id_go)].remove(lid_go)
                nodes['conn'][nodes['id'].index(node_id_go)].append(lid_stay)
                
                # Update link connectivity of stay link
                conn_go.remove(nid)
                conn_stay.remove(nid)
                # Add the "go" link connectivity to the "stay" link
                conn_stay.extend(conn_go)
                
                # Update the indices of the link
                idcs_go = links['idx'][lgo_idx]
                idcs_stay = links['idx'][lstay_idx]
                if idcs_go[0] == idcs_stay[-1]:
                    new_idcs = idcs_stay[:-1] + idcs_go
                elif idcs_go[-1] == idcs_stay[0]:
                    new_idcs = idcs_go[:-1] + idcs_stay
                elif idcs_go[0] ==  idcs_stay[0]:
                    new_idcs = idcs_stay[::-1][:-1] + idcs_go
                elif idcs_go[-1] == idcs_stay[-1]:
                    new_idcs = idcs_stay[:-1] + idcs_go[::-1]
                links['idx'][lstay_idx] = new_idcs
                
                # Delete the "go" link
                lidx_go = links['id'].index(lid_go)
                for lk in linkkeys:
                    if lk == 'id': # have to treat orderedset differently
                        links[lk].remove(lid_go)
                    elif type(links[lk]) is np.ndarray: # have to treat numpy arrays differently
                        links[lk] = np.delete(links[lk], lid_go)
                    else:
                        links[lk].pop(lidx_go)
                                
                ct = ct + 1

Traceback (most recent call last):

  File "<ipython-input-297-4805c3adaf0e>", line 22, in <module>
    node_id_go = (set(conn_go) - set([nid])).pop()

KeyError: 'pop from an empty set'

This happens to node 29223:

nid
Out[309]: 29223

conn
Out[310]: [16301, 16301]

lid_go
Out[311]: 16301

lid_stay
Out[313]: 16301

It seems the prior call to remove self loops should have removed it. This seems like a RG bug, although not sure. Originally this corresponded to a node with 3 links, one of which is a spur... my python skills are not up there to get this taken care of quickly, and you may be able to get a solution to this faster.

I previously uploaded some links/nodes, not sure if those had been partially trimmed. I am going to upload again:
links, nodes, skeleton (all before pruning):
yukon_beforepruning.zip
the mask used:
compoun_binary.zip

I have attempted to reproduce the example provided for the Colville using a binary mask of the Yukon (can attach if necessary). When I reach yukon.plot('network'), having run the example code modified for the Yukon, I get an error about axis not being found. Code and error reproduced below. This occurred both in a modified jupyter notebook and in a spyder environment.

from rivgraph.rivgraph import delta
import matplotlib.pyplot as plt
import os

# Define the path to the georeferenced binary image.
mask_path = r"C:\Users\Lawrence Vulis\Documents\deltas\r_work\UTM\Yukon\compoun_binary.tif"

# Results will be saved with this name
name = 'Yukon' 

# A folder called Yukon will be created within this path for storing outputs
results_folder = r"C:\Users\Lawrence Vulis\Documents\deltas\RivGraph-master\LV_extraction\Results" 

# Boot up the delta class! We set verbose=True to see progress of processing.
yukon = delta(name, mask_path, results_folder=results_folder, verbose=True) 

# The mask has been re-binarized and stored as an attribute of yukon:
plt.imshow(yukon.Imask)


# Simply use the skeletonize() method.
yukon.skeletonize()

# After running, yukon has a new attribute: Iskel. Let's take a look.
plt.imshow(yukon.Iskel)


# The skeleton is hard to see; perhaps we'd like to look at it closer?
# One option is to save it as a geotiff and pull it up in a GIS (like QGIS).
# We use the write_geotiff() method with the "skeleton" option.
yukon.to_geotiff('skeleton')



# Simply use the compute_network() method.
yukon.compute_network()


# Now we can see that the "links" and "nodes" dictionaries have been added as yukon attributes:
links = yukon.links
nodes = yukon.nodes
print('links: {}'.format(links.keys()))
print('nodes: {}'.format(nodes.keys()))


yukon.plot('network')

yukon.plot('network')
Traceback (most recent call last):

  File "<ipython-input-7-23d8c36cc499>", line 1, in <module>
    yukon.plot('network')

  File "C:\Anaconda3\envs\geo_py37\lib\site-packages\rivgraph-0.2-py3.7.egg\rivgraph\rivgraph.py", line 263, in plot
    f = lnu.plot_network(self.links, self.nodes, self.Imask, self.name, axis=axis)

TypeError: plot_network() got an unexpected keyword argument 'axis'

Hi all, I'm trying to create a binary channel mask for a river in Alaska to run through RivGraph. I have created two channel masks - initially, I created one using the NAD 1927 Albers coordinate system, but then decided that using UTM Zone 6N was more appropriate, so created a mask with this CRS too. The binary channel masks covers an area of approximately 25 square km and are at 10m resolution, in 2-Bit TIFF format.

The issue I'm having is slightly different for both of the masks. For the original NAD Albers one, the code runs successfully up until the mesh creation section, where I get a shapely error (despite being on version 1.7.1). The error says "NotImplementedError: Multi-part geometries do not provide a coordinate sequence". I also get this error with the UTM mask, although there seems to be other issues with this one too. There is an issue with the scales where 'Instantiate River Mask' shows the river reaching a scale of 500 x 400, which then shrinks to 500 x 300 at the 'Skeletonization' step. The network and nodes are also stretched across the binary mask and in totally the wrong place. I suspected this was due to NoData cells in the binary mask, but I can't seem to find any so I'm not sure what is causing all these scale problems. See images below.

I'm currently on a Windows computer, but also have a Mac I can test things on. I'm using the following package versions:

• conda=4.10.1
• python=3.9.15
• geopandas=0.9.0 (I was getting a GeoDataFrame error with the latest version)
• scikit-image=0.19.3
• opencv=4.5.5
• networkx=2.6.3
• shapely=1.7.1
• matplotlib=3.6.2
• fastdtw=0.3.4
• loguru=0.6.0
• scipy=1.8.1
• numpy=1.22.4

I have tried creating multiple masks in the UTM Zone 6N coordinate system, some with 2m resolution, and still getting the same issues.

Any advice or help on this would be hugely appreciated.

Thank you!

I encounter the following error from pruning the network with the shoreline & inlet nodes I have. I am documenting this now in case you may know, but will look at it more further (bad shoreline, bad inlet nodes, etc.). If one of those, maybe this error type can give rise to good documentation. Apologies for all the late night notifications.

...: ### Pruning the network

yukon.prune_network(path_shoreline='C:/Users/Lawrence Vulis/Documents/deltas/r_work/utm/yukon/outlines/Yukon_shoreline_v1.shp', path_inletnodes='C:/Users/Lawrence Vulis/Documents/deltas/r_work/utm/yukon/outlines/Yukon_inlet_nodes.shp')
Geotiff written to C:\Users\Lawrence Vulis\Documents\deltas\RG_experiments\LV_extraction\Results\Yukon\Yukon_skel.tif.
Resolving links and nodes...done.
links: dict_keys(['idx', 'conn', 'id', 'n_networks'])
nodes: dict_keys(['idx', 'conn', 'id'])
C:\Anaconda3\envs\geo_py37\lib\site-packages\pyproj\crs\crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))
C:\Users\Lawrence Vulis\Documents\deltas\RG_experiments\LV_extraction\Results\Yukon\Yukon_links.json
C:\Users\Lawrence Vulis\Documents\deltas\RG_experiments\LV_extraction\Results\Yukon\Yukon_nodes.json
Traceback (most recent call last):

  File "<ipython-input-356-11a37adc6754>", line 59, in <module>
    yukon.prune_network(path_shoreline='C:/Users/Lawrence Vulis/Documents/deltas/r_work/utm/yukon/outlines/Yukon_shoreline_v1.shp', path_inletnodes='C:/Users/Lawrence Vulis/Documents/deltas/r_work/utm/yukon/outlines/Yukon_inlet_nodes.shp')

  File "C:\Anaconda3\envs\geo_py37\lib\site-packages\rivgraph-0.2-py3.7.egg\rivgraph\rivgraph.py", line 513, in prune_network
    self.links, self.nodes = du.prune_delta(self.links, self.nodes, path_shoreline, path_inletnodes, self.gdobj)

  File "C:\Anaconda3\envs\geo_py37\lib\site-packages\rivgraph-0.2-py3.7.egg\rivgraph\deltas\delta_utils.py", line 56, in prune_delta
    links, nodes = clip_by_shoreline(links, nodes, shoreline_shp, gdobj)

  File "C:\Anaconda3\envs\geo_py37\lib\site-packages\rivgraph-0.2-py3.7.egg\rivgraph\deltas\delta_utils.py", line 174, in clip_by_shoreline
    lidx = links['id'].index(clid)

  File "C:\Anaconda3\envs\geo_py37\lib\site-packages\rivgraph-0.2-py3.7.egg\rivgraph\ordered_set.py", line 194, in index
    return self.map[key]

KeyError: -1

At the end of the About section in the README, there is link to a "delta channel network" that links to https://github.com/jonschwenk/RivGraph/blob/master/examples/delta_example.py.ipynb. This leads to a 404 on my end.