python: Add slope.py.

* python/slope.py: New file.

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.")