In [1]:
import python_subgrid.wrapper
import matplotlib.cm
import matplotlib.colors
import matplotlib.pyplot as plt
import numpy as np
In [2]:
subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/home/fedor/Checkouts/models/kaapstad_centrum/kaapstad_centrum.mdu")
In [3]:
subgrid.start()
INFO:python_subgrid.wrapper:Using subgrid fortran library /home/fedor/.local/lib/libsubgrid.so
INFO:python_subgrid.wrapper:Loading library from path /home/fedor/.local/lib/libsubgrid.so
INFO:python_subgrid.wrapper:Current version: Deltares, 3Di - subgrid Version 0.1.1.979M, Jan 3 2014, 21:18:50
INFO:python_subgrid.wrapper:Startup complete.
INFO:python_subgrid.wrapper:Loading model /home/fedor/Checkouts/models/kaapstad_centrum/kaapstad_centrum.mdu in directory /home/fedor/Checkouts/models/kaapstad_centrum
INFO:python_subgrid.wrapper:Start loadModel for '/home/fedor/Checkouts/models/kaapstad_centrum/kaapstad_centrum.mdu'.
WARNING:python_subgrid.wrapper:Ignored 'Teta' in MDU file, please use 'IntegrationMethod' instead.
INFO:python_subgrid.wrapper:Wrong NumLayers. Ground water calculations switched off.
INFO:python_subgrid.wrapper:Reading bathmetry file ' subgrid/kaapstad_centrum.tif '...
DEBUG:python_subgrid.wrapper: 4 4.39809990 dps(0:imax+1, 0:jmax+1)
INFO:python_subgrid.wrapper:Completed reading bathymetry file.
INFO:python_subgrid.wrapper:Opened file : polygons/verfijning.pol
INFO:python_subgrid.wrapper:Start initialization of grid administration.
DEBUG:python_subgrid.wrapper: 2034 2030.04309 ip,jp
DEBUG:python_subgrid.wrapper: 2035 1.05267596 ls(k)
DEBUG:python_subgrid.wrapper: 2035 0.264196008 ls(k)
DEBUG:python_subgrid.wrapper: 2035 6.65640011E-02 ls(k)
DEBUG:python_subgrid.wrapper: 2035 1.68999992E-02 ls(k)
DEBUG:python_subgrid.wrapper: 2035 4.35600011E-03 ls(k)
DEBUG:python_subgrid.wrapper: 2035 1.15599995E-03 ls(k)
INFO:python_subgrid.wrapper: Internal nodes = 54111
INFO:python_subgrid.wrapper: Boundary nodes = 0
DEBUG:python_subgrid.wrapper: 2035 5.41139990E-02 nodk,nodm,nodn
DEBUG:python_subgrid.wrapper: 2035 5.41110002E-02 ks(nodall)
DEBUG:python_subgrid.wrapper: 2036 0.108222000 xz,yz
INFO:python_subgrid.wrapper: Internal 1D nodes = 0
INFO:python_subgrid.wrapper: Boundary 1D nodes = 0
DEBUG:python_subgrid.wrapper: 2036 0.216260001 line
DEBUG:python_subgrid.wrapper: 2036 0.432520002 lik,lim,lin,kcu
DEBUG:python_subgrid.wrapper: 2037 0.432895988 nru,nrv
INFO:python_subgrid.wrapper:Initializing boundaries took: 0.7000D-01 seconds.
DEBUG:python_subgrid.wrapper: 2037 0.108130001 flod
DEBUG:python_subgrid.wrapper: 2037 0.108130001 flou
DEBUG:python_subgrid.wrapper: 2037 0.324389994 dpumin
DEBUG:python_subgrid.wrapper: 2038 0.324389994 dpumax
ERROR:python_subgrid.wrapper:No network found, cross sections are ignored.
INFO:python_subgrid.wrapper:Grid Admin initialisation took: 0.3980D+00 seconds.
INFO:python_subgrid.wrapper:Completed initialization of grid administration.
DEBUG:python_subgrid.wrapper: 2038 5.41110002E-02 dmax(1:nodall)
DEBUG:python_subgrid.wrapper: 2038 5.41110002E-02 dmin(1:nodall)
DEBUG:python_subgrid.wrapper: 2038 5.41110002E-02 sumax(1:nodtot)
DEBUG:python_subgrid.wrapper: 2038 0.216444001 dqmin(1:nodtot)
DEBUG:python_subgrid.wrapper: 2038 0.216444001 dqmax(1:nodtot)
DEBUG:python_subgrid.wrapper: 2038 0.216444001 suqmax(1:nodtot)
DEBUG:python_subgrid.wrapper: 2039 0.216444001 ksq(1:nodtot)
DEBUG:python_subgrid.wrapper: 2039 0.162336007 s1(0:nodall)
DEBUG:python_subgrid.wrapper: 2039 0.162336007 si(0:nodall)
DEBUG:python_subgrid.wrapper: 2039 0.324393004 u1(0:linall)
DEBUG:python_subgrid.wrapper: 2039 0.108130999 kf(0:linall)
DEBUG:python_subgrid.wrapper: 2039 5.41119985E-02 vol1(0:nodall)
DEBUG:python_subgrid.wrapper: 2039 0.108223997 vol2(0:nodall)
DEBUG:python_subgrid.wrapper: 2040 0.108223997 su(0:nodall)
DEBUG:python_subgrid.wrapper: 2040 0.108223997 volErr(0:nodall)
DEBUG:python_subgrid.wrapper: 2040 0.108130001 conv(linall)
DEBUG:python_subgrid.wrapper: 2040 0.216447994 vq1(0:4*nodall)
DEBUG:python_subgrid.wrapper: 2040 5.41119985E-02 implis(0:nodall)
DEBUG:python_subgrid.wrapper: 2040 0.270559996 aii,bi,d0,a00,b0(0:nodall)
DEBUG:python_subgrid.wrapper: 2041 0.216447994 vol0(0:nodall),vq0(0:4*nodall)
DEBUG:python_subgrid.wrapper: 2041 0.865779996 vqt,denom(0:4*nodall),fbq(0:4*nodall),umq(0:4*nodall)
DEBUG:python_subgrid.wrapper: 2042 0.324669003 qs,qs1,qs2(0:2*nodall)
DEBUG:python_subgrid.wrapper: 2042 0.324393004 u0,uh(0:linall)
DEBUG:python_subgrid.wrapper: 2042 5.41119985E-02 uc(0:nodall,2)
DEBUG:python_subgrid.wrapper: 2042 5.41119985E-02 vnorm(0:nodall)
DEBUG:python_subgrid.wrapper: 2043 0.540650010 cu,ru(linall)
DEBUG:python_subgrid.wrapper: 2043 0.324393004 au1(linall,0:2)
DEBUG:python_subgrid.wrapper: 2043 0.324393004 au0(linall,0:2)
DEBUG:python_subgrid.wrapper: 2044 0.216261998 as0(linall,1:2)
DEBUG:python_subgrid.wrapper: 2044 0.216261998 as1(linall,1:2)
DEBUG:python_subgrid.wrapper: 2044 0.108130001 cz(linall)
DEBUG:python_subgrid.wrapper: 2044 0.216261998 qh1,qh2(0:linall)
DEBUG:python_subgrid.wrapper: 2044 0.216261998 ade,adi(0:linall)
INFO:python_subgrid.wrapper:Calculation of bathymetry statistics took: 0.4266D+01 seconds.
INFO:python_subgrid.wrapper:Setting initial water levels from WaterLevelFile took: 0.0000D+00 seconds.
INFO:python_subgrid.wrapper:Initial floodfill from FloodFill file took: 0.0000D+00 seconds.
INFO:python_subgrid.wrapper:loadModel complete.
In [4]:
subgrid.update(-1)
INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init= 3
DEBUG:python_subgrid.wrapper: 2045 0.216260001 tables
DEBUG:python_subgrid.wrapper: 2045 0.216444001 tables
DEBUG:python_subgrid.wrapper: 2045 0.216444001 tables
DEBUG:python_subgrid.wrapper: 2045 0.216444001 tables
DEBUG:python_subgrid.wrapper: 2045 0.216260001 tables
DEBUG:python_subgrid.wrapper: 2046 0.216444001 tables
DEBUG:python_subgrid.wrapper: 2046 0.216444001 tables
DEBUG:python_subgrid.wrapper: 2046 5.41110002E-02 tables
DEBUG:python_subgrid.wrapper: 2046 5.41110002E-02 tables
DEBUG:python_subgrid.wrapper: 2065 19.0422955 autb
DEBUG:python_subgrid.wrapper: 2088 23.2193604 vltb
DEBUG:python_subgrid.wrapper: 2143 54.8972969 vq1tb,dentb
INFO:python_subgrid.wrapper:Initialization of Flow Values took: 0.4353D+01 seconds.
DEBUG:python_subgrid.wrapper: 2143 5.41110002E-02 ij(nodtot)
DEBUG:python_subgrid.wrapper: 2143 5.41110002E-02 noel(nodtot)
DEBUG:python_subgrid.wrapper: 2143 5.41110002E-02 noddgr(nodtot)
DEBUG:python_subgrid.wrapper: 2143 5.41110002E-02 nodstk(nodtot)
DEBUG:python_subgrid.wrapper: 2143 5.41110002E-02 row(nodtot)
DEBUG:python_subgrid.wrapper: 2143 0.108130001 ijl(lintot)
DEBUG:python_subgrid.wrapper: 2144 5.41110002E-02 ia(nodtot)
INFO:python_subgrid.wrapper:reduce time: 0.871000
INFO:python_subgrid.wrapper: noactive= 54111 nogauss= 27652 nocg 26459
DEBUG:python_subgrid.wrapper: 2144 5.41110002E-02 pk(nodtot)
DEBUG:python_subgrid.wrapper: 2144 5.41110002E-02 rk(nodtot)
DEBUG:python_subgrid.wrapper: 2144 5.41110002E-02 zk(nodtot)
DEBUG:python_subgrid.wrapper: 2144 5.41110002E-02 apk(nodtot)
DEBUG:python_subgrid.wrapper: 2144 0.220502004 aij(ijtot)
DEBUG:python_subgrid.wrapper: 2144 0.220502004 aij0(ijtot)
DEBUG:python_subgrid.wrapper: 2144 0.108130001 aijl(lintot)
DEBUG:python_subgrid.wrapper:Completed initFlowAdmin.
INFO:python_subgrid.wrapper:finding rain data for input
INFO:python_subgrid.wrapper:Timestep 1 number of iterations 1 epsmax 0.1066E-13 idthalf 1 np = 1
INFO:python_subgrid.wrapper:Timestep complete.
Out[4]:
0
In [6]:
python_subgrid.wrapper.DOCUMENTED_VARIABLES
Out[6]:
{u'FlowElemContour_x': u'List of x points forming flow element',
u'FlowElemContour_y': u'List of y points forming flow element',
u'FlowElem_xcc': u'Cell center x coordinate (pressure point) for all quadtree cells (i.e. nodes).',
u'FlowElem_ycc': u'Cell center y coordinate (pressure point) for all quadtree cells (i.e. nodes).',
u'FlowLink': u'Flow links/lines between coarse grid cells: line(L,1) = nod1, line(L,2) = nod2',
u'FlowLink_xu': u'Velocity x coordinate for 1d, 2d, including boundaries',
u'FlowLink_yu': u'Velocity y coordinate for 1d, 2d, including boundaries',
u'LandUse': u'percentage of a crop per node',
u'MaxInterception': u'max thickness of interception layer on fine base grid',
u'ade': u'advection',
u'adi': u'advection',
u'cropType': u'crop factors on fine grid',
u'culverts': u'Culvert',
u'dps': u'bathymetry pixel values on fine base grid',
u'ds1d': u'gridsize in 1d channels',
u'dsnop': u'Dummy missing value for bathymetry pixels, single precision.',
u'dt': u'delta t',
u'dtmax': u'maximum delta t',
u'dtmin': u'minumum delta t',
u'dx': u'gridsize for all coarse quadtree levels 1:kmax',
u'dxp': u'pixel dimensions',
u'dyp': u'pixel dimensions',
u'imax': u'number of pixels in x directions',
u'infiltrationRate': u'infiltration rate values on fine base grid',
u'ip': u'subgrid pixel numbers in coarse grid cells, see subroutine couplegrids(). E.g., ip(1:kmax,1:mmax(k),0:3).',
u'jmax': u'number of pixels in y directions',
u'jp': u'subgrid pixel numbers in coarse grid cells, see subroutine couplegrids(). E.g., ip(1:kmax,1:mmax(k),0:3).',
u'lg': u'Quadtree refinement level for each pixel (kmax=coarsest, 1=finest)',
u'lh': u'Quadtree node number to which each pixel belongs.',
u'link_branchid': u'Original link number',
u'link_chainage': u'Distance along the branch',
u'link_idx': u'Index in the u vector',
u'link_type': u'Type of link',
u'lu1dmx': u'number of u points per channel (for embedded: nr of 2D cell interfaces crossed by 1D channel)',
u'mbndry': u'coordinates of East boundaries',
u'mmax': u'quadtree grid administration',
u'nFlowElem1d': u'total number of nodal points in 1d without boundary points',
u'nFlowElem1dBounds': u'number of 1d boundary points',
u'nFlowElem2d': u'number of 2d nodes',
u'nFlowElem2dBounds': u'number of nodal points 2d boundary points',
u'nbndry': u'coordinates of North boundaries',
u'nmax': u'quadtree grid administration',
u'nod_type': u'Type of node',
u'nodk': u'inverse indirect adressing of coarse grid cells: nod = ls(k)$mn(m,n) --> m = nodm(nod)',
u'nodm': u'inverse indirect adressing of coarse grid cells: nod = ls(k)$mn(m,n) --> m = nodm(nod)',
u'nodn': u'inverse indirect adressing of coarse grid cells: nod = ls(k)$mn(m,n) --> m = nodm(nod)',
u'orifices': u'Orifice',
u'pumps': u'Pump',
u'q': u'discharge on full grid on current timestep',
u'qh1': u'discharge on half grid on current timestep',
u'qh2': u'discharge on half grid on next timestep',
u'qrain': u'rainfall',
u'rain': u'rainfall',
u'rootLength': u'root lengths on fine grid',
u's0': u'waterlevel at previous timestep',
u's1': u'waterlevel at current timestep',
u's2': u'waterlevel at next timestep',
u'sg': u'ground water level measured from "ground level" upwards (time dependent)',
u'si': u'thickness of interception layer per node',
u'si0': u'thickness of interception layer per node at previous timestep',
u'soilType': u'soil types on fine grid',
u't0': u'time at previous timestep',
u't1': u'time at current timestep',
u'tend': u'stop time of simulation',
u'tstart': u'start time of simulation',
u'u0': u'velocity at previous timestep',
u'u1': u'velocity at current timestep',
u'u2': u'velocity at next timestep',
u'weirs': u'Weir',
u'x0p': u'origin of pixel grid',
u'x1p': u'origin of pixel grid',
u'y0p': u'origin of pixel grid',
u'y1p': u'origin of pixel grid'}
Out[6]:
{u'FlowElemContour_x': u'List of x points forming flow element',
u'FlowElemContour_y': u'List of y points forming flow element',
u'FlowElem_xcc': u'Cell center x coordinate (pressure point) for all quadtree cells (i.e. nodes).',
u'FlowElem_ycc': u'Cell center y coordinate (pressure point) for all quadtree cells (i.e. nodes).',
u'FlowLink': u'Flow links/lines between coarse grid cells: line(L,1) = nod1, line(L,2) = nod2',
u'FlowLink_xu': u'Velocity x coordinate for 1d, 2d, including boundaries',
u'FlowLink_yu': u'Velocity y coordinate for 1d, 2d, including boundaries',
u'LandUse': u'percentage of a crop per node',
u'MaxInterception': u'max thickness of interception layer on fine base grid',
u'ade': u'advection',
u'adi': u'advection',
u'cropType': u'crop factors on fine grid',
u'culverts': u'Culvert',
u'dps': u'bathymetry pixel values on fine base grid',
u'ds1d': u'gridsize in 1d channels',
u'dsnop': u'Dummy missing value for bathymetry pixels, single precision.',
u'dt': u'delta t',
u'dtmax': u'maximum delta t',
u'dtmin': u'minumum delta t',
u'dx': u'gridsize for all coarse quadtree levels 1:kmax',
u'dxp': u'pixel dimensions',
u'dyp': u'pixel dimensions',
u'imax': u'number of pixels in x directions',
u'infiltrationRate': u'infiltration rate values on fine base grid',
u'ip': u'subgrid pixel numbers in coarse grid cells, see subroutine couplegrids(). E.g., ip(1:kmax,1:mmax(k),0:3).',
u'jmax': u'number of pixels in y directions',
u'jp': u'subgrid pixel numbers in coarse grid cells, see subroutine couplegrids(). E.g., ip(1:kmax,1:mmax(k),0:3).',
u'lg': u'Quadtree refinement level for each pixel (kmax=coarsest, 1=finest)',
u'lh': u'Quadtree node number to which each pixel belongs.',
u'link_branchid': u'Original link number',
u'link_chainage': u'Distance along the branch',
u'link_idx': u'Index in the u vector',
u'link_type': u'Type of link',
u'lu1dmx': u'number of u points per channel (for embedded: nr of 2D cell interfaces crossed by 1D channel)',
u'mbndry': u'coordinates of East boundaries',
u'mmax': u'quadtree grid administration',
u'nFlowElem1d': u'total number of nodal points in 1d without boundary points',
u'nFlowElem1dBounds': u'number of 1d boundary points',
u'nFlowElem2d': u'number of 2d nodes',
u'nFlowElem2dBounds': u'number of nodal points 2d boundary points',
u'nbndry': u'coordinates of North boundaries',
u'nmax': u'quadtree grid administration',
u'nod_type': u'Type of node',
u'nodk': u'inverse indirect adressing of coarse grid cells: nod = ls(k)$mn(m,n) --> m = nodm(nod)',
u'nodm': u'inverse indirect adressing of coarse grid cells: nod = ls(k)$mn(m,n) --> m = nodm(nod)',
u'nodn': u'inverse indirect adressing of coarse grid cells: nod = ls(k)$mn(m,n) --> m = nodm(nod)',
u'orifices': u'Orifice',
u'pumps': u'Pump',
u'q': u'discharge on full grid on current timestep',
u'qh1': u'discharge on half grid on current timestep',
u'qh2': u'discharge on half grid on next timestep',
u'qrain': u'rainfall',
u'rain': u'rainfall',
u'rootLength': u'root lengths on fine grid',
u's0': u'waterlevel at previous timestep',
u's1': u'waterlevel at current timestep',
u's2': u'waterlevel at next timestep',
u'sg': u'ground water level measured from "ground level" upwards (time dependent)',
u'si': u'thickness of interception layer per node',
u'si0': u'thickness of interception layer per node at previous timestep',
u'soilType': u'soil types on fine grid',
u't0': u'time at previous timestep',
u't1': u'time at current timestep',
u'tend': u'stop time of simulation',
u'tstart': u'start time of simulation',
u'u0': u'velocity at previous timestep',
u'u1': u'velocity at current timestep',
u'u2': u'velocity at next timestep',
u'weirs': u'Weir',
u'x0p': u'origin of pixel grid',
u'x1p': u'origin of pixel grid',
u'y0p': u'origin of pixel grid',
u'y1p': u'origin of pixel grid'}
In [11]:
vars = {'FlowElem_xcc', 'FlowElem_ycc', 'FlowElemContour_x', 'FlowElemContour_y', 'dx', 'nmax', 'mmax',
'mbndry', 'nbndry', 'ip', 'jp', 'nodm', 'nodn', 'nodk', 'nod_type'}
v = {var:subgrid.get_nd(var)
for var
in vars}
nodslice = slice(np.where(v['nod_type'] == 2)[0].min(), np.where(v['nod_type'] == 2)[0].max())
v['xcc'] = v['FlowElem_xcc'][nodslice]
v['ycc'] = v['FlowElem_ycc'][nodslice]
v['xc'] = v['FlowElemContour_x'][nodslice]
v['yc'] = v['FlowElemContour_y'][nodslice]
v['areas'] = v['dx'][nodslice]**2
In [194]:
# Create lookup index
n = v['nodn'][nodslice] # column lookup
m = v['nodm'][nodslice] # row lookup
k = v['nodk'][nodslice] # level lookup
# pixel boundary lookup
Y = v['ip'][:,m-1,k-1] # rows
X = v['jp'][:,n-1,k-1] # columns
minpx = 5 # is this fixed?
# slices
Xmin = (X[0,:] - 1)/minpx
Xmax = (X[3,:] + 1)/minpx
Ymin = (Y[0,:] - 1)/minpx
Ymax = (Y[3,:] + 1)/minpx
# total image size
mmax = Xmax.max()
nmax = Ymax.max()
# quad lookup
quad_grid = np.ma.empty((mmax,nmax), dtype='int32')
quad_grid.mask = True
slices = np.c_[Xmin, Xmax, Ymin, Ymax]
indices = np.arange(k.shape[0])
# fill single pixel quads separately
quad_grid[slices[k==1,0], slices[k==1,2]] = indices[k==1]
# fill the bigger quads
for i, (xmin, xmax, ymin, ymax) in zip(indices[k>1], slices[k>1]):
quad_grid[xmin:xmax, ymin:ymax] = i
In [195]:
plt.imshow(quad_grid, cmap='Set2', vmax=200)
quad_grid.dtype
Out[195]:
dtype('int32')
In [196]:
def C(x):
N = matplotlib.colors.Normalize(vmin=-10, vmax=300)
colors = matplotlib.cm.Blues(N(x), alpha=1)
# add a row completely transparent
colors = np.r_[colors, np.array([0,0,0,0])[np.newaxis,:]]
return colors
In [197]:
s1 = subgrid.get_nd('s1')[1:]
colors = C(s1)
plt.imshow(colors[quad_grid.filled(colors.shape[0]-1)], cmap='Blues')
Out[197]:
<matplotlib.image.AxesImage at 0x18b05450>
In [ ]:
Content source: nens/python-subgrid
Similar notebooks: