In [4]:
import os
import knmi
import scipy
import netCDF4
import waterbase
import pandas as pd
import cPickle as pickle
from datetime import datetime
root = '/Users/hoonhout/GitHub/aeolis-models/sandmotor/'
t0 = datetime(2011,8,1)
t1 = datetime(2012,8,1)
key = '20110803'
In [2]:
# read bathy
with open('/Users/hoonhout/Checkouts/Zandmotor/morphology/JETSKI/bathy_gridded.pkl', 'rb') as fp:
grid = pickle.load(fp)
In [5]:
def rotate(x, y, phi, origin=(72502.0, 452074.0)):
x = np.asarray(x)
y = np.asarray(y)
phi = phi / 180. * np.pi
R = np.asarray([[np.cos(phi), -np.sin(phi)],
[np.sin(phi), np.cos(phi)]])
XY = np.concatenate((x.reshape(-1,1) - origin[0],
y.reshape(-1,1) - origin[1]), axis=1)
XY = np.dot(XY, R)
xr = XY[:,0].reshape(x.shape)
yr = XY[:,1].reshape(y.shape)
return xr, yr
# select and rotate bathy
alpha = 48.
X, Y = rotate(grid['X'], grid['Y'], alpha)
Z = grid['Z'][key]
In [33]:
# interpolate bathy to computational grid and write to file
d = 50
x = np.arange(-500, 1000, d)
y = np.arange(-1500, 2500, d)
Xr, Yr = np.meshgrid(x, y)
np.savetxt(os.path.join(root, 'x.txt'), Xr)
np.savetxt(os.path.join(root, 'y.txt'), Yr)
for i, key in enumerate(sorted(grid['Z'].keys())):
Z = grid['Z'][key]
ix = ~np.isnan(Z.flatten())
Zr = scipy.interpolate.griddata((Y.flatten()[ix],
X.flatten()[ix]),
Z.flatten()[ix],
(Xr.flatten(),
Yr.flatten()), method='nearest').reshape(Xr.shape)
np.savetxt(os.path.join(root, 'z_%d.txt' % i), Zr)
print Zr.shape
In [7]:
# plot computational bathy
fig, axs = subplots()
axs.pcolormesh(Yr, Xr, Zr, cmap='Spectral', vmin=-10, vmax=10)
axs.set_xlim((-1500, 2500))
axs.set_ylim((-500, 1000))
axs.set_aspect('equal')
In [11]:
# read wind time series
meteo, legend = knmi.read_uurgeg(os.path.join(root, 'uurgeg_330_2011-2020.txt'))
ix = (meteo['DD'] > 360.) | (meteo['FH'] >= 900)
meteo = meteo.ix[~ix]
In [51]:
# write wind time series
m = meteo[t0:t1].filter(['FH','DD'])
t = np.asarray([[x.total_seconds() for x in m.index - m.index[0]]]).T
m['FH'] = m['FH'].divide(10.)
m['DD'] = m['DD'].subtract(alpha)
np.savetxt(os.path.join(root, 'wind.txt'), np.concatenate((t, m.as_matrix()), axis=1))
In [14]:
# read and write water levels
wb = waterbase.Waterbase()
wl = wb.get_data('1', 'HOEKVHLD', t0, t1)
t = np.asarray([[(x - wl['datetime'][0]).total_seconds() for x in wl['datetime']]]).T
wl = np.concatenate((t, wl['value'].reshape((-1,1))/100.), axis=1)
np.savetxt(os.path.join(root, 'tide.txt'), wl)
In [18]:
# read and write wave heights
Hs = wb.get_data('22', 'Euro platform', t0, t1)
t = np.asarray([[(x - Hs['datetime'][0]).total_seconds() for x in Hs['datetime']]]).T
Hs = np.concatenate((t, Hs['value'].reshape((-1,1))/100.), axis=1)
np.savetxt(os.path.join(root, 'waves.txt'), Hs)
In [50]:
def create_plot(par, data):
# remove inf's
ix = np.isinf(data)
data[ix] = np.nan
# remove extra dimensions
if data.ndim == 4:
data = data[:,:,0,:] # pick top layer
if data.ndim == 3:
data = data.sum(axis=-1) # sum over fractions
# determine data range
mn = np.nanmin(data[1:-1,1:-1])
mx = np.nanmax(data[1:-1,1:-1])
mx = np.maximum(np.abs(mn), np.abs(mx))
if mn > -1e-10:
mn = 0.
cmap = 'Reds'
else:
mn = -mx
cmap = 'bwr'
# create plot
fig, axs = subplots(figsize=(10,4))
p = axs.pcolormesh(y, x, data, cmap=cmap, vmin=mn, vmax=mx)
axs.set_aspect('equal')
axs.set_title(par)
fig.colorbar(p)
# save plot
fig.savefig(os.path.join(root, 'plots', '%s.png' % par))
n = 10 # number of time steps
with netCDF4.Dataset(os.path.join(root, 'aeolis.nc'), 'r') as ds:
t = ds.variables['time'][:] / 3600.
x = ds.variables['x'][:,:]
y = ds.variables['y'][:,:]
f = ds.variables['fractions'][:]
nx = x.shape[1]
ny = y.shape[0]
for t in np.linspace(1, len(t)-1, n):
t = int(np.round(t))
# bathy difference
data = ds.variables['zb'][...]
create_plot('dz_t=%d' % t, data[t,:,:] - data[0,:,:])
# mean grain size
data = ds.variables['mass'][...]
gs = ds.groups['settings'].getncattr('grain_size') * 1e6
create_plot('d50_t=%d' % t,
np.asarray([np.average(gs, weights=ds.variables['mass'][t,j,i,0,:])
for j in range(ny) for i in range(nx)]).reshape((ny, nx)))
# other variables
for par in ['Ct', 'Cu', 'pickup', 'zb', 'zs']:
create_plot('%s_t=%d' % (par, t), ds.variables[par][t,...])