In [47]:

    

import python_subgrid.wrapper
import logging


    

    




/home/fedor/.virtualenvs/main/local/lib/python2.7/site-packages/pkg_resources.py:991: 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 [48]:

    

subgrid = python_subgrid.wrapper.SubgridWrapper('/home/fedor/Checkouts/models/brouwersdam/brouwersdam.mdu')
python_subgrid.wrapper.logger.setLevel(logging.WARN)


    

In [49]:

    

subgrid.start()


    

    




WARNING:python_subgrid.wrapper:reading *.ext forcings file  getting unknown QUANTITY initialwaterlevel
WARNING:python_subgrid.wrapper:reading *.ext forcings file  getting unknown QUANTITY initialwaterlevel

In [50]:

    

import python_subgrid.plotting


    

In [51]:

    

rows = []
for i in range(500):
    s1 = subgrid.get_nd('s1')[320].copy()
    t1 = subgrid.get_nd('t1').copy()
    rows.append(dict(s1=s1, t1=t1, i=i))
    subgrid.update(-1)


    

In [56]:

    

import pandas
df = pandas.DataFrame.from_records(rows)
plt.plot(df['t1']/3600.0, df['s1'])


    

    
Out[56]:





[<matplotlib.lines.Line2D at 0xa9b0290>]






    

In [5]:

    

names = {'FlowElem_xcc', 'FlowElem_ycc', 
         'FlowLink_xu', 'FlowLink_yu', 's1', 
         'nod_type', 'link_type',
         'FlowElemContour_x', 'FlowElemContour_y', 'x0p', 'x1p', 'imax', 
         'jmax', 'y0p', 'y1p', 'dxp', 'dps', 'dx'}
vars = {var: subgrid.get_nd(var) for var in names}


    

In [6]:

    

print vars['x1p'], vars['x0p'] + vars['imax']*vars['dxp']


    

    




87955.0 87955.0

In [7]:

    

# compute an index for all 2d cells
idx2d = vars['nod_type'] == 1
idx1d = vars['nod_type'] == 2
python_subgrid.wrapper.NODE_TYPES, idx2d.sum(), idx1d.sum()


    

    
Out[7]:





({0: 'internal boundary',
  1: '2d',
  2: '1d',
  3: '2d boundary',
  4: '1d boundary'},
 6608,
 0)

In [8]:

    

# compute the coordinates for the bathymetry
slice = np.s_[::4] # don't plot everything
x = np.arange(vars['x0p'], vars['x1p'], vars['dxp'])
y = np.arange(vars['y0p'], vars['y1p'], vars['dxp'])


    

In [34]:

    

orifices = subgrid.get_nd('orifices')

orifice = orifices.set_index('id').ix['3']


    

In [35]:

    

# plot everything
fig, ax = plt.subplots(figsize=(15,8))
X,Y=np.meshgrid(x[slice], y[slice])
print(X.shape, Y.shape)
#ax.pcolormesh(X,Y, -vars['dps'][:,:-1][slice, slice], vmin=-20, vmax=20, cmap='gist_earth', alpha=0.5)
ax.pcolorfast(x[slice],y[slice], -vars['dps'][:,:-1][slice, slice], vmin=-20, vmax=20, cmap='gist_earth', alpha=0.5)
ax.plot(vars['FlowElem_xcc'][idx2d], vars['FlowElem_ycc'][idx2d], 'k.', alpha=0.3)
#ax.plot(vars['FlowElem_xcc'][idx1d], vars['FlowElem_ycc'][idx1d], 'ro', alpha=0.3)

ax.plot(vars['FlowElemContour_x'].ravel(), vars['FlowElemContour_y'].ravel(), 'k+', alpha=0.3)
ax.axis('equal')
ax.set_xlim((45000,51000))
ax.set_ylim(417000, 424000)
for i, orifice in orifices.iterrows():
    x0 = vars['FlowElem_xcc'][orifice['left_calc_point']]
    y0 = vars['FlowElem_ycc'][orifice['left_calc_point']]
    x1 = vars['FlowElem_xcc'][orifice['right_calc_point']]
    y1 = vars['FlowElem_ycc'][orifice['right_calc_point']]
    lx = vars['FlowLink_xu'][orifice['link_number']]
    ly = vars['FlowLink_yu'][orifice['link_number']]
    ax.plot([x0, x1], [y0, y1], 'r-', linewidth=3 )
    #ax.annotate(orifice['id'],(lx, ly))
    ax.text(lx,  ly, orifice['id'].strip())


    

    




