In [2]:

    

import python_subgrid.wrapper
import python_subgrid.plotting
import matplotlib.cm
import matplotlib.colors
import matplotlib.pyplot as plt
import numpy as np


    

    




/home/fedor/.virtualenvs/main/local/lib/python2.7/site-packages/pkg_resources.py:1054: UserWarning: /home/fedor/.python-eggs is writable by group/others and vulnerable to attack when used with get_resource_filename. Consider a more secure location (set with .set_extraction_path or the PYTHON_EGG_CACHE environment variable).
  warnings.warn(msg, UserWarning)

In [3]:

    

subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/home/fedor/Checkouts/models/kaapstad_centrum/kaapstad_centrum.mdu")
# subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/media/fedor/AC203/Test gebiedje WF/test_wf_snelheid_quad.mdu")


    

In [4]:

    

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.1032M, Jan 26 2014, 13:29:46
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'.
INFO:python_subgrid.wrapper:NetworkFile is not specified. Channel network is disabled.
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:           3   3.94224596      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:        2033   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 =     16473
INFO:python_subgrid.wrapper: Boundary nodes =   0
DEBUG:python_subgrid.wrapper:        2035   1.64759997E-02  nodk,nodm,nodn
DEBUG:python_subgrid.wrapper:        2035   1.64730009E-02  ks(nodall)
DEBUG:python_subgrid.wrapper:        2035   3.29460017E-02  xz,yz
INFO:python_subgrid.wrapper: Internal 1D nodes =         0
INFO:python_subgrid.wrapper: Boundary 1D nodes =   0
DEBUG:python_subgrid.wrapper:        2035   6.66259974E-02  line
DEBUG:python_subgrid.wrapper:        2035  0.133251995      lik,lim,lin,kcu,dpumin,dpumax
DEBUG:python_subgrid.wrapper:        2035  0.131791994      nru,nrv
INFO:python_subgrid.wrapper:Initializing boundaries took:   0.2000D-02 seconds.
DEBUG:python_subgrid.wrapper:        2035   3.33129987E-02  flod
DEBUG:python_subgrid.wrapper:        2035   3.33129987E-02  flou
INFO:python_subgrid.wrapper:Grid Admin initialisation took:  0.2400D-01 seconds.
INFO:python_subgrid.wrapper:Completed initialization of grid administration.
DEBUG:python_subgrid.wrapper:        2035   1.64730009E-02  dmax(1:nodall)
DEBUG:python_subgrid.wrapper:        2035   1.64730009E-02  dmin(1:nodall)
DEBUG:python_subgrid.wrapper:        2035   1.64730009E-02  sumax(1:n2dtot)
DEBUG:python_subgrid.wrapper:        2035   6.58920035E-02  dqmin(1:nodtot)
DEBUG:python_subgrid.wrapper:        2036   6.58920035E-02  dqmax(1:nodtot)
DEBUG:python_subgrid.wrapper:        2036   6.58920035E-02  suqmax(1:nodtot)
DEBUG:python_subgrid.wrapper:        2036   6.58920035E-02  ksq(1:nodtot)
DEBUG:python_subgrid.wrapper:        2036   4.94219996E-02  s1(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   4.94219996E-02  si(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   9.99419987E-02  u1(0:linall)
DEBUG:python_subgrid.wrapper:        2036   3.33140008E-02  kf(0:linall)
DEBUG:python_subgrid.wrapper:        2036   1.64739992E-02  vol1(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   3.29479985E-02  vol2(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   3.29479985E-02  su(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   3.29479985E-02  volErr(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   3.33129987E-02  conv(linall)
DEBUG:python_subgrid.wrapper:        2036   6.58959970E-02  vq1(0:4*nodall)
DEBUG:python_subgrid.wrapper:        2036   1.64739992E-02  implis(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   8.23699981E-02  aii,bi,d0,a00,b0(0:nodall)
DEBUG:python_subgrid.wrapper:        2036   6.58959970E-02  vol0(0:nodall),vq0(0:4*nodall)
DEBUG:python_subgrid.wrapper:        2037  0.263572007      vqt,denom(0:4*nodall),fbq(0:4*nodall),umq(0:4*nodall)
DEBUG:python_subgrid.wrapper:        2037   9.88409966E-02  qs,qs1,qs2(0:2*nodall)
DEBUG:python_subgrid.wrapper:        2037   9.99419987E-02  u0,uh(0:linall)
DEBUG:python_subgrid.wrapper:        2037   1.64739992E-02  uc(0:nodall,2)
DEBUG:python_subgrid.wrapper:        2037   1.64739992E-02  vnorm(0:nodall)
DEBUG:python_subgrid.wrapper:        2037  0.166565001      cu,ru(linall)
DEBUG:python_subgrid.wrapper:        2037   9.99419987E-02  au(linall,0:2)
DEBUG:python_subgrid.wrapper:        2037   3.33129987E-02  cz(linall)
DEBUG:python_subgrid.wrapper:        2037   6.66280016E-02  qh1,qh2(0:linall)
DEBUG:python_subgrid.wrapper:        2037   6.66280016E-02  ade,adi(0:linall)
INFO:python_subgrid.wrapper:Calculation of bathymetry statistics took:   0.6160D+00 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 [5]:

    

subgrid.update(-1)


    

    




INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init=              3
DEBUG:python_subgrid.wrapper:        2037   6.66259974E-02  tables
DEBUG:python_subgrid.wrapper:        2037   6.58920035E-02  tables
DEBUG:python_subgrid.wrapper:        2037   6.58920035E-02  tables
DEBUG:python_subgrid.wrapper:        2037   6.58920035E-02  tables
DEBUG:python_subgrid.wrapper:        2038   6.66259974E-02  tables
DEBUG:python_subgrid.wrapper:        2038   6.58920035E-02  tables
DEBUG:python_subgrid.wrapper:        2038   6.58920035E-02  tables
DEBUG:python_subgrid.wrapper:        2038   1.64730009E-02  tables
DEBUG:python_subgrid.wrapper:        2038   1.64730009E-02  tables
DEBUG:python_subgrid.wrapper:        2043   4.79712820      autb
DEBUG:python_subgrid.wrapper:        2047   4.66408014      vltb
DEBUG:python_subgrid.wrapper:        2060   13.0855522      vq1tb,dentb
INFO:python_subgrid.wrapper:Initialization of Flow Values took:   0.5530D+00 seconds.
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  ij(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  noel(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  noddgr(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  nodstk(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  row(nodtot)
DEBUG:python_subgrid.wrapper:        2060   3.33129987E-02  ijl(lintot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  ia(nodtot)
INFO:python_subgrid.wrapper:reduce time:       0.071000
INFO:python_subgrid.wrapper: noactive=     16473 nogauss=      8390 nocg      8083
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  pk(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  rk(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  zk(nodtot)
DEBUG:python_subgrid.wrapper:        2060   1.64730009E-02  apk(nodtot)
DEBUG:python_subgrid.wrapper:        2061   6.83180019E-02  aij(ijtot)
DEBUG:python_subgrid.wrapper:        2061   6.83180019E-02  aij0(ijtot)
DEBUG:python_subgrid.wrapper:        2061   3.33129987E-02  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.0000E+00 idthalf 1  np = 1
INFO:python_subgrid.wrapper:Timestep complete.






    
Out[5]:





0

In [6]:

    

quad_grid = python_subgrid.plotting.make_quad_grid(subgrid)

vars = {var: subgrid.get_nd(var ,sliced=True) for var in {'s1', 'vol1'}}
plt.imshow(vars['s1'][quad_grid][::4,::4], cmap ='Blues')


    

    
Out[6]:





<matplotlib.image.AxesImage at 0xd642190>






    

In [7]:

    

grid = {}
for var in {'nodn', 'nodm', 'nodk', 'dxp', 'dx', 'nmax', 'mmax', 'imaxk', 'jmaxk', 'imax', 'jmax', 'dps', 'x0p', 'y0p'}:
    grid[var] = subgrid.get_nd(var, sliced=True)

m = (grid['nodm']-1)*grid['imaxk'][grid['nodk']-1]
n = (grid['nodn']-1)*grid['jmaxk'][grid['nodk']-1]
size = grid['imaxk'][grid['nodk']-1]


quad_grid_shape = (grid['jmax'], grid['imax'])
quad_grid = np.ma.empty(quad_grid_shape, dtype='int32')
quad_grid.mask = True
quad_grid.fill_value = -1
for i, (m_i, n_i, size_i) in enumerate(zip(m,n,size)):
    quad_grid[n_i:(n_i+ size_i), m_i:(m_i+size_i)] = i


fig, ax = plt.subplots(figsize=(13,8))

ax.axis('equal')
ax.imshow(quad_grid, vmax=30)


    

    
Out[7]:





<matplotlib.image.AxesImage at 0xc6eff50>






    

In [8]:

    

quad_x = grid['mmax'][0]*grid['dx'][0]
quad_y = grid['nmax'][0]*grid['dx'][0]
perc_x = quad_grid.shape[1]/(grid['mmax'][0]*grid['dxp'])
perc_y = quad_grid.shape[0]/(grid['nmax'][0]*grid['dxp'])


    

In [9]:

    

fig, ax = plt.subplots(figsize=(10,6))
dps_extent=(grid['x0p'], grid['x0p']+grid['dxp']*grid['imax'], grid['y0p']+grid['dxp']*grid['jmax'], grid['y0p'] )
print(dps_extent)

quad_extent= (grid['x0p'], grid['x0p']+grid['dxp']*quad_grid.shape[1], grid['y0p']+grid['dxp']*quad_grid.shape[0], grid['y0p'] )
print(quad_extent)
ax.imshow(np.ma.masked_less(grid['dps'], -10000), extent=dps_extent)
ax.imshow(quad_grid, cmap='Set2', alpha=0.5,  extent=quad_extent)
ax.set_ylim(quad_extent[3], quad_extent[2])
#ax.plot(mc, nc, 'k+', alpha=0.5)


    

    




(array(255580.19380886), 267760.19380886003, 6247234.4410560122, array(6239159.441056012))
(array(255580.19380886), 267760.19380886003, 6247234.4410560122, array(6239159.441056012))






    
Out[9]:





(array(6239159.441056012), 6247234.4410560122)






    

In [10]:

    

import scipy.interpolate
mc = m + size/2.0
nc = n + size/2.0

N, M = np.mgrid[0:grid['dps'].shape[0], 0:grid['dps'].shape[1]]
fig, axes = plt.subplots(1, 2, figsize=(20,8), sharey=True, sharex=True)
axes[0].set_title('Non interpolated')
axes[0].imshow(python_subgrid.plotting.colors(vars['s1'], cmap='jet')[quad_grid.filled(-1)])
    L = scipy.interpolate.LinearNDInterpolator(np.c_[mc, nc], vars['s1'])
axes[0].triplot(L.points[:,0], L.points[:,1], L.tri.vertices, alpha=0.5)
axes[1].set_title('Linear cell center')
axes[1].imshow(L(M, N), alpha=0.3)
axes[1].triplot(L.points[:,0], L.points[:,1], L.tri.vertices, alpha=0.5)


    

In [11]:

    

L = scipy.interpolate.LinearNDInterpolator(np.c_[mc, nc], vars['s1'])


    

In [12]:

    

dps = grid['dps']
waterlevel = L(M, N)


    

In [13]:

    

dps.shape, waterlevel.shape, quad_grid.shape


    

    
Out[13]:





((1615, 2436), (1615, 2436), (1615, 2436))

In [14]:

    

dps_ma = np.ma.masked_less(dps, -10000)
fig, ax= plt.subplots(figsize=(20,10))
ax.imshow(np.ma.masked_less(waterlevel - -dps_ma, 0), cmap='Blues')


    

    
Out[14]:





<matplotlib.image.AxesImage at 0x164e5410>






    

In [15]:

    

fig, axes = plt.subplots(1,5, figsize=(20,6))
im = axes[0].imshow(np.ma.masked_less(waterlevel - -dps_ma, 0))
plt.colorbar(im, ax=axes[0])
im = axes[1].imshow(waterlevel)
plt.colorbar(im, ax=axes[1])
im = axes[2].imshow(-dps_ma)
plt.colorbar(im, ax=axes[2])
slice = np.s_[:dps_ma.shape[0]-10, :dps_ma.shape[1]]
im = axes[3].imshow(np.ma.masked_less(vars['s1'][quad_grid[slice].filled(-1)] - -dps_ma[slice],0))
plt.colorbar(im, ax=axes[3])
im = axes[4].imshow(quad_grid[slice], vmax=20)


    

In [16]:

    

fig, axes = plt.subplots(1,3,figsize=(15,7), sharex=True, sharey=True)
slice = np.s_[:1500, :2000]
slice = np.s_[300:500,300:500]
slice = np.s_[1200:1400, 1000:1200]
axes[0].imshow(np.ma.masked_less(waterlevel[slice] - -dps_ma[slice], 0))
im = axes[1].imshow(waterlevel[slice]- -dps_ma[slice], vmin=0)
plt.colorbar(im, ax=axes[1])

im = axes[2].imshow(np.ma.masked_less(vars['s1'][quad_grid[slice].filled(-1)] - -dps_ma[slice],0))


    

In [1]:

    

import scipy.interpolate


    

In [3]:

    

x = [0,0,1,1]
y = [0,1,0,1]
v = [1,1,9,1]
L = scipy.interpolate.LinearNDInterpolator(np.c_[x,y], v)
print(L(0.5,  0.5))
plt.imshow(L(*np.mgrid[0:1:10j,0:1:10j]), interpolation='none')
plt.colorbar()
triplot(L.points[:,0], L.points[:,1], L.tri.vertices)


    

    




1.0






    

In [20]:

    

L = scipy.interpolate.LinearNDInterpolator(np.c_[mc, nc], vars['s1'])
wl = np.ma.masked_less(L(M, N) - -dps_ma, 0)

L = scipy.interpolate.LinearNDInterpolator(np.c_[mc, nc], vars['vol1'])
wl_ma = np.ma.masked_array(wl, mask=np.logical_or(L(M, N) <= 0.1, np.isnan(L(M, N))))
plt.imshow(wl_ma, cmap='Blues')

plt.colorbar()


    

    
Out[20]:





<matplotlib.colorbar.Colorbar instance at 0xc703098>






    

In [27]:

    

iswet = np.ma.masked_array(vars['vol1'][quad_grid.filled(0)], quad_grid.mask) > 0


    

In [28]:

    

plt.imshow(iswet)
plt.colorbar()


    

    
Out[28]:





<matplotlib.colorbar.Colorbar instance at 0x1720cb48>






    

In [129]:

    

%%timeit

N = matplotlib.colors.Normalize(vars['s1'].min(), vars['s1'].max())
C = matplotlib.cm.Blues
colors = colors = C(N(vars['s1']))
colors = np.r_[colors, np.array([0,0,0,0])[np.newaxis,:]]

kmin = grid['jmaxk'][0]
img = colors[quad_grid[::kmin, ::kmin].filled()]


    

    




10 loops, best of 3: 22.3 ms per loop

In [202]:

    

import pandas
mask = np.logical_or(quad_grid.mask, dps_ma.mask)
idx = quad_grid[~mask].ravel()
dps_flat = dps_ma[~mask].ravel()
df = pandas.DataFrame(data=dict(quad=idx, topo=-dps_flat))
quad_min = df.groupby('quad').min()
quad_min.sort_index(inplace=True)
# validity check
assert quad_min.index.tolist() == range(quad_grid.min(), quad_grid.max()+1)
topo_min = np.ma.masked_array(np.array(quad_min['topo'])[quad_grid.filled()])
topo_min.mask = quad_grid.mask


    

In [264]:

    

topo_diff = (-dps_ma - topo_min)
np.percentile(topo_diff[~topo_diff.mask].flatten(), 80)


    

    
Out[264]:





165.85341186523439

In [266]:

    

plt.imshow(waterlevel - topo_min)


    

    
Out[266]:





<matplotlib.image.AxesImage at 0x13bf80f90>






    

In [ ]: