notebook.community

Edit and run



In [1]:

    
import osgeo.gdal
import osgeo.osr

try:
    from PIL import Image
except ImportError:
    import Image

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
logger = logging.getLogger(__name__)









    



/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 [2]:

    
subgrid = python_subgrid.wrapper.SubgridWrapper(mdu="/home/fedor/Checkouts/models/kaapstad_centrum/kaapstad_centrum.mdu")
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.1079:1081M, Feb  3 2014, 10:04:16
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 '...
INFO:python_subgrid.wrapper: Projection information: PROJCS["Cape / UTM zone 34S",GEOGCS["Cape",DATUM["Cape",SPHEROID["Clarke 1880 (Arc)",6378249.145,293.4663076999908,AUTHORITY["EPSG","7013"]],TOWGS84[-136,-108,-292,0,0,0,0],AUTHORITY["EPSG","6222"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4222"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",21],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","22234"]]
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.5000D-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.2900D-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.5440D+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 [25]:

    
for i in range(5):
    subgrid.update(-1)









    



INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init=              0
INFO:python_subgrid.wrapper:finding rain data for input
INFO:python_subgrid.wrapper:Timestep 11 number of iterations 1 epsmax   0.0000E+00 idthalf 1  np = 1
INFO:python_subgrid.wrapper:Timestep complete.
INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init=              0
INFO:python_subgrid.wrapper:finding rain data for input
INFO:python_subgrid.wrapper:Timestep 12 number of iterations 1 epsmax   0.0000E+00 idthalf 1  np = 1
INFO:python_subgrid.wrapper:Timestep complete.
INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init=              0
INFO:python_subgrid.wrapper:finding rain data for input
INFO:python_subgrid.wrapper:Timestep 13 number of iterations 1 epsmax   0.0000E+00 idthalf 1  np = 1
INFO:python_subgrid.wrapper:Timestep complete.
INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init=              0
INFO:python_subgrid.wrapper:finding rain data for input
INFO:python_subgrid.wrapper:Timestep 14 number of iterations 1 epsmax   0.0000E+00 idthalf 1  np = 1
INFO:python_subgrid.wrapper:Timestep complete.
INFO:python_subgrid.wrapper:Start timestep-update.
DEBUG:python_subgrid.wrapper:Start initFlowAdmin, flow_init=              0
INFO:python_subgrid.wrapper:finding rain data for input
INFO:python_subgrid.wrapper:Timestep 15 number of iterations 1 epsmax   0.0000E+00 idthalf 1  np = 1
INFO:python_subgrid.wrapper:Timestep complete.



In [4]:

    
drivers = {osgeo.gdal.GetDriver(i).ShortName:osgeo.gdal.GetDriver(i) for i in range(osgeo.gdal.GetDriverCount())}



In [5]:

    
driver = drivers['GTiff']



In [6]:

    
import python_subgrid.plotting
quad_grid = python_subgrid.plotting.make_quad_grid(subgrid)
grid = {}
for var in {'x0p', 'x1p', 'y0p', 'y1p', 'dps', 'dxp', 'dyp', 'imax', 'jmax', 'imaxk', 'nodk', 'nodm', 'nodn'}:
    grid[var] = subgrid.get_nd(var)
    
quad_transform= (grid['x0p'],  # xmin
                 grid['dxp'], # xmax
                 0,            # for rotation
                 grid['y0p'],
                 0,
                 grid['dyp'])
print quad_transform









    



(array(255580.19380886), array(5.0), 0, array(6239159.441056012), 0, array(5.0))



In [32]:

    
vars = {var: subgrid.get_nd(var)[1:] for var in {'s1', 'vol1', 'su'}}
vars['dmax'] = subgrid.get_nd('dmax')
vars['nodk'] = grid['nodk']
import pandas
df = pandas.DataFrame(vars)
df[np.logical_and.reduce([df.vol1>0])].head()









    Out[32]:






  
    
      
      dmax
      nodk
      s1
      su
      vol1
    
  
  
    
      1 
       37.668747
       8
      -10
       200
       2475.376415
    
    
      3 
       41.190121
       8
      -10
        25
        779.753017
    
    
      25
       13.222917
       6
      -10
        50
        153.018069
    
    
      38
       92.939850
       6
      -10
        50
       4089.230919
    
    
      40
       16.450367
       6
      -10
        25
        161.259174



In [7]:

    
grid['y1p'], grid['y0p']









    Out[7]:





(array(6247234.441056012), array(6239159.441056012))



In [8]:

    
src_ds = driver.Create("input.tiff", quad_grid.shape[1], quad_grid.shape[0], 4, eType=osgeo.gdal.GDT_Byte)



In [9]:

    
src_ds.SetGeoTransform(quad_transform)
src_srs = osgeo.osr.SpatialReference()
epsg = 22234
src_srs.ImportFromEPSG(epsg)
src_ds.SetProjection(src_srs.ExportToWkt())









    Out[9]:





0



In [31]:

    
mask = quad_grid.mask
s1 = vars['s1'][quad_grid.filled(0)]
vol1 = vars['vol1'][quad_grid.filled(0)]
vol1 = np.ma.masked_array(vol1, quad_grid.mask)
dps_ma = np.ma.masked_less(grid['dps'][1:-1,1:-1], -9000)
waterlevel = s1 - (-dps_ma)
mask = np.logical_or.reduce([quad_grid.mask, dps_ma.mask, vol1 == 0])

area = grid['imaxk'][grid['nodk']-1]*grid['dxp']*grid['dyp']
volumelevel = vars['vol1']/area

fig, axes = plt.subplots(1,5, figsize=(30,6))
im1 = axes[0].imshow(-dps_ma, vmin=-1000, vmax=1000, cmap='gist_earth')
axes[0].set_title('-dps')
plt.colorbar(im1, ax=axes[0])
im = axes[1].imshow(np.ma.masked_array(s1, mask=quad_grid.mask), vmin=-1000, vmax=1000, cmap='gist_earth')
plt.colorbar(im, ax=axes[1])
axes[1].set_title('s1')
im = axes[2].imshow(s1 - (-dps_ma), vmin=0, cmap='Blues')
plt.colorbar(im, ax=axes[2])
axes[2].set_title('s1--dps')
im = axes[3].imshow(vol1, cmap='Blues', vmin=0)
plt.colorbar(im, ax=axes[3])
axes[3].set_title('vol1')


N = matplotlib.colors.Normalize(0, waterlevel.max())
img = matplotlib.cm.Blues(N(waterlevel), bytes=True)
axes[-1].imshow(img.astype('uint8'))
axes[-1].set_title('waterlevel')
scalar = matplotlib.cm.ScalarMappable(norm = N, cmap=matplotlib.cm.Blues)
scalar.set_array(waterlevel)
plt.colorbar(scalar, ax=axes[-1])
img_rolled = np.rollaxis(img, 2, 0)



In [33]:

    
import pandas 
df = pandas.DataFrame(data=dict(idx=quad_grid.ravel(), dps=np.ma.masked_less(grid['dps'][1:-1,1:-1],-9000).ravel()))
maxs = df.groupby('idx').max().reset_index()
dps_max = np.take(np.array(maxs['dps']), quad_grid)
cell_topo = - np.ma.masked_less(grid['dps'][1:-1,1:-1], -9000)  + dps_max
#plt.imshow( )


#plt.colorbar()
mask = (volumelevel[quad_grid.filled(0)] - cell_topo) < 0
level = volumelevel[quad_grid.filled(0)] - cell_topo
level.mask = np.logical_or(level.mask, mask)
plt.imshow(level, cmap='Blues')
plt.colorbar()









    Out[33]:





<matplotlib.colorbar.Colorbar instance at 0x11598170>



In [ ]:

    
for i in range(img_rolled.shape[0]):
    band = src_ds.GetRasterBand(i+1)
    band.WriteArray(img_rolled[i])
src_ds.FlushCache()

img2_rolled  = src_ds.ReadAsArray()
img2 = np.rollaxis(img2_rolled, 0, 3)
plt.imshow(img2)



In [ ]:

    
width= 256
height= 256
srs= "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over"
xmin= 2044843.38069
ymin= -4021199.18403
xmax= 2054627.32031
ymax= -4011415.24441
zoom= 12



In [ ]:

    
area_ds = driver.Create('output.tiff', width, height, 4, eType = osgeo.gdal.GDT_Byte)
if area_ds is None:
    raise Exception('uh oh.')

merc = osgeo.osr.SpatialReference()
merc.ImportFromProj4(srs)

area_ds.SetProjection(merc.ExportToWkt())

w = xmax - xmin
h = ymax - ymin
print(h)
gtx = [xmin, w/width, 0, ymin, 0, h/height]
area_ds.SetGeoTransform(gtx)


# Adjust resampling method -----------------------------------------

resample = osgeo.gdal.GRA_NearestNeighbour

# Create rendered area ---------------------------------------------


osgeo.gdal.ReprojectImage(src_ds, area_ds, src_ds.GetProjection(), area_ds.GetProjection(), resample, 100000, 0.1 )
area_ds.FlushCache()
data = area_ds.ReadAsArray()
print data.shape
try:
    #driver.Delete('output.tiff')
    pass
except:
    logger.exception("can't delete output.tiff")
try:
    #driver.Delete('input.tiff')
    pass
except:
    logger.exception("can't delete input.tiff")
import PIL
area = PIL.Image.fromarray(np.rollaxis(data, 0,3))
area.save('test.png')



In [ ]:

    
plt.imshow(np.rollaxis(data, 0, 3))



In [ ]:

    
from IPython.display import Image
Image('test.png')



In [ ]:

    
subgrid.get_nd('dxp')



In [ ]:

    
grid['imax']



In [ ]:

	dmax	nodk	s1	su	vol1
1	37.668747	8	-10	200	2475.376415
3	41.190121	8	-10	25	779.753017
25	13.222917	6	-10	50	153.018069
38	92.939850	6	-10	50	4089.230919
40	16.450367	6	-10	25	161.259174