1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
|
#!/usr/bin/env python3
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy as np
import math
import os
with open('./topography.txt', "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
xtp, ytp, ztp = np.loadtxt(data, delimiter=' ', unpack=True)
xmax = (xtp[0]+xtp[-1]) # domain extent in x-direction
ymax = (ytp[0]+ytp[-1]) # domain extent in y-direction
NXCELL = int(math.sqrt(len(xtp)*xmax/ymax)) # number of cells in x-direction
NYCELL = int(len(xtp)/NXCELL) # number of cells in y-direction
# 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(xtp, (NXCELL,NYCELL)))) # x co-ordinates
y_co = np.flipud(np.rot90(np.reshape(ytp, (NXCELL,NYCELL)))) # y co-ordinates
elev = np.flipud(np.rot90(np.reshape(ztp, (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)
# Note: dz/dy index and co-ordinates are related by:
# index = int(co-ord/DX - 0.5)
#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',
cmap=plt.get_cmap('rainbow'))
# norm=colors.Normalize(vmin=0, vmax=0.1))
ax2.set_title('slope ($dz/dy$)')
# 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(MY, PX):
"""
Fit straight line to channel elevation data. Pixel size PX is
used to convert pixel indices to co-ordinate values. The maximum
co-ordinate value in the y-direction is MY.
"""
# 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([MY - y*PX for y in ych], zch, 1)
print('gradient =', m, 'intercept =', c)
print('MY =', MY)
fig, ax = plt.subplots()
line1, = ax.plot([MY - y*PX for y in ych], zch)
line2, = ax.plot([MY - y*PX for y in ych],
[c + m*MY - m*y*PX 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(MY, PX):
"""
Write channel elevation data to file. Pixel size PX is used to
convert pixel indices to co-ordinate values. The maximum
co-ordinate value in the y-direction is MY.
"""
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(PX/2.0+xitem*PX,
MY-yitem*PX,
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(ytp[-1], DY)
elif(input_str == "s"):
save_xyz(ytp[-1], DY)
#print("End.")
|