((2377, 2701), (2377, 2701))






    

In [36]:

    

orifices[['left_calc_point', 'right_calc_point', 'link_number']]


    

    
Out[36]:






  

    

      

      
left_calc_point
      
right_calc_point
      
link_number
    
  
  

    

      
0
      
 5258
      
 5284
      
 12927
    
    

      
1
      
 5230
      
 5256
      
 12934
    
    

      
2
      
 5964
      
 6000
      
 12940
    
    

      
3
      
 5927
      
 5961
      
 12945
    
    

      
4
      
 6397
      
 6413
      
 12951
    
    

      
5
      
 6382
      
 6396
      
 12955
    
    

      
6
      
 5263
      
 5289
      
 12961
    
    

      
7
      
 5156
      
 6527
      
 12966
    
  

In [37]:

    

#this should work 
subgrid.set_structure_field('orifices', orifice['id'], 'crest_level', -1)
subgrid.get_nd('orifices').set_index('id').ix[orifice['id']]['crest_level']


    

    
Out[37]:





-1.0

In [38]:

    

orifices['crest_level']


    

    
Out[38]:





0   -10
1   -10
2   -30
3   -10
4   -10
5   -10
6   -10
7   -10
Name: crest_level, dtype: float64

In [91]:

    

V = {}
for i in range(0, 300):
    subgrid.update(-1)
    V[i] = {var: subgrid.get_nd(var).copy() for var in {'q', 's1'}}


    

In [40]:

    

print(orifice)
# Lookup waterlevel and discharge for this structure
V[4]['s1'][orifice['left_calc_point']], V[4]['s1'][orifice['right_calc_point']], V[4]['q'][orifice['link_number']]


    

    




allowed_flow_dir                 1
branchid                         8
chainage                  369.9438
contraction_coeff                1
crest_level                    -10
crest_width                     20
id                               8
lat_contr_coeff                0.9
left_calc_point               5156
limit_flow_neg       6.953179e-310
limit_flow_pos                   0
link_number                  12966
open_level                    -5.2
right_calc_point              6527
x                                0
y                                0
Name: 7, dtype: object






    
Out[40]:





(-1.0608190854262996e-13, 2.3079331074216172, 0.0)

In [41]:

    

quad_grid = python_subgrid.plotting.make_quad_grid(subgrid)


    

In [42]:

    

s1 = V[4]['s1'].copy()
s1[0] = -3

plt.imshow(s1[quad_grid[::4,::4].filled(0)], origin='lower', vmin=-2, vmax=2, cmap='gist_earth')
plt.colorbar()


    

    
Out[42]:





<matplotlib.colorbar.Colorbar instance at 0xcee3950>






    

In [77]:

    

fig, axes = plt.subplots(1,3, figsize=(13,4))
nod_type = subgrid.get_nd('nod_type')
link_type = subgrid.get_nd('link_type')
axes[0].plot([v['s1'][nod_type == 1][orifice['left_calc_point']-1] for v in V.values()])
axes[1].plot([v['s1'][nod_type == 1][orifice['right_calc_point']-1] for v in V.values()])
axes[2].plot([v['q'][link_type >= 1][orifice['link_number']-1] for v in V.values()])


    

    
Out[77]:





[<matplotlib.lines.Line2D at 0x6b80910>]






    

In [45]:

    

import netCDF4
ds = netCDF4.Dataset('/home/fedor/Checkouts/models/brouwersdam/subgrid_map.nc')
x = ds.variables['FlowElem_xcc'][:]
y = ds.variables['FlowElem_ycc'][:]
plt.plot(x,y, 'k.')

plt.xlim(45000, 50000)
plt.ylim(400000, 430000)
for i, (x,y) in enumerate(zip(x,y)):
    if 45000 < x < 50000 and 400000 < y < 430000:
        plt.text(x,y, str(i))


    

In [46]:

    

t = netCDF4.num2date(ds.variables['time'][:], ds.variables['time'].units)
plt.plot(t, ds.variables['s1'][:,319])


    

    
Out[46]:





[<matplotlib.lines.Line2D at 0xdb9bed0>]






    

In [ ]: