In [1]:
%matplotlib inline
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
import shapely.geometry
In [2]:
subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/Users/baart_f/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 [4]:
# Get the variables needed or the levee update
grid = {}
for var in {'FlowElem_xcc', 'FlowElem_ycc', 'nod_type', 'FlowElemContour_x',
'FlowElemContour_y', 'dmax', 'dps', 'x0p', 'y0p', 'x1p', 'y1p',
'dxp', 'dyp',
'levnd0', 'levnd1', 'xleve', 'yleve', 'dlev', 'xlevnd', 'ylevnd'}:
value = subgrid.get_nd(var, sliced=True)
if value is not None:
grid[var] = value.copy()
quad_grid = python_subgrid.plotting.make_quad_grid(subgrid)
In [5]:
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 [5]:
In [6]:
x = grid['FlowElem_xcc'][grid['nod_type']==1]
y = grid['FlowElem_ycc'][grid['nod_type']==1]
In [7]:
# if we have a model with levees
if 'xlevnd' in grid:
# we want to have separate levees
# lookup the split points
leveesplits, = np.where(np.diff(grid['levnd0']) != 1)
xlev = grid['xlevnd']
ylev = grid['ylevnd']
points = np.r_[
np.c_[x,y],
np.c_[xlev, ylev]
]
# split up the levees, we now have a list of arrays with start end points
lev_segment_start = np.split(grid['levnd0']-1, leveesplits+1)
lev_segment_end = np.split(grid['levnd1']-1, leveesplits+1)
# reconstruct the separate lines
lines = []
# match up begin and end
for (seg_begin, seg_end) in zip(lev_segment_start, lev_segment_end):
line_x = xlev[seg_begin]
line_y = ylev[seg_begin]
# add the final point
prodigal = (set(seg_end) - set(seg_begin)).pop()
# it should be the last point
assert prodigal == seg_end[-1]
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)
levees.to_wkt()
Out[7]:
In [8]:
def render(subgrid, grid, quad_grid, patch):
fig, ax = plt.subplots(figsize=(20,13))
origin= (grid['x0p'], grid['y0p'])
extent=(grid['x0p'], grid['x1p'], grid['y0p'], grid['y1p'])
wl = subgrid.get_nd('s1', sliced=True)[quad_grid.filled(-1)]- -grid['dps']
vol_mask = (subgrid.get_nd('vol1', sliced=True)==0)[quad_grid.filled(-1)]
wl_masked = np.ma.masked_array(wl, vol_mask)
im = ax.imshow(wl_masked, cmap='Blues', vmin=0.0, vmax=5, origin=origin, extent=extent)
#im = ax.imshow(subgrid.get_nd('s1', sliced=True)[quad_grid.filled(-1)], cmap='Blues', vmin=-5, vmax=5, origin=origin, extent=extent)
xc = grid['FlowElemContour_x']
yc = grid['FlowElemContour_y']
verts = [np.c_[xc_i, yc_i] for xc_i, yc_i in zip(xc, yc)]
poly_collection = matplotlib.collections.PolyCollection(verts, closed=True, linewidth=0.1, transOffset=ax.transData, facecolor='none')
if 'xlevnd' in grid:
for levee in levees:
plt.plot(levee.xy[0], levee.xy[1], '-+')
sc = ax.scatter(grid['xleve'].mean(0), grid['yleve'].mean(0), c=grid['dlev'], s=3*(-grid['dlev']-10), edgecolor='none', vmin=-1,vmax=1)
plt.colorbar(sc, ax=ax)
plt.colorbar(im, ax=ax)
#trans = fig.dpi_scale_trans + matplotlib.transforms.Affine2D().scale(1.0/72.0)
#poly_collection.set_transform(trans) # the points to pixels transform
ax.add_collection(poly_collection)
ax.add_patch(patch)
ax.autoscale_view()
In [9]:
# let the simulation run for a bit
for i in range(20):
subgrid.update(-1)
# create an area where we want to update the levees
# assume that if you lower the rectangle, that the area south of it will flood
patch = matplotlib.patches.Rectangle([147790,527080], 70,50, alpha=0.3, facecolor='green')
# plot the model
render(subgrid, grid, quad_grid, patch)
In [10]:
# Convert the rectangle to a polygon
# We want it to work for any polygon (or at least the exterior)
def bbox2poly(bbox):
coords = [
bbox.min,
[bbox.max[0], bbox.min[1]],
bbox.max,
[bbox.min[0], bbox.max[1]]
]
return shapely.geometry.Polygon(shell=coords)
bbox = patch.get_bbox()
polygon = bbox2poly(bbox)
# print the geojson
shapely.geometry.mapping(polygon)
Out[10]:
In [11]:
# intersect the levees with the polygon
shapely.geometry.mapping(levees.intersection(polygon))
Out[11]:
In [12]:
import shapely.prepared
# create a prepared polygon for performance
polygon_prepared = shapely.prepared.prep(polygon)
In [13]:
# lookup the levees that are at least partially in the polygon
start_points = np.c_[grid['xleve'][0], grid['yleve'][0]]
end_points = np.c_[grid['xleve'][1], grid['yleve'][1]]
lines = []
for start, end in zip(start_points, end_points):
line = shapely.geometry.LineString([start, end])
lines.append(line)
# Use contains if you only want the levees that are fully contained
in_levee = [polygon_prepared.intersects(line) for line in lines]
levee_idx, = np.where(in_levee)
levee_idx
Out[13]:
In [14]:
# lower the levees
grid['dlev'][levee_idx] = -40
grid['dlev']
# update tables, nodes are ignored for now
subgrid.update_tables('dlev', levee_idx + 1)
# let the simulation run for a bit
for i in range(200):
subgrid.update(-1)
# plot the model
patch = matplotlib.patches.Rectangle([147790,527080], 70,50, alpha=0.3, facecolor='red')
render(subgrid, grid, quad_grid, patch)
In [15]:
# mark the point that we'll check for waterlevel
fig, ax = plt.subplots(figsize=(20,13))
#im = ax.imshow(subgrid.get_nd('s1', sliced=True)[quad_grid.filled(-1)], cmap='Blues', vmin=-5, vmax=5, origin=origin, extent=extent)
xc = grid['FlowElemContour_x']
yc = grid['FlowElemContour_y']
verts = [np.c_[xc_i, yc_i] for xc_i, yc_i in zip(xc, yc)]
poly_collection = matplotlib.collections.PolyCollection(verts, closed=True, linewidth=0.1, transOffset=ax.transData, facecolor='none')
#trans = fig.dpi_scale_trans + matplotlib.transforms.Affine2D().scale(1.0/72.0)
#poly_collection.set_transform(trans) # the points to pixels transform
ax.add_collection(poly_collection)
poly_collection = matplotlib.collections.PolyCollection([verts[231]], closed=True, linewidth=0.1, transOffset=ax.transData, facecolor='red')
ax.add_collection(poly_collection)
if 'xlevnd' in grid:
for levee in levees:
plt.plot(levee.xy[0], levee.xy[1], '-+')
sc = ax.scatter(grid['xleve'].mean(0), grid['yleve'].mean(0), c=grid['dlev'], s=3*(-grid['dlev']-10), edgecolor='none', vmin=-1,vmax=1)
plt.colorbar(sc, ax=ax)
In [15]: