<jupyter_start><jupyter_text>Channel profiles and working with the `_MChiSegmented` file *This lesson has been written by <NAME> at the University of Edinburgh**Last update 01/03/2021* In a few of the notebooks we have used `lsdtopotools` to extract a channel profile. The software has a few different ways to extract profiles, but one produces a file with `_MChiSegmented` in the filename. This file has a selection of useful data in it. Lets load the file and have a look.  Load packages and then load the data<jupyter_code>import pandas as pd
import geopandas as gpd
import cartopy as cp
import cartopy.crs as ccrs
import rasterio as rio
import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>We are going to load the csv data, which we generated in the previous notebook. We first use a package called `pandas` to do this, and then import that into `geopandas`. `pandas` is for dealing with data that you might use excel to look at. So data that might be in a spreadsheet. But it has many powerful data processing options. It was originally developed by a "quant" financial analyst to look at stock trends! But it is a super useful python package for many applications.The dataset you load is in the filename in quotes in line 1 of the cell below (if you can't see line numbers, go to `view` -> `toggle line numbers` in the menu at the top of the page. So in this example the filename is "SorbasChi_MC_MChiSegmented.csv"<jupyter_code>df = pd.read_csv("SorbasChi_MC_MChiSegmented.csv")
gdf = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df.longitude, df.latitude))

# We have to tell the geopandas data what geographic system we are in by using something called an EPSG code. 
# All major geographic projection and transformation system have this code. 
gdf.crs = "EPSG:4326" 

# The head command shows you what is in the file.
print(gdf.head())<jupyter_output><empty_output><jupyter_text>Okay, I am not going to explain "chi", or "M_chi" or "b_chi" yet. Most other things should be self explanatory. The exceptions are `source_key` and `basin_key`:* Each basin is tagged by a `basin_key` so if you pick this, then you isolate the basin. * Each channel has a "source" at the tip of the channel. These have a source key. So if you pick this you have the flow path from that source to the outlet.   Making new data columns: slope and smoothed slope Okay, we have flow distance and elevation in this file, but we also want to look at the slope of the channel. To get the slope, we need to calculate the change in elevation over the change in flow distance. The mathematical operation for this is called the gradient (or, if you want to use the notation of derivatives it is `dz/dx`).The python package `numpy` has a built in function for calculating the gradient (`np.gradient`), which we use below to get the slope along the channel.<jupyter_code>z = gdf.elevation
x = gdf.flow_distance
S = np.gradient(np.asarray(z),np.asarray(x))
gdf["slope"] = S
gdf.head()<jupyter_output><empty_output><jupyter_text>In this notebook, I want to isolate one of the basins. `pandas` has lots of easy ways to isolate data.In line 2 below, you use this syntax to isolate data in your data set. The `basin_key` is one of the data columns in the dataset (you can always see the data columns by using the `.head()` command on your dataset.<jupyter_code># First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]
gdf_b1.head()<jupyter_output><empty_output><jupyter_text>Okay, now lets plot the channel profile and the channel slope. Much of the code below is for plotting, but there are some key lines (lines 7-8):     The main stem channel is the one with the minimum source key in this basin    min_source = np.amin(gdf_b1.source_key)    gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]    In thse three lines, we are isolating a `source_key` which is a single channel in the DEM. The basin has lots of channels, but we only want one. The minimm source key in a basin is the longest channel in this dataset.<jupyter_code>plt.rcParams['figure.figsize'] = [10, 10]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
S = gdf_b2.slope

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1, ax2) = plt.subplots(2, 1)
ax1.scatter(x_locs, z,s = 0.2)
ax2.scatter(x_locs, S,s = 1)


ax1.set_xlabel("Distance from outlet ($m$)")
ax1.set_ylabel("elevation (m)")

ax2.set_xlabel("Distance from outlet ($m$)")
ax2.set_ylabel("Slope (m/m)")

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>This slope (bottom figure) is very noisy. One way to deal with this is to smooth the data. We can smooth the data by running a mobing window over it and doing some averaging inside the window. Python has lots of tools for this. In this case I use a `rolling` window and I have picked various settings. You don't need to worry about this too much, the only number that you might wanty to play with is the first number after `rolling` which is the number of datapoints in the window. The bigger this number, the more smoothed the data becomes.<jupyter_code>gdf['slope_rolling'] = gdf.slope.rolling(40,win_type='hamming').mean()
gdf.head()<jupyter_output><empty_output><jupyter_text>This plot will show the slope and the rolling slope, so you can see how the rolling window smooths the data.<jupyter_code>plt.rcParams['figure.figsize'] = [10, 10]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
# To get a single data column from a pandas dataframe (in this case called gdf_b2) you just put
# a full stop and then the name of the column
# If your column has spaces or funny characters in the name you need to use the square brackets like this:
# z = gdf_b2["elevation"]
# Which is an alternative way of isolating data
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
S = gdf_b2.slope
SR = gdf_b2.slope_rolling

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1, ax2) = plt.subplots(2, 1)
ax1.scatter(x_locs, S,s = 1)
ax2.scatter(x_locs, SR,s = 1)


ax1.set_xlabel("Distance from outlet ($m$)")
ax1.set_ylabel("Slope (m/m)")

ax2.set_xlabel("Distance from outlet ($m$)")
ax2.set_ylabel("Rolling Slope (m/m)")

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Looking at the gradient and where the high gradient channels are along the channel profile Okay, the rolling slope allows us to see some spikes in the gradient. Can we see this in the right places along the channel profile?<jupyter_code>plt.rcParams['figure.figsize'] = [10, 5]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
# If you want to play with this a bit you can change the source number to look at different channels
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
S = gdf_b2.slope
SR = gdf_b2.slope_rolling

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1) = plt.subplots(1, 1)
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

# Make the scatter plots
ax1.scatter(x_locs, z,s = 1, label='Longitudinal profile')
ax2.scatter(x_locs, SR,s = 1,c="r", label='Channel slope')

# Some code to make sure the legend renders on the same axis
lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2, loc=0)

ax1.set_xlabel("Distance from outlet ($m$)")
ax1.set_ylabel("Elevation (m)")

ax2.set_xlabel("Distance from outlet ($m$)")
ax2.set_ylabel("Rolling Slope (m/m)")

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Well, it looks like those really steep bits in the channel profile are showing up gradient data. That is good: at least we know the data looks sensible. But, in the headwaters, there are some steep sections. Can we tell if these are steeper than the rest of the channel if we control for drainage area?The traditional way to do this was to use a plot of the logarithm of slope against the logarithm of drainage area. The theory behind that is explained in the lessons in the `Channel_incision` directory. Here I am just going to show you the data.   Slope and area<jupyter_code>plt.rcParams['figure.figsize'] = [10, 10]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
S = gdf_b2.slope
SR = gdf_b2.slope_rolling
A = gdf_b2.drainage_area

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1, ax2) = plt.subplots(2, 1)
ax1.scatter(x_locs, SR,s = 1)
ax2.scatter(A, SR,s = 1)


ax1.set_xlabel("Distance from outlet ($m$)")
ax1.set_ylabel("Slope (m/m)")

ax2.set_xlabel("Distance from outlet ($m$)")
ax2.set_ylabel("Rolling Slope (m/m)")
ax2.set_xscale("log")
ax2.set_yscale("log")
ax2.set_ylim(0.001,0.2)

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>So the plot above is the slope-area plot. As you can see, it is very messy. If you want to get a sensible S-A relationship out of this you need to do some more work of binning and averaging the data. Remember we have already smoothed the data! So getting a clean slope-area plot requires massaging the data quite a lot. The problem here is that slope data varies a lot along the channel. Slope is noisy because channels are naturally bumpy, and also because there are errors in the topographic data. In addition, drainage area jumps every time there is a junction (you add the two drainage areas of the tributaries together) so your area data has big gaps. So in the traditional way to analyse channels, figuring out what part of the channel is relatively steep (that is what is steeper after you control for drainage area) is quite challenging. <NAME>, with collaborators, developed a method to get around this noise: she presented it at a conference in 2000 and then never published the work in a journal. She then published it with another author 13 years later:https://onlinelibrary.wiley.com/doi/abs/10.1002/esp.3302<NAME>., <NAME>., 2013. Solutions of the stream power equation and application to the evolution of river longitudinal profiles. Journal of Geophysical Research: Earth Surface 118, 497–518. https://doi.org/10.1002/jgrf.20031The theory is explained in other notebooks (`channel_incision/Lesson_06_integral_transform.ipynb`). We are going to concentrate on the outcome of that work.   Steepness data with chi profiles<jupyter_code>plt.rcParams['figure.figsize'] = [10, 10]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
chi = gdf_b2.chi

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1, ax2) = plt.subplots(2, 1)
ax1.scatter(x_locs, z,s = 0.2)
ax2.scatter(chi, z,s = 0.2)


ax1.set_xlabel("Distance from outlet (m)")
ax1.set_ylabel("elevation (m)")

ax2.set_xlabel("$\chi$ (m)")
ax2.set_ylabel("elevation (m)")

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>What is this? In short, the $\chi$ coordinate is _transformed_ (a fancy way of saying we applied some maths to it) and the transformation included drainage area. So the $\chi$ coordinate is like the flow distance but we have wrapped drainage area into it.This coordinate was **specifically designed** so that the gradient in the "chi plot" (that is, the bottom plot which is a plot a $\chi$ vs elevation) corresponds to the relative steepness of the channel. And by relative steepness I mean steepness normalised for drainage area.As it turns out, the program that produced the data also calculates the gradient in the $\chi$ vs elevation plot. This gradient is called `m_chi`. There is a bunch of statistics behind the way it does it, and if you are in need of some bedtime reading in order to battle insomnia, you can read the details here:*<NAME>., <NAME>., <NAME>., <NAME>., <NAME>., 2014. A statistical framework to quantify spatial variation in channel gradients using the integral method of channel profile analysis. Journal of Geophysical Research: Earth Surface 119, 138–152. https://doi.org/10.1002/2013JF002981*In the plot below we look at that data:<jupyter_code>plt.rcParams['figure.figsize'] = [10, 5]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
# If you want to play with this a bit you can change the source number to look at different channels
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
m_chi = gdf_b2.m_chi
chi = gdf_b2.chi

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1) = plt.subplots(1, 1)
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

# Make the scatter plots
ax1.scatter(x_locs, z,s = 1, label='Longitudinal profile')
ax2.scatter(x_locs, m_chi,s = 1,c="r", label='Channel steepness')

# Some code to make sure the legend renders on the same axis
lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2, loc=0)

# These are just the labels for the figure
ax1.set_xlabel("Distance from outlet ($m$)")
ax1.set_ylabel("Elevation (m)")

ax2.set_xlabel("Distance from outlet ($m$)")
ax2.set_ylabel("Channel steepness ($k_{sn}$)")

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Okay, this plot shows the channel profile alongside the channel steepness. This resembles the slope plot, but these steepnesses are normalised for drainage area. What this is showing is that1. The top of the channel is not steep relative to the bottom of the channel!2. The are some very steep bits in the middle. We call the knickpoints. 3. There seem to be three prominent knickpoints in the middle of the channel (at around 30km from the outlet and also at around 45 km from the outlet). 4. The channels below the knickpoints are a lot steeper than channels above the knickpoints. We can also skip having two plots and just colour the points in the channel profile by their steepness, which I do below:<jupyter_code>plt.rcParams['figure.figsize'] = [10, 5]

# First lets isolate just one of these basins. There is only basin 0 and 1
gdf_b1 = gdf[(gdf['basin_key'] == 0)]

# The main stem channel is the one with the minimum source key in this basin
# If you want to play with this a bit you can change the source number to look at different channels
min_source = np.amin(gdf_b1.source_key)
gdf_b2 = gdf_b1[(gdf_b1['source_key'] == min_source)]
#gdf_b2 = gdf_b1

# Now make channel profile plots
z = gdf_b2.elevation
x_locs = gdf_b2.flow_distance
m_chi = gdf_b2.m_chi
chi = gdf_b2.chi

# Create two subplots and unpack the output array immediately
plt.clf()
f, (ax1) = plt.subplots(1, 1)
#ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

# Make the scatter plots
cs = ax1.scatter(x_locs, z,s = 2, c=m_chi, cmap="viridis",label='Longitudinal profile')

cbar = f.colorbar(cs)
cbar.set_label('Channel steepness ($k_{sn}$)')

ax1.set_xlabel("Distance from outlet ($m$)")
ax1.set_ylabel("Elevation (m)")

#ax2.set_xlabel("Distance from outlet ($m$)")
#ax2.set_ylabel("Channel steepness ($k_{sn}$)")

plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08645790005347911, 'fraction_numerical': 0.02038148434250401, 'mean_word_length': 4.878449179182676, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 59, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '240313', 'license': 'MIT', 'n_tokens_mistral': 5766, 'n_tokens_neox': 5401, 'n_words': 2295, 'path': 'Southern_spain/Lesson_02_looking_at_the_channel_profile.ipynb', 'repo_name': 'simon-m-mudd/smm_teaching_notebooks', 'size': '16877'}