In [1]:
%matplotlib inline
import netCDF4
import matplotlib.pyplot as plt
import numpy as np
In [2]:
ds = netCDF4.Dataset('/Users/baart_f/data/delft3d/WinterCycle3_DisSeabed_TSS.nc')
In [3]:
import shapely.affinity
import shapely.geometry
In [4]:
var = ds.variables['suspsedconc']
depth = ds.variables['depth'][1:-1,1:-1]
waterlevel = ds.variables['waterlevel'][0, 1:-1, 1:-1]
x = ds.variables['x'][1:-1,1:-1]
y = ds.variables['y'][1:-1,1:-1]
grid_x = ds.variables['grid_x'][1:-1,1:-1,:]
grid_y = ds.variables['grid_y'][1:-1,1:-1,:]
layer_interf = ds.variables['LayerInterf'][:].squeeze()
layer = ds.variables['Layer'][:].squeeze()
# sigma positive upward
Z = -np.tile(layer_interf, x.shape + (1,)) \
* np.repeat((depth + waterlevel)[:,:,np.newaxis], layer_interf.shape[0], 2) \
- np.repeat(depth[:,:,np.newaxis], layer_interf.shape[0],2)
# total water depth
# total water height from bottom
# convert from structured (curvilinear) to rectilinear grid
pts = [shapely.geometry.Point(x_, y_) for x_, y_ in np.c_[x[0], y[0]]]
distances_x = -np.array([pts[0].distance(pts[i]) for i in range(len(pts))])
pts = [shapely.geometry.Point(x_, y_) for x_, y_ in np.c_[x[:,0], y[:,0]]]
distances_y = np.array([pts[0].distance(pts[i]) for i in range(len(pts))])
distances_x.shape, distances_y.shape, x.shape, depth.shape, Z.shape
Out[4]:
In [5]:
fig, axes = plt.subplots(1,3, sharex=True, figsize=(16,5))
ax0, ax1, ax2 = axes
pc = ax0.pcolormesh(x, y, Z[...,0])
plt.colorbar(pc, ax=ax0)
ax0.set_title('depth layer 0')
pc = ax1.pcolormesh(x, y, Z[...,15])
ax1.set_title('depth layer 15')
plt.colorbar(pc, ax=ax1)
# top interface should follow sealevel
ax2.set_title('depth layer -1')
pc = ax2.pcolormesh(x, y, Z[:-1,:-1,-1])
plt.colorbar(pc, ax=ax2)
fig.tight_layout()
In [6]:
data = var[60,...]
data = data.squeeze()[:,1:-1,1:-1]
# roll the first dimension to the end
data = np.rollaxis(data, 0, 3)
# data[data>=0.00015] = 0.00016
#data[data<0] = 0.0
#data[data<=0.00001] = 0.00001
#data = np.log10(data)
x.shape, data[0].shape
Out[6]:
In [69]:
centers = np.c_[x.ravel(), y.ravel()]
x_ll = x - np.gradient(x)[0]/2.0 - np.gradient(x)[1]/2.0
x_lr = x - np.gradient(x)[0]/2.0 + np.gradient(x)[1]/2.0
x_ur = x + np.gradient(x)[0]/2.0 + np.gradient(x)[1]/2.0
x_ul = x + np.gradient(x)[0]/2.0 - np.gradient(x)[1]/2.0
y_ll = y - np.gradient(y)[0]/2.0 - np.gradient(y)[1]/2.0
y_lr = y - np.gradient(y)[0]/2.0 + np.gradient(y)[1]/2.0
y_ur = y + np.gradient(y)[0]/2.0 + np.gradient(y)[1]/2.0
y_ul = y + np.gradient(y)[0]/2.0 - np.gradient(y)[1]/2.0
cell_x = np.dstack([x_ll, x_lr, x_ur, x_ul])
cell_y = np.dstack([y_ll, y_lr, y_ur, y_ul])
nodes_x = np.empty((x.shape[0]+1, x.shape[1]+1), dtype='double')
nodes_y = np.empty_like(nodes_x)
nodes_x[:-1,:-1] = x_ul
nodes_x[-1,:-1] = x_ll[-1,:]
nodes_x[:-1,-1] = x_ur[:,-1]
nodes_x[-1,-1] = x_lr[-1,-1]
nodes_y[:-1,:-1] = y_ul
nodes_y[-1,:-1] = y_ll[-1,:]
nodes_y[:-1,-1] = y_ur[:,-1]
nodes_y[-1,-1] = y_lr[-1,-1]
#nodes_x = nodes_x.ravel()
#nodes_y = nodes_y.ravel()
polys = []
n_rows = x.shape[0]
n_cols = x.shape[1]
n_cells = np.prod(x.shape)
print nodes_x.shape, n_cells
maxidx =0
t maxidx=max(maxidx, max(*idx))
x_i = nodes_x.ravel()[idx]
y_i = nodes_y.ravel()[idx]
xy_i = np.c_[x_i,y_i]
polys.append(xy_i)
print maxidx
import matplotlib.collections
patches = matplotlib.collections.PatchCollection((matplotlib.patches.Polygon(xy, closed=True) for xy in polys), alpha=0.5, facecolor='none', edgecolor='black')
fig, ax = plt.subplots()
ax.add_collection(patches)
ax.autoscale()
ax.set_ylim(5162000, 5167000)
ax.set_xlim(690000, 695000)
Out[69]:
In [101]:
plt.pcolormesh(distances_x, distances_y, data[:-1, :-1, -10])
plt.colorbar()
Out[101]:
In [8]:
import mayavi.mlab as mlab
from tvtk.api import tvtk
In [55]:
pts = np.c_[X.ravel(), Y.ravel(), Z.ravel()]
grid = tvtk.
grid.print_traits()
grid.point_data.scalars = data.ravel()
grid.point_data.scalars.name = 'sediment'
dataset = mlab.pipeline.add_dataset(grid)
gx = mlab.pipeline.grid_plane(dataset)
gy = mlab.pipeline.grid_plane(dataset)
gy.grid_plane.axis = 'y'
gz = mlab.pipeline.grid_plane(dataset)
gz.grid_plane.axis = 'z'
iso = mlab.pipeline.iso_surface(dataset)
iso.contour.maximum_contour = 0.001
vol = mlab.pipeline.volume(grid)
In [9]:
mlab.show()
In [58]:
mlab.contour3d(np.swapaxes(np.swapaxes(data, 0,2), 0,1), extent=[0, x.shape[0], 0, x.shape[1],0,5])
mlab.surf(depth, extent=[0, x.shape[0], 0, x.shape[1],0,5])
mlab.outline()
engine = mlab.get_engine()
scene = mlab.gcf()
from mayavi.filters.point_to_cell_data import PointToCellData
point_to_cell_data = PointToCellData()
array_source = engine.scenes[0].children[0]
engine.add_filter(point_to_cell_data, array_source)
from mayavi.modules.volume import Volume
volume = Volume()
volume.print_traits()
engine.add_filter(volume, point_to_cell_data)
mlab.outline()
#volume.update_ctf = 0
Out[58]:
In [59]:
mlab.show()
In [43]:
X.min(), X.max(), Y.min(), Y.max(), Z.min(), Z.max()
Out[43]:
In [63]:
import IPython.display
IPython.display.Image('plume.png')
Out[63]:
In [ ]: