aboutsummaryrefslogtreecommitdiff
path: root/python
diff options
context:
space:
mode:
Diffstat (limited to 'python')
-rwxr-xr-xpython/slope.py206
1 files changed, 206 insertions, 0 deletions
diff --git a/python/slope.py b/python/slope.py
new file mode 100755
index 0000000..c53b811
--- /dev/null
+++ b/python/slope.py
@@ -0,0 +1,206 @@
+#!/usr/bin/env python3
+
+import matplotlib.colors as colors
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
+with open('../topography/LiDAR_DTM_cropped.xyz', "r") as data:
+ while True:
+ p = data.tell()
+ line = data.readline()
+ col1, col2, col3 = line.strip().split()
+ if col2 == 'ncols': NXCELL = int(col3) #number of cells in x direction
+ if col2 == 'nrows': NYCELL = int(col3) #number of cells in y direction
+ if not line.startswith('#'):
+ data.seek(p) # go back one line
+ break
+
+ x, y, z = np.loadtxt(data, delimiter=' ', unpack=True)
+
+# first reshape to 2-D array then rotate by ninety degrees and flip in
+# vertical direction to conform to FullSWOF indexing convention
+x_co = np.flipud(np.rot90(np.reshape(x, (NXCELL,NYCELL)))) # x co-ordinates
+y_co = np.flipud(np.rot90(np.reshape(y, (NXCELL,NYCELL)))) # y co-ordinates
+elev = np.flipud(np.rot90(np.reshape(z, (NXCELL,NYCELL)))) # elevation map
+
+# calculate cell size in x and y directions
+DX = x_co[0,1] - x_co[0,0]
+DY = y_co[1,0] - y_co[0,0]
+#print(DX)
+
+#print(y_co[0])
+#plt.plot(x,z, label='elevation')
+#plt.xlabel('x / m')
+#plt.ylabel('z / m')
+
+# initialize boundary cells
+bb = np.zeros((1, NXCELL)) # bottom boundary
+tb = np.zeros((1, NXCELL)) # top boundary
+lb = np.zeros((1, NYCELL + 2)) # left boundary
+rb = np.zeros((1, NYCELL + 2)) # right boundary
+# use forward/backward differencing to calculate boundary elevation values
+for i in range(0, NXCELL):
+ bb[0,i] = 2.0*elev[0,i] - elev[1,i]
+ tb[0,i] = 2.0*elev[NYCELL-1,i] - elev[NYCELL-2,i]
+
+# add boundary cells
+z2 = np.concatenate((np.concatenate((bb, elev)), tb))
+z3 = np.concatenate((np.concatenate((lb.T, z2), axis=1), rb.T), axis=1)
+#print(z2[NYCELL + 1])
+#print(z3.shape)
+
+# initialize gradient dz/dy
+dzdy = np.zeros((NYCELL, NXCELL))
+# use central differencing to calculate nodal values
+for i in range(1, NYCELL + 1):
+ for j in range(1, NXCELL + 1):
+ dzdy[i-1, j-1] = (z3[i+1, j] - z3[i-1, j])/(2.0*DY)
+# array index and co-ordinates are related by:
+# index = (co-ord - 0.25)*2
+#print("dzdy", dzdy[NYCELL-1, 280:301])
+#print("elev", elev[NYCELL-1, 280:301])
+
+xMarker = []
+yMarker = []
+
+xch = [] # /pixels
+ych = [] # /pixels
+zch = [] # /metres
+
+def detach_display():
+ fig, (ax1,ax2) = plt.subplots(1, 2, sharey='row')
+ #fig, ax1 = plt.subplots()
+
+ left = ax1.imshow(elev,
+ interpolation='nearest',
+ origin='lower',
+ cmap=plt.get_cmap('terrain_r'))
+ ax1.set_title('elevation / m')
+ ax1.set_xticks([0, 100, 200, 300, 400, 500])
+ ax1.plot(xMarker, yMarker, 'ro')
+ fig.colorbar(left, ax=ax1)
+
+ right = ax2.imshow(dzdy,
+ interpolation='nearest',
+ origin='lower',
+ norm=colors.Normalize(vmin=0, vmax=0.1))
+ ax2.set_title('slope')
+ ax2.set_xticks([0, 100, 200, 300, 400, 500])
+ fig.colorbar(right, ax=ax2)
+
+ # press a suitable key (one that is not bound already) to mark bank.
+ def on_key(event):
+ if event.inaxes is not None:
+ # print(np.rint(event.xdata).astype(int),
+ # np.rint(event.ydata).astype(int))
+ xMarker.append(np.rint(event.xdata).astype(int))
+ yMarker.append(np.rint(event.ydata).astype(int))
+ ax1.plot(event.xdata, event.ydata, 'ro')
+ fig.canvas.draw()
+ else:
+ print('Key pressed ouside axes bounds but inside plot window')
+ # TODO: check xMarker, yMarker have even number of items
+ # (all left bank and right bank markers are present).
+
+
+ #print(elev.min())
+ # Note: row-index is plotted vertically, column-index plotted
+ # horizontally.
+ #print(np.unravel_index(elev.argmin(), elev.shape))
+
+ cid = fig.canvas.callbacks.connect('key_press_event', on_key)
+ #plt.figure()
+ plt.show()
+ fig.canvas.callbacks.disconnect(cid)
+ #print(xMarker)
+
+def plot_curve():
+ # sort markers by y value
+ #yx = list(zip(yMarker, xMarker))
+ #yx.sort()
+ #x_sorted = [x for y, x in yx]
+ #y_sorted = [y for y, x in yx]
+ #print(x_sorted)
+ #print(y_sorted)
+
+ xch.clear()
+ ych.clear()
+ zch.clear()
+ for i in range(0, len(xMarker), 2):
+ length = int(np.hypot(xMarker[i+1]-xMarker[i],
+ yMarker[i+1]-yMarker[i]))
+ #print(length)
+ x = np.linspace(xMarker[i+1], xMarker[i], length)
+ y = np.linspace(yMarker[i+1], yMarker[i], length)
+ # Extract the values along the line. Note index order: row,
+ # column.
+ zi = elev[y.astype(np.intp), x.astype(np.intp)]
+ # identify minimum value
+ zch.append(zi.min())
+ ind = np.unravel_index(zi.argmin(), zi.shape)
+ xch.append(x[ind].astype(np.intp))
+ ych.append(y[ind].astype(np.intp))
+ #print('xch = ', xch[i//2], ' ych = ', ych[i//2], ' zch = ', zch[i//2])
+
+ # TODO: Define starting point in channel. Find nearest neighbour to
+ # starting point. Extract from xch, ych, zch. Create new list with
+ # extracted point. Find nearest neighbour to extracted point. Extract.
+ # Add to list. Repeat.
+
+ # Fit with polyfit
+ m, c = np.polyfit([499.75 - y*0.5 for y in ych], zch, 1)
+ print('gradient =', m, 'intercept =', c)
+
+ fig, ax = plt.subplots()
+ line1, = ax.plot([499.75 - y*0.5 for y in ych], zch)
+ line2, = ax.plot([499.75 - y*0.5 for y in ych],
+ [c + m*499.75 - m*y*0.5 for y in ych], '--')
+ ax.set_xlabel('Distance / m')
+ ax.set_ylabel('Height / m')
+ ax.set_title('Channel profile')
+ #plt.figure()
+ plt.show()
+ #fig.clear()
+ #plt.clf()
+ #plt.cla()
+ #plt.close()
+ #fig.canvas.draw()
+ #fig.canvas.flush_events()
+
+def save_xyz():
+ with open('1D.txt', 'w') as f:
+ for xitem, yitem, zitem in zip(reversed(xch),
+ reversed(ych),
+ reversed(zch)):
+ f.write('{:.2f} {:.2f} {:.4f}\n'.format(0.25+xitem*0.5,
+ 499.75-yitem*0.5,
+ zitem))
+
+
+EXIT_COMMAND = "q"
+
+# if os.fork():
+# # parent
+# pass
+# else:
+# # child
+# detach_display()
+
+#detach_display()
+
+while (True):
+ input_str = input('Locate markers (m), plot channel profile (p), save profile (s) or exit (q): ')
+ #print("input_str = {}".format(input_str))
+
+ if (input_str == EXIT_COMMAND):
+ print("End.")
+ break
+ elif(input_str == "m"):
+ detach_display()
+ elif(input_str == "p"):
+ plot_curve()
+ elif(input_str == "s"):
+ save_xyz()
+
+#print("End.")