notebook.community

Edit and run

Import stuff





In [1]:



    


import os.path
import logging

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


import python_subgrid.wrapper

Input model and reference point





In [2]:



    


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

Start the model





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.
WARNING:python_subgrid.wrapper:File

Read the pixel grid





In [4]:



    


pixels = {}
pixels['dps'] = subgrid.get_nd('dps')[1:-1,1:-1].T # contains ghost cells? # why transposed
pixels['dps'] = np.ma.masked_array(-pixels['dps'], mask=pixels['dps']==pixels['dps'].min())
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'])

Plot the pixel grid





In [5]:



    


# Don't bother rendering everything
thin = 10
# Big picture
fig, ax = plt.subplots(figsize=(30,20))
# Add the grid
ax.pcolorfast(pixels['x'][::thin], 
              pixels['y'][::thin], 
              pixels['dps'][::thin,::thin], 
              vmin=-20, vmax=20, cmap='gist_earth')
# show the molen location
ax.plot(molen['x'], molen['y'], 'bo', markersize=50, alpha=0.3)
ax.plot(molen['x'], molen['y'], 'wo', markersize=10, alpha=1)



    


















    
Out[5]:








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










    
























In [6]:



    


# Access the data from the regenradar
url = 'http://opendap.nationaleregenradar.nl/thredds/dodsC/radar/TF0005_A/2013/10/01/RAD_TF0005_A_20131001000000.h5'
ds = netCDF4.Dataset(url)



    











In [7]:



    


!ncdump -h $url



    


















    






netcdf RAD_TF0005_A_20131001000000 {
dimensions:
	east = 500 ;
	north = 490 ;
	time = 8928 ;
variables:
	byte available(time) ;
		available:_lastModified = "1970-01-25T08:16:08Z" ;
		available:_FillValue = 0b ;
		available:_Unsigned = "true" ;
		available:_ChunkSize = 8928 ;
	double east(east) ;
		east:_lastModified = "1970-01-25T08:16:08Z" ;
		east:_ChunkSize = 500 ;
	double north(north) ;
		north:_lastModified = "1970-01-25T08:16:08Z" ;
		north:_ChunkSize = 490 ;
	int time(time) ;
		time:long_name = "time" ;
		time:calendar = "gregorian" ;
		time:units = "seconds since 2013-10-01" ;
		time:standard_name = "time" ;
		time:_lastModified = "1970-01-06T15:16:55Z" ;
		time:_Unsigned = "true" ;
		time:_ChunkSize = 2232 ;
	float precipitation(north, east, time) ;
		precipitation:_lastModified = "1970-01-25T08:16:08Z" ;
		precipitation:_FillValue = -9999.f ;
		precipitation:_ChunkSize = 20, 20, 24 ;
}

Get data for location of interest





In [8]:



    


rain = {}
rain['time'] = netCDF4.num2date(ds.variables['time'][:], ds.variables['time'].units)
rain['x'] = ds.variables['east'][:]
rain['y'] = ds.variables['north'][:]
rain['P'] = ds.variables['precipitation']
molen['idx'] = (np.abs((rain['y'] - molen['y'])).argmin(), np.abs((rain['x'] - molen['x'])).argmin())
print molen['idx']
print rain['P'].dimensions, rain['P'].shape, rain['y'].shape
molen['p'] = rain['P'][molen['idx'][0], molen['idx'][1],:]



    


















    






(249, 197)
(u'north', u'east', u'time') (490, 500, 8928) (490,)

Show the timeseries





In [9]:



    


# use this period
rain['t0'] = datetime.datetime(2013,10,12,12, 0)
rain['t1'] = datetime.datetime(2013,10,14, 0, 0)
# create a timeseries object
series = pandas.Series(data=molen['p'].filled(0), index=rain['time'])

# Create a few plots
fig, axes = plt.subplots(1,2, figsize=(20,4))
# show the whole series
ax = axes[0]
ax.plot(series.index, series)
pandas.rolling_sum(series, window=24*60/5.0, center=True).plot(title="rolling day sum [mm/day]", ax=ax)
ax.fill_betweenx(ax.get_ylim(), [rain['t0']]*2, [rain['t1']]*2, alpha=0.3)

# show the zoomed in series
ax = axes[1]
ax.plot(series.index, series*60/5.0)
pandas.rolling_sum(series, window=24*60/5.0, center=True).plot(title="rolling day sum [mm/day]", ax=ax)
ax.set_xlim(rain['t0'], rain['t1'])
rain['time_idx'] = np.where(np.logical_and(rain['time'] >= rain['t0'], rain['time'] < rain['t1']))[0]
rain['t0_idx'] = rain['time_idx'].min()
rain['t1_idx'] = rain['time_idx'].max()
rain['thalf_idx'] = int(np.floor((rain['t0_idx'] + rain['t1_idx'])/2))
rain['p'] = rain['P'][:,:,rain['thalf_idx']] # x, y expected by some routines

Create some colormaps that correspond with rain





In [10]:



    


# looking up to the sky
rgb = matplotlib.colors.hex2color('#87CEEB')
segmentdata = {'red': [], 'green': [], 'blue': []}
for value, color in zip((0, 0.2, 0.9, 1), ('#87CEEB', '#FFFFFF', '#888888', '#666666')):
    rgb = matplotlib.colors.hex2color(color)
    segmentdata['red'].append((value, rgb[0], rgb[0]))
    segmentdata['green'].append((value, rgb[1], rgb[1]))
    segmentdata['blue'].append((value, rgb[2], rgb[2]))
Rain = matplotlib.colors.LinearSegmentedColormap('rain', segmentdata) 

# looking down
segmentdata = {'red': [], 'green': [], 'blue': [], 'alpha': []}
for value, color, alpha in zip((0, 0.2, 0.9, 1), ('#FFFFFF', '#FFFFFF', '#888888', '#666666'), (0,0.5,0.7,0.9)):
    rgb = matplotlib.colors.hex2color(color)
    segmentdata['red'].append((value, rgb[0], rgb[0]))
    segmentdata['green'].append((value, rgb[1], rgb[1]))
    segmentdata['blue'].append((value, rgb[2], rgb[2]))
    segmentdata['alpha'].append((value, alpha, alpha))
Rain_alpha = matplotlib.colors.LinearSegmentedColormap('rain', segmentdata)



    











In [11]:



    


fig, axes = plt.subplots(1,2, figsize=(15,5))
ax = axes[0]
thin = 10
# show the sky with the grid on top
# 0.2mm /5 min == 0.2*12mm/h == 2.4 mm/h == 2.4*24mm/day == 57.6mm/day 
im = ax.pcolorfast(rain['x'], rain['y'], rain['p'].filled(0), vmin=0, vmax=0.2, cmap=Rain)
# plot orangy to show the location
ax.pcolorfast(pixels['x'][::thin], pixels['y'][::thin], np.zeros_like(pixels['dps'][::thin,::thin]), vmin=0, vmax=1, cmap='Oranges')
ax.set_xlabel('projected x coordinate [m] epsg:28992')
ax.set_ylabel('projected y coordinate [m] epsg:28992')
cb = plt.colorbar(im, ax=ax)
cb.set_label('precipitation [mm/5min]')

# now zoom in on the location
ax = axes[1]
thin = 10
ax.pcolorfast(pixels['x'][::thin], pixels['y'][::thin], pixels['dps'][::thin,::thin], vmin=-20, vmax=20, cmap='gist_earth')
im = ax.pcolorfast(rain['x'], rain['y'], rain['p'].filled(0), vmin=0, vmax=0.2, cmap=Rain_alpha)
ax.set_xlabel('projected x coordinate [m] epsg:28992')
ax.set_ylabel('projected y coordinate [m] epsg:28992')
cb = plt.colorbar(im, ax=ax)
cb.set_label('precipitation [mm/5min]')
ax.plot(molen['x'], molen['y'], 'bo', markersize=50, alpha=0.3)
ax.plot(molen['x'], molen['y'], 'wo', markersize=10, alpha=1)
ax.set_xlim(pixels['x'].min(), pixels['x'].max())
ax.set_ylim(pixels['y'].min(), pixels['y'].max())



    


















    
Out[11]:








(435200.94391317002, 460345.94391317002)

Interpolate





In [12]:



    


import scipy.interpolate
print rain['x'].shape, rain['y'].shape, rain['p'].shape
interp = {}


F = scipy.interpolate.RectBivariateSpline(
                                          rain['x'].ravel(), 
                                          rain['y'][::-1], 
                                          # reverse y axis and swap x,y
                                          np.swapaxes(rain['p'][::-1,:].filled(0), 0,1)
                                          )
# Create a new grid
interp['x'] = np.linspace(pixels['x'].min(), pixels['x'].max(), num=500)
interp['y'] = np.linspace(pixels['y'].min(), pixels['y'].max(), num=300)
interp['X'], interp['Y'] = np.meshgrid(interp['x'], interp['y'])
# Evaluate the interpolation function on the new grid
interp['z'] = F.ev(interp['X'].ravel(), interp['Y'].ravel())
fig, ax = plt.subplots()

interp['Z'] = interp['z'].reshape(interp['X'].shape)
im = ax.pcolorfast(interp['x'], interp['y'], interp['Z'], cmap=Rain, vmax=0.2)
plt.colorbar(im, ax=ax)
interp['Z'].shape



    


















    






(500,) (490,) (490, 500)










    
Out[12]:








(300, 500)










    
























In [13]:



    


# Now store the variable in the model
memcdf = netCDF4.Dataset("precipitation.nc", mode="w", diskless=True, persist=False)

memcdf.createDimension("nx", interp['x'].shape[0])
memcdf.createDimension("ny", interp['y'].shape[0])

var = memcdf.createVariable("x", datatype="double", dimensions=("nx",))
var[:] = interp['x']
var.standard_name = 'projected_x_coordinate'
var.units = 'm'

var = memcdf.createVariable("y", datatype="double", dimensions=("ny", ))
var[:] = interp['y']
var.standard_name = 'projected_y_coordinate'
var.units = 'm'

var = memcdf.createVariable("rainfall", datatype="double", dimensions=("ny", "nx"), fill_value=-9999)
var[:,:] = interp['Z']*(1/5.0)*(1/1000.0) #  mm/5min * 5min/min * m/mm -> m/min 
var.standard_name = 'precipitation'
var.coordinates = 'y x'
var.units = 'm/min'

s0 = subgrid.get_nd('s1').copy()



    











In [14]:



    


#memcdf.close()
subgrid.subscribe_dataset("precipitation.nc")



    


















    
Out[14]:








0

















In [15]:



    


print(subgrid.get_nd('t1'))
for i in range(5):
    subgrid.update(-1)
print(subgrid.get_nd('t1'))



    


















    






0.0
300.0

















In [16]:



    


s1 = subgrid.get_nd('s1')



    











In [17]:



    


np.allclose(s1, s0) # some water should have appeared



    


















    
Out[17]:








False

Show the distribution of qrain and lookup contours





In [18]:



    


qrain = subgrid.get_nd('qrain').copy()
xz = subgrid.get_nd('FlowElem_xcc').copy()
yz = subgrid.get_nd('FlowElem_ycc').copy()
xc = subgrid.get_nd('FlowElemContour_x').copy()
yc = subgrid.get_nd('FlowElemContour_y').copy()

bins = plt.hist(qrain)



    


















    
























In [19]:



    


xz.shape, qrain.shape, s1.shape, s0.shape



    


















    
Out[19]:








((44207,), (44206,), (44207,), (44207,))

















In [49]:



    


verts = [np.c_[x_i,y_i] for x_i,y_i in zip(xc.T, yc.T)]
polys = matplotlib.collections.PolyCollection(verts, alpha=0.6)
N = matplotlib.colors.Normalize(vmin=0,vmax=1e-6)

area = (xc[1,:]- xc[0,:]) * (yc[2,:]- yc[0,:] )
val = (s1-s0)[1:]
colors = matplotlib.cm.Blues(N(val))

Show the effect (s1-s0)





In [50]:



    


fig, axes = plt.subplots(1,2, figsize=(16,5))

ax = axes[0]
ax.pcolorfast(pixels['x'][::thin], pixels['y'][::thin], pixels['dps'][::thin,::thin], vmin=-20, vmax=20, cmap='gist_earth')
im = ax.pcolorfast(interp['x'], interp['y'], interp['Z'], vmin=0, vmax=0.2, cmap=Rain_alpha, alpha=0.5)


ax = axes[1]
polys.set_facecolors(colors)
polys.set_edgecolor('none')
polys.set_linewidth(0.1)
ax.pcolorfast(pixels['x'][::thin], pixels['y'][::thin], pixels['dps'][::thin,::thin], vmin=-20, vmax=20, cmap='gist_earth')
ax.add_collection(polys)
#im = ax.pcolorfast(interp['x'], interp['y'], interp['Z'], vmin=0, vmax=0.2, cmap=Rain_alpha, alpha=0.5)

ax.autoscale()
#ax.autoscale()



    


















    
























In [ ]:



    


# Don't close the dataset before the model has done the timestep
#subgrid.subscribe_dataset('precipitation.nc')
#subgrid.update(60.0)
# Clears the memory
memcdf.close()



    











In [ ]:

Content source: nens/python-subgrid

Similar notebooks:

notebook.community | gallery | about