In [1]:
import python_subgrid.wrapper
import python_subgrid.plotting
import matplotlib.cm
import matplotlib.colors
import matplotlib.pyplot as plt
import numpy as np
import logging
import shapely.ops
import scipy.spatial
import scipy.interpolate
In [2]:
subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/home/fedor/Checkouts/models/testcase_1d_levee/testcase_1d_levee.mdu")
#subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/home/fedor/Checkouts/models/delfland-model-voor-3di/hhdlipad.mdu")
python_subgrid.wrapper.logger.setLevel(logging.WARN)
In [3]:
subgrid.start()
In [44]:
grid = {}
for var in {'FlowElem_xcc', 'FlowElem_ycc', 'nod_type', 'xlevnd', 'ylevnd', 'FlowElemContour_x',
'FlowElemContour_y', 'dmax', 'dps', 'x0p', 'y0p', 'x1p', 'y1p',
'dxp', 'dyp',
'levnd0', 'levnd1',
'line'}:
grid[var] = subgrid.get_nd(var).copy()
quad_grid = python_subgrid.plotting.make_quad_grid(subgrid)
In [45]:
for i in range(20):
subgrid.update(-1)
vars = {}
for var in {'s1', 'q', 'vol1'}:
vars[var] = subgrid.get_nd(var, sliced=True).copy()
vars['wl'] = vars['s1'] - -grid['dmax']
colors = matplotlib.cm.Blues(matplotlib.colors.Normalize(vmin=vars['wl'].min(), vmax=vars['wl'].max())(vars['wl']))
In [46]:
x = grid['FlowElem_xcc'][grid['nod_type']==1]
y = grid['FlowElem_ycc'][grid['nod_type']==1]
leveesplits = np.where(np.diff(grid['levnd0']) != 1)[0]
xlev = grid['xlevnd']
ylev = grid['ylevnd']
points = np.r_[
np.c_[x,y],
np.c_[xlev, ylev]
]
In [46]:
In [47]:
lev_segment_start = np.split(grid['levnd0']-1, leveesplits+1)
lev_segment_end = np.split(grid['levnd1']-1, leveesplits+1)
lines = []
for (seg_begin, seg_end) in zip(lev_segment_start, lev_segment_end):
line_x = xlev[seg_begin]
line_y = ylev[seg_begin]
prodigal = (set(seg_end) - set(seg_begin)).pop()
line_x = np.append(line_x, xlev[prodigal])
line_y = np.append(line_y, ylev[prodigal])
lines.append([(line_x_i, line_y_i) for line_x_i, line_y_i in zip(line_x, line_y)])
levees = shapely.geometry.MultiLineString(lines)
In [48]:
for levee in levees:
plt.plot(levee.xy[0], levee.xy[1], '-')
In [49]:
tri0 = scipy.spatial.Delaunay(points)
tri1 = scipy.spatial.Delaunay(np.c_[x,y], incremental=True)
tri1.add_points(np.c_[xlev, ylev])
In [50]:
fig, axes =plt.subplots(1,2, figsize=(20,7))
axes[0].triplot(points[:,0], points[:,1], tri0.simplices.copy())
for levee in levees:
axes[0].plot(levee.xy[0], levee.xy[1], '-')
axes[1].triplot(points[:,0], points[:,1], tri1.simplices.copy())
for levee in levees:
axes[1].plot(levee.xy[0], levee.xy[1], '-')
In [51]:
import shapely.geometry
In [52]:
#levee = shapely.geometry.LineString(np.c_[xlev, ylev])
#levee_start = shapely.geometry.Point(xlev[0], ylev[0])
#levee_end = shapely.geometry.Point(xlev[-1], ylev[-1])
In [53]:
xc = grid['FlowElemContour_x']
yc = grid['FlowElemContour_y']
polys = [shapely.geometry.Polygon(np.c_[xc_i, yc_i]) for xc_i, yc_i in zip(xc, yc)]
points = [shapely.geometry.Point(x_i, y_i) for x_i, y_i in zip(x,y)]
In [54]:
boundaries = shapely.ops.unary_union([levee.boundary for levee in levees])
{boundaries.relate(poly) for poly in polys[::10]}
Out[54]:
In [55]:
intersects = [levees.intersects(poly) for poly in polys]
inside = np.array([not poly.disjoint(boundaries) for poly in polys])
idx = np.logical_and(intersects, np.logical_not(inside))
(splitidx,) = np.where(idx)
for i in splitidx:
poly.relate(levee), levee.relate(poly), i
In [56]:
splitidx
Out[56]:
In [57]:
def split(poly, line):
"""split a polygon by a line, for convex polys"""
edges = poly.exterior.union(line)
polys = shapely.ops.polygonize(edges)
for poly_split in polys:
if poly.contains(poly_split):
yield poly_split
In [58]:
# idx, polys, levee, points, grid,
def search(idx, polys, levees, points, grid):
"""search for waterlevel cells that intersect with the cell idx and contain a water level point"""
poly = polys[idx]
geoms = split(poly, levees)
for cell in geoms:
if cell.contains(points[idx]):
continue
line = grid['line'] -1
connected = set(line[:,(line == idx).any(0)].ravel()) - set([idx])
found = False
for i in connected:
for geom in split(polys[i], levees):
if geom.contains(points[i]):
if geom.relate(cell)[4] not in {'F', '0'}:
found = True
yield i, geom
if not found:
yield idx, cell
connections = {}
for idx in splitidx:
connections[idx] = list(search(idx, polys, levees, points, grid))
In [59]:
print connections[167]
for poly in list(split(polys[167], levees)):
plt.fill(*poly.exterior.xy)
In [66]:
import shapely.affinity
connectedpoints = []
connectedidx = []
for idx, links in connections.items():
poly = polys[idx]
geoms = split(poly, levees)
for geom in geoms:
geom = shapely.affinity.scale(geom, 0.9, 0.9)
pts = np.array(geom.exterior)
connectedpoints.append(pts)
if geom.contains(points[idx]):
connectedidx.extend([idx]*pts.shape[0])
else:
l = [link for link, cell in links if not cell.relate(geom)[4] in {'F', '0'}]
if not l:
l = links[0][0]
connectedidx.extend([l]*pts.shape[0])
connectedpoints = np.concatenate(connectedpoints, axis=0)
connectedidx = np.array(connectedidx)
print connectedpoints.shape, connectedidx.shape
otherpoints = np.array([[point.x, point.y] for i, point in enumerate(points) if i not in connections])
otheridx = np.array([i for i in range(len(points)) if i not in connections])
tripoints = np.r_[otherpoints, connectedpoints]
triidx = np.r_[otheridx, connectedidx]
L = scipy.interpolate.LinearNDInterpolator(tripoints, vars['s1'][triidx])
#L = scipy.interpolate.Rbf(tripoints[:,0], tripoints[:,1], vars['wl'][triidx])
X, Y = np.mgrid[grid['x0p']:grid['x1p']:grid['dxp'], grid['y0p']:grid['y1p']:grid['dyp']]
In [71]:
plt.imshow((L(X,Y).T) - -dps , cmap='Blues', vmin=0, vmax=1, origin='lower')
Out[71]:
In [74]:
tri = scipy.spatial.Delaunay(tripoints)
fig, axes= plt.subplots(1,2, figsize=(20,7), sharex=True ,sharey=True)
N = matplotlib.colors.Normalize(0,1)
cmap = matplotlib.cm.Blues
for xc_i, yc_i, wl in zip(xc, yc, vars['wl']):
alpha= 0.5
axes[0].fill(xc_i, yc_i, alpha=alpha, facecolor=cmap(N(wl)), edgecolor='green', linewidth=3)
axes[1].fill(xc_i, yc_i, alpha=alpha, facecolor='none', edgecolor='green', linewidth=3)
#axes[0].triplot(tripoints[:,0], tripoints[:,1], tri.simplices.copy(), alpha=0.3, color='black', linewidth=1)
for levee in levees:
axes[0].plot(levee.xy[0], levee.xy[1], '-')
axes[1].plot(levee.xy[0], levee.xy[1], '-')
axes[0].plot(x, y, 'k.', alpha=0.3)
origin= (grid['x0p'], grid['y0p'])
extent=(grid['x0p'], grid['x1p'], grid['y0p'], grid['y1p'])
#cmap(N(L(X,Y)))
im= axes[1].imshow((L(X,Y).T - - dps), extent=extent, vmax=1, vmin=0,cmap='Blues', interpolation='none')
axes[1].set_xlim(147250,148000)
axes[1].set_ylim(527000, 527500)
axes[1].triplot(tripoints[:,0], tripoints[:,1], tri.simplices.copy(), alpha=0.3)
for i, links in connections.items():
poly = polys[i]
#axes[0].fill(*poly.exterior.xy)
geoms = split(poly, levees)
for cell in geoms:
x1, y1 = cell.centroid.x, cell.centroid.y
if not cell.contains(points[i]):
l = [link for link, g in links if not g.relate(cell)[4] in {'F', '0'}][0]
if not l:
l = links[0][0]
x0, y0 = points[l].x, points[l].y
axes[0].text(x1, y1, l)
dx = x1 - x0
dy = y1 - y0
axes[0].arrow(x0, y0, dx, dy, head_width=10, length_includes_head=True, color='black', antialiased=True)
plt.colorbar(im, ax=axes[1] )
Out[74]:
In [42]:
for poly in split(polys[29], levees):
plt.fill(*poly.exterior.xy)
connections[29]
Out[42]:
In [36]:
dps = np.ma.masked_outside( grid['dps'][1:-1,1:-1], -100,100)
fig, ax = plt.subplots(figsize=(20,13))
#ax.imshow(vars['s1'][quad_grid] - -dps , cmap='Blues', vmin=0,
# origin=(grid['x0p'], grid['y0p']), extent=(grid['x1p']-grid['x0p'], grid['y1p']-grid['y0p']))
for xc_i, yc_i, color in zip(xc, yc, colors):
alpha= 0.5
ax.fill(xc_i, yc_i, alpha=alpha, facecolor=color)
for i, poly in enumerate(polys):
if i not in splitidx:
continue
for geom in split(poly, levee):
geom.type == 'Polygon'
pts = np.array(geom.exterior)
alpha = 1 if points[i].intersects(geom) else 0.3
ax.fill(pts[:,0], pts[:,1], alpha=alpha, facecolor=colors[i])
for geom in split(polys[splitidx[0]], levee):
pts = np.array(geom.exterior)
ax.fill(pts[:,0], pts[:,1])
for i, links in connections.items():
poly = polys[i]
geoms = split(poly, levee)
cell = [geom for geom in geoms if not geom.contains(points[i])][0]
x0, y0 = cell.centroid.x, cell.centroid.y
ax.text(x0, y0, i)
for l in links:
dx = points[l].x- x0
dy = points[l].y- y0
ax.arrow(x0, y0, dx, dy, head_width=10)
ax.axis('equal')
vars['wl'].min(), vars['wl'].max()
In [ ]:
list(split(polys[264], levee))
In [ ]:
poly.