Plotting smoothed stream object results

[fluviotect]

I'm trying to wrap my brain around the dimensionality of the CRS output and what it represents in the real-world.

The entire dem (a 2D object) is the CRS argument run on the stream object (also 2D) to produce a 1D array.

So is the dimension reduction of the CRS just stacking the non-null values either row-wise or column-wise? Something else?

How best to unpack it, and plot stream profiles?

Cheers,
Will

[wkearn]

The output of the smoothing functions is what we call a "node attribute list," which is a 1D array with one element for every node in the stream network. These nodes correspond to pixels of the source digital elevation model, but not every pixel is in the stream network. The StreamObject.shape attribute is the shape of this DEM. The easiest way to figure out how many nodes are in the stream network is to call s.stream.shape on your StreamObject, but any function that returns a node attribute list will return an array of this length. This includes crs as well as downstream_distance and streamorder among others.

To plot profiles, there is a method StreamObject.plotdz that takes either a GridObject like a DEM or a node attribute list like the output of crs and plots it against upstream distance:

dem = tt3.load_dem('bigtujunga')
fd = tt3.FlowObject(dem)
s = tt3.StreamObject(dem)
s_smth = s.crs(dem)
ax = s.plotdz(dem, color='grey') # dem is a GridObject
ax = s.plotdz(s_smth, ax=ax, color='black') # s_smth is a node attribute list
ax.autoscale_view()

There are also other ways to plot the node attribute lists such as plotting them on the stream network in map view with colored values. This is currently a lot more complicated because of limitations of matplotlib, but I can provide an example if you'd like.

[fluviotect]

Thank-you, that's great.

The relationship between the (various) node attribute lists and the parent nodes (s.streams) is presumably the ordering (index values) of the lists/arrays, yes?

Eventually I will want to plot different colours, but my more immediate question is how to identify and query different node subsets?

For example:

• by exporting the entire network to a GIS format, adding a binary field to the attribute table that indicates interest/non-interest, and then importing updated shp back to TT3 for use as a mask
• by specifying a downstream pourpoint

Once I understand that, I should be able to sort-out differential colouring, etc.

I'm happy to start a different Discussion thread if you would like to keep items more discrete.

[wkearn]

The relationship between the (various) node attribute lists and the parent nodes (s.streams) is presumably the ordering (index values) of the lists/arrays, yes?

s.stream is an array of the linear indices of pixels in the DEM of each node, but this array is best thought of as an implementation detail -- the values in the array and their ordering might change in the future.

Eventually I will want to plot different colours, but my more immediate question is how to identify and query different node subsets?

This depends on the format of the data that you want to use to query the network. I've given some examples below, but let me know if you need something else.

For example:

* by exporting the entire network to a GIS format, adding a binary field to the attribute table that indicates interest/non-interest, and then importing updated shp back to TT3 for use as a mask

This is possible, but might not be as intuitive as it could be. The StreamObject.to_geodataframe and StreamObject.to_shapefile methods will export the StreamObject to a vector dataset that you could use in a GIS, but right now they export the StreamObject as a collection of line segments very similar to what we do when we plot a StreamObject in Python. It is more useful for visualization than analysis. We're currently working on providing better export functionality for StreamObject.

You can also extract the x and y coordinates of the nodes:

import numpy as np
j, i = s.node_indices
xs, ys = self.transform * np.vstack((i, j))

xs and ys are node attribute lists with the x- and y- coordinates of nodes in the stream network, which you can export as you need.

* by specifying a downstream pourpoint

If your pour point is a set of geographic coordinates, you'll need to find the pixel in the DEM or node in the stream network that is closest to those coordinates. MATLAB TopoToolbox has a few functions like snap2stream and coord2ind that can help with this, but we haven't ported these to pytopotoolbox yet. If that would be helpful, we can prioritize it.

It is easier if you have a pour point or some other data labeled in a raster dataset that has the same size and georeferencing as the source DEM for your StreamObject. We do have a few functions in pytopotoolbox for interfacing between the gridded datasets and the node attribute lists.

StreamObject.ezgetnal takes a gridded data set and gives you a node attribute list with the values of the grid at each of the nodes. For example:

dem = tt3.load_dem('bigtujunga')

fd = tt3.FlowObject(dem)
db = fd.drainagebasins()

s = tt3.StreamObject(fd)
db_s = s.ezgetnal(db)

db is a GridObject matching the size of dem containing drainage basin labels, and db_s is a node attribute list of drainage basin labels. Now you could filter the stream network by a specific drainage basin label using StreamObject.subgraph, which takes a Boolean node attribute list and returns a stream network containing only those nodes where the node attribute list is True. For example, you could extract the network from the basin with the largest area by counting unique values in the drainage basins raster:

values, counts = np.unique(db, return_counts=True)
s2 = s.subgraph(db_s == values[np.argmax(counts)])

You can do this for any property that you can compute on a node attribute list. If you do go the route of exporting the coordinates and labeling them elsewhere, you can also import that back into Python and use StreamObject.subgraph to extract the desired portion of the stream network.

There is also an StreamObject.upstreamto and a StreamObject.downstreamto method. These take Boolean node attribute lists or rasters and compute return the portion of the stream network upstream or downstream of any True values. If you create a 2D array where your pour points are True and everything else is False, s.upstreamto(pourpoints) would give you a new StreamObject that only has points upstream of your pour points.

[fluviotect]

What a great explanation. Thank-you.

Definitely not as straight forward as I'd hoped. But lots of options and all achievable.

I see a new issue to expose StreamObject coordinates more explicitly. That will be good.

For now, am having a play with the geodataframe and node export approach as I'm already familiar with that approach. Appreciate your time in the various questions/issues I've raised the past few days. Cheers