In [1]:

    

import os.path
import logging

import osgeo.gdal
import netCDF4
import numpy as np
import matplotlib.pyplot as plt
import pandas

# for points in poly (faster than shapely)
import matplotlib.path


import python_subgrid.wrapper


    

In [2]:

    

mdu = '../../models/delfland-model-voor-3di/hhdlipad.mdu'
molen = {'x': 87519, 'y': 450328}


    

In [3]:

    

subgrid = python_subgrid.wrapper.SubgridWrapper(mdu=mdu)
python_subgrid.wrapper.logger.setLevel(logging.WARNING)
subgrid.start()


    

    




WARNING:python_subgrid.wrapper:Ignored 'Teta' in MDU file, please use 'IntegrationMethod' instead.
No handlers could be found for logger "progress"
ERROR:python_subgrid.wrapper:No network found, 1d boundaries are ignored.
WARNING:python_subgrid.wrapper:File 
ERROR:python_subgrid.wrapper:No 1d network found. Structures are disabled.

In [4]:

    

subgrid.get_var_shape('dps')


    

    
Out[4]:





(5032, 5322)

In [13]:

    

pixels = {}
pixels['dps'] = subgrid.get_nd('dps')[1:-1,1:-1] # contains ghost cells? # why transposed
pixels['dsnop'] = subgrid.get_nd('dsnop')
pixels['dps'] = np.ma.masked_array(-pixels['dps'], mask=pixels['dps']==pixels['dsnop'])
pixels['x0p'] = subgrid.get_nd('x0p')
pixels['dxp'] = subgrid.get_nd('dxp')
pixels['y0p'] = subgrid.get_nd('y0p')
pixels['dyp'] = subgrid.get_nd('dyp')
pixels['imax'] = subgrid.get_nd('imax')
pixels['jmax'] = subgrid.get_nd('jmax')
pixels['x'] = np.arange(pixels['x0p'], pixels['x0p'] + pixels['imax']*pixels['dxp'], pixels['dxp'])
pixels['y'] = np.arange(pixels['y0p'], pixels['y0p'] + pixels['jmax']*pixels['dyp'], pixels['dyp'])
pixels['Y'], pixels['X'] = np.meshgrid(pixels['y'], pixels['x'])


    

In [14]:

    

fig, ax = plt.subplots()
thin = 10
# Add the grid
ax.pcolorfast(pixels['x'][::thin], 
              pixels['y'][::thin], 
              pixels['dps'][::thin,::thin], 
              vmin=-20, vmax=20, cmap='gist_earth')


    

    
Out[14]:





<matplotlib.image.AxesImage at 0x65efb50>






    

In [15]:

    

# parallellogram
coords = [[77500,445000], [80000,447500], [80000,450000], [77500,447500]]


    

In [16]:

    

poly = matplotlib.path.Path(vertices =coords, closed=False)


    

In [17]:

    

ax.add_patch(matplotlib.patches.PathPatch(poly))
fig


    

    
Out[17]:

In [18]:

    

pts = np.c_[pixels['X'].ravel(), pixels['Y'].ravel()]


    

In [19]:

    

%%time
# keep this below a second
inpoly = poly.contains_points(pts)
pixels['inpoly'] = inpoly.reshape(pixels['X'].shape)


    

    




CPU times: user 561 ms, sys: 58.3 ms, total: 619 ms
Wall time: 619 ms

In [20]:

    

subgrid.get_nd('dps')[1:-1,1:-1][pixels['inpoly'].T] += 10


    

In [21]:

    

subgrid.update(-1)


    

    




WARNING:python_subgrid.wrapper:Progress stopped but no current progress message
WARNING:python_subgrid.wrapper:Progress stopped but no current progress message






    
Out[21]:





0

In [25]:

    

pixels['dps'] = subgrid.get_nd('dps')[1:-1,1:-1] # contains ghost cells? # why transposed
pixels['dsnop'] = subgrid.get_nd('dsnop')
pixels['dps'] = np.ma.masked_array(-pixels['dps'], mask=pixels['dps']==pixels['dsnop'])


    

In [26]:

    

fig, ax = plt.subplots()
thin = 10
# Add the grid
ax.pcolorfast(pixels['x'][::thin], 
              pixels['y'][::thin], 
              pixels['dps'][::thin,::thin], 
              vmin=-20, vmax=20, cmap='gist_earth')


    

    
Out[26]:





<matplotlib.image.AxesImage at 0x69dcd10>






    

In [24]:

    




    

In [ ]: