Skip to content

Commit

Permalink
Merge pull request #12 from josiearcuri/bendtracking
Browse files Browse the repository at this point in the history
Bend Tracking
josiearcuri authored Mar 30, 2021
2 parents 0adabc9 + 3ff83dd commit 2c876d4
Showing 3 changed files with 239 additions and 11 deletions.
178 changes: 167 additions & 11 deletions HKplus/HKplus.py
Original file line number Diff line number Diff line change
@@ -43,7 +43,7 @@ def update_progress(progress, start_time):

class Channel:
"""class for Channel objects"""
def __init__(self,x,y,W,D):
def __init__(self,x,y,W,D,MR):
"""initialize Channel object
x, y, z - coordinates of centerline
W - channel width
@@ -52,6 +52,7 @@ def __init__(self,x,y,W,D):
self.y = y
self.W = W
self.D = D
self.MR = MR



@@ -122,13 +123,14 @@ def migrate_years(self,nit,saved_ts,deltas,pad, crdist,Cf,kl,dt,dens=1000):
omega = -1.0 # constant in curvature calculation (Howard and Knutson, 1984)
gamma = 2.5 # from Ikeda et al., 1981 and Howard and Knutson, 1984
ne = np.zeros_like(x) #array to keep track of nonlocal effects

MR = 0
DS = 0
for itn in range(nit): # main loop
update_progress(itn/nit, start_time)
ne = update_nonlocal_effects(ne, s, self.decay_rate, self.bump_scale, cut_dist, cut_len) #update array of ne with last itn's cutoff(s) and decay old ne
klarray = nominal_rate(kl, ne)## compute array of nominal migration rate in m/s with nonlocal effects accounted for
curv = compute_curvature(x,y)#get curvature of bends before cutoffs happen
x, y = migrate_one_step(x,y,W,klarray,dt,k,Cf,D,pad,omega,gamma)
x, y, R1 = migrate_one_step(x,y,W,klarray,dt,k,Cf,D,pad,omega,gamma)
x,y,xc,yc,cut_dist, cut_len, ind1, ind2 = cut_off_cutoffs(x,y,s,crdist,deltas) # find and execute cutoffs
x,y,dx,dy,ds,s = resample_centerline(x,y,deltas) # resample centerline

@@ -142,7 +144,7 @@ def migrate_years(self,nit,saved_ts,deltas,pad, crdist,Cf,kl,dt,dens=1000):
self.cutoffs.append(cutoff)
# saving centerlines:
if np.mod(itn,saved_ts)==0 or itn == nit-1:
channel = Channel(x,y,W,D) # create channel object, save year
channel = Channel(x,y,W,D, MR) # create channel object, save year
self.cl_times.append(last_cl_time+(itn)*dt/(365*24*60*60.0))
self.channels.append(channel)
self.sinuosity.append(Sin)
@@ -185,7 +187,8 @@ def migrate_cuts(self,saved_ts,deltas,pad, crdist,Cf,kl,dt,dens=1000):
ne = update_nonlocal_effects(ne, s, self.decay_rate, self.bump_scale, cut_dist, cut_len) #update array of ne with last itn's cutoff(s) and decay old ne
curv = compute_curvature(x,y)
klarray = nominal_rate(kl, ne)## compute array of nominal migration rate in m/s with nonlocal effects accounted for
x, y = migrate_one_step(x,y,W,klarray,dt,k,Cf,D,pad,omega,gamma)
x, y, R1 = migrate_one_step(x,y,W,klarray,dt,k,Cf,D,pad,omega,gamma)
MR, DS = segmented_MR(curv, R1, ds)
x,y,xc,yc,cut_dist, cut_len,ind1, ind2 = cut_off_cutoffs(x,y,s,crdist,deltas) # find and execute cutoffs
x,y,dx,dy,ds,s = resample_centerline(x,y,deltas) # resample centerline
Sin = get_sinuosity(x,s)
@@ -199,10 +202,68 @@ def migrate_cuts(self,saved_ts,deltas,pad, crdist,Cf,kl,dt,dens=1000):

# saving centerlines:
if np.mod(itn,saved_ts)==0 or len(self.cutoffs)>=self.cut_thresh:
channel = Channel(x,y,W,D) # create channel object, save year
channel = Channel(x,y,W,D,MR, DS) # create channel object, save year
self.cl_times.append(last_cl_time+(itn)*dt/(365*24*60*60.0))
self.channels.append(channel)
self.sinuosity.append(Sin)

def migrate_bendtracking(self,saved_ts,deltas,pad, crdist,Cf,kl,dt,dens=1000):
"""function for computing migration rates along channel centerlines and moving them, limited by number of cutoffs the channel experiences
inputs:
saved_ts - which time steps will be saved
deltas - distance between nodes on centerline
pad - padding for upstream bc (number of nodepoints along centerline)
crdist - threshold distance at which cutoffs occur
Cf - dimensionless Chezy friction factor
kl - migration rate constant (m/s)
dt - time step (s)"""
start_time = time.time()

channel = self.channels[-1] # first channel is the same as last channel of input
x = channel.x; y = channel.y; W = channel.W; D = channel.D;

k = 1.0 # constant in HK equation
xc = [] # initialize cutoff coordinates
yc = []
cut_dist = []# initialize cutoff distance ds array
cut_len = []# initialize cutoff length removal array
# determine age of last channel:
if len(self.cl_times)>0:
last_cl_time = self.cl_times[-1]
else:
last_cl_time = 0
dx, dy, ds, s = compute_derivatives(x,y)
omega = -1.0 # constant in curvature calculation (Howard and Knutson, 1984)
gamma = 2.5 # from Ikeda et al., 1981 and Howard and Knutson, 1984
ne = np.zeros_like(x) #array to keep track of nonlocal effects
ymax = self.bump_scale*kl*2
itn = 0

while len(self.cutoffs)<self.cut_thresh: # main loop
itn = itn+1
update_progress(len(self.cutoffs)/self.cut_thresh, start_time)
ne = update_nonlocal_effects(ne, s, self.decay_rate, self.bump_scale, cut_dist, cut_len) #update array of ne with last itn's cutoff(s) and decay old ne
curv = compute_curvature(x,y)
klarray = nominal_rate(kl, ne)## compute array of nominal migration rate in m/s with nonlocal effects accounted for
x, y, MR = migrate_one_step(x,y,W,klarray,dt,k,Cf,D,pad,omega,gamma)

x,y,xc,yc,cut_dist, cut_len,ind1, ind2 = cut_off_cutoffs(x,y,s,crdist,deltas) # find and execute cutoffs
x,y,dx,dy,ds,s = resample_centerline(x,y,deltas) # resample centerline

if len(xc)>0: # save cutoff data
rad = get_radii(curv, ind1, ind2, W)
cutoff = Cutoff(xc,yc,W,cut_dist,last_cl_time+(itn)*dt/(365*24*60*60.0), cut_len, rad) # create cutoff object
#keep track of year cutoff occurs, where it occurs, and save an object.
self.cutoff_times.append(last_cl_time+(itn)*dt/(365*24*60*60.0))
self.cutoff_dists.append(cut_dist)
self.cutoffs.append(cutoff)
# saving centerlines:
if np.mod(itn,saved_ts)==0 or len(self.cutoffs)>=self.cut_thresh:

channel = Channel(x,y,W,D,MR) # create channel object, save year
self.cl_times.append(last_cl_time+(itn)*dt/(365*24*60*60.0))
self.channels.append(channel)


def plot_channels(self):
cot = np.array(self.cutoff_times)
@@ -276,6 +337,19 @@ def cutoff_distributions(self, year, filepath, mode):
newcuts = cuts.to_csv(filepath+mode+str(len(cuts['time']))+"_cutoffs_distribution.csv", index_label = "Cutoff")
plot_cuts(cuts,self.channels[-1].W, filepath)
return cuts
def MR_time(self, filepath):
MR = [[bend for bend in i.MR] for i in self.channels[1:]]
clt = np.array(self.cl_times[1:])

MRdf= pd.DataFrame(MR).dropna(axis=1, how = 'all')
print(MRdf.head())
MRdf.to_csv(filepath, index_label = "Cutoff")
MRdf = pd.read_csv(filepath, sep = ',', index_col = 0)
plot_segmented_MR(MRdf)




def plot_cuts(cuts,W, filepath):
fig = plt.figure(figsize = (5,5))
plt.rcParams.update({'font.size': 10})
@@ -309,13 +383,16 @@ def migrate_one_step(x,y,W,klarray,dt,k,Cf,D,pad,omega,gamma):
R0 = klarray*curv #nominal migration rate with local curvature
alpha = k*2*Cf/D # exponent for convolution function G
R1 = compute_migration_rate(pad,len(x),ds,alpha,omega,gamma,R0)

# calculate bend-by-bend migration rate
MR = segmented_MR(curv, R1*dt, s)
# calculate new centerline coordinates:
dy_ds = dy/ds
dx_ds = dx/ds
# move x and y coordinates:
x = x + R1*dy_ds*dt
y = y - R1*dx_ds*dt
return x,y
return x,y, MR

def generate_initial_channel(W,D,deltas,pad):
"""generate straight Channel object with some noise added that can serve
@@ -333,7 +410,9 @@ def generate_initial_channel(W,D,deltas,pad):
x = np.linspace(0, cl_length+(2*pad1)*deltas, int(cl_length/deltas+(2*pad1))+1) # x coordinate
y = 10.0 * (2*np.random.random_sample(int(cl_length/deltas)+1,)-1)
y = np.hstack((np.zeros((pad1),),y,np.zeros((pad1),))) # y coordinate
return Channel(x,y,W,D)
MR = np.zeros_like(x)

return Channel(x,y,W,D, MR)

def load_initial_channel(filepath, W, D, deltas):
"""generate initial channel from centerline csv that can serve
@@ -345,7 +424,9 @@ def load_initial_channel(filepath, W, D, deltas):
df = pd.read_csv(filepath, sep = ',', header=None).values
x = df[:,0]
y = df[:,1]
return Channel(x,y,W,D)
MR = np.zeros(int(len(x)/30))

return Channel(x,y,W,D,MR)
def generate_channel_from_file(filelist, D_in= 10, smooth_factor=.25, matlab_corr= -1):
"""function for creating a MeanderPy Channel object from an externally-sourced centerline in .csv file format.
inputs:
@@ -373,7 +454,7 @@ def generate_channel_from_file(filelist, D_in= 10, smooth_factor=.25, matlab_cor

#average over widths to get a reach-constant width scalar
W = np.mean(varlist[1][:,0])*30
print(W)

## water depth scalar#
D = D_in
# Linear length along the line, add a zero for first point:
@@ -389,7 +470,9 @@ def generate_channel_from_file(filelist, D_in= 10, smooth_factor=.25, matlab_cor

## z-dim array, interpolated with constant slope along points of centerline. assumes centerline points are equidistantly placed along original centerline.
#deltas = round(distance[-1]/(len(points_fitted[0])-1))
return Channel(points_fitted[0],points_fitted[1],W,D)
MR = np.zeros_like(point_fitted[0])

return Channel(points_fitted[0],points_fitted[1],W,D, MR)

@numba.jit(nopython=True)
def compute_migration_rate(pad,ns,ds,alpha,omega,gamma,R0):
@@ -627,5 +710,78 @@ def scatter_hist(x, y, ax, ax_histx, ax_histy):
ybins = np.arange(0, ylim + ybinwidth, ybinwidth)
ax_histx.hist(x, bins=xbins)
ax_histy.hist(y, bins=ybins, orientation='horizontal')
def segmented_MR(curv, R1, s, n=90):
"""
approximate a bend-by-bend nth percentile lateral migration rate.
Inputs:
curve: array of curvature for every node along the centerline.
R1: array of already computed migration distances for every node along the centerline.
ds: array of cumulative distance downstream between nodes.
n: percentile
Output:
MR: array of nth percentile migration rate for each segment.
upstream: distance downstream of each segment start
downstream: distance downstream of each segment end
"""
R1 = np.array(R1)
#where curvature changes direction, =1
nodes = np.array([(curv[i-1]*curv[i])< 0 for i in range(1,len(curv)-1)])
idx = np.where(nodes==1)[0]
idx = [idx[i] for i in range(len(idx)-1) if idx[i+1]-idx[i] > 5]


MR = [max(np.abs(R1[idx[i]:idx[i+1]])) for i in range(len(idx)-1)]

return MR
def moving_average(matrix, window):
"""
averages migration rates on each bend over a set time window upwind.
inputs
matrix: n years by m bend migration rate numpy array
window: how many indices to average over, as years in the past
output
mid: moving average migration rate over window years for each bend
"""

years = len(matrix[:,0])
bends = np.min([np.count_nonzero(~np.isnan(matrix[i, :])) for i in range(years)])
print(years)
print(bends)
mid = np.zeros(shape = (years, bends))

for year in range(years):
if year < window:
mid[year, :] = np.mean(matrix[0:year, :bends], axis = 0)
else:
mid[year, :] = np.mean(matrix[year-window:year, :bends], axis = 0)

return mid

def plot_segmented_MR(MR):

#mean_mr_per_year = np.nanmean(MR, axis = 1)
# mean_mr_per_bend = np.nanmean(MR, axis = 0)

fig, axs = plt.subplots(1,2, figsize=(20,10), sharey= True)
#heatmap = ax1.imshow(MR, cmap = 'gist_heat')
heatmap = axs[1].imshow(MR, cmap = 'cividis', vmin = 0, vmax =np.nanmax(MR), aspect = 'auto', origin = 'lower')
axs[1].set_xlabel('distance downstream (bend #)')

cb = fig.colorbar(heatmap, ax=axs[1])
cb.set_label("maximum migration rate (m/yr)")
mid = moving_average(MR, 2)
#min
axs[0].plot(mid, range(len(mid)), alpha=.2, c = 'k')
#mean
axs[0].plot(np.nanmean(mid, axis = 1), range(len(mid)), alpha=1, c = 'r')
axs[0].set_xlim((0,np.nanmax(mid)))
axs[0].set_ylim((0,len(MR[:, 0])))
axs[0].set_ylabel('time (yr)')
axs[0].set_xlabel('max mr along bend (m/yr)')

return fig

56 changes: 56 additions & 0 deletions TrackMigrationThroughTime.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
"""
This script uses the HKplus functions to load an existing centerline, run it for a certain number of years, then analyize how migration rate changes through time on a bend-by-bend resultion
"""
import math
import matplotlib.pyplot as plt
import HKplus as hkp
import numpy as np
import pandas as pd

#Set Variables for centerline and curvature calculation

D=3.4
W = 100 #constant width
deltas = W//2; #spacing of nodes along centerline
nit = 1000 # number of iterations
Cf = 0.005 # dimensionless Chezy friction factor
kl =22.7/(365*24*60*60.0) # migration rate constant (m/s)
dt = .5*365*24*60*60.0 # time step (s)
pad= 200 # dont change
saved_ts = 2 # which time steps centerline will be saved at
crdist = 4*W # how close banks get before cutoff in m


#Set Variables fro nonlocal efects
decay_rate = dt/(5*(365*24*60*60.0)); #ranges between 1/3 to 1/10, to be developed
bump_scale = 0 #to multiple kl by,amplitude of ne bump, range between 1 and 4, set to 0 for no nonlocal effects
cut_thresh = 10 #how many cutoffs to simulate, arbitrary if running for time

#Set Result Directory
result_dir = "sample_results/case2/7/" ##change this to wherever you want to save your results

#Load Existing Centerline

filepath =result_dir + "InitialCL_result.csv"
ch= hkp.load_initial_channel(filepath, W, D, deltas)

#Initialize Channel Belt for Migration
chb = hkp.ChannelBelt(channels=[ch],cutoffs=[],cl_times=[0.0],cutoff_times=[], cutoff_dists = [], decay_rate = decay_rate, bump_scale = bump_scale, cut_thresh = cut_thresh, sinuosity = [3])

#Plot initial channel
chb.plot_channels()
plt.title("Initial Centerline")
plt.show()

#Migrate
chb.migrate_bendtracking(saved_ts,deltas,pad,crdist,Cf,kl,dt)

#Plot resulting centerline
chb.plot_channels()
plt.title(str(int(nit*dt/(365*24*60*60.0)))+ " years at "+ str(kl*(365*24*60*60.0))+ "m/yr")
plt.savefig(result_dir+"channelfrombendbybend.png", dpi = 500)
plt.close()

chb.MR_time(result_dir+"bendbybendmr.csv")
plt.show()
16 changes: 16 additions & 0 deletions plotmigrationtroughtime.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@

import math
import matplotlib.pyplot as plt
import HKplus as hkp
import numpy as np
import pandas as pd



path = "sample_results/case2/17/"
name = "bendbybendmr.csv"
MR = pd.read_csv(path+name, sep = ',', index_col = 0)

hkp.plot_segmented_MR(MR.to_numpy())
plt.savefig(path+"bendbybend.png")
plt.show()

0 comments on commit 2c876d4

Please sign in to comment.