In [2]:
import netCDF4
from tvtk.api import tvtk
import mayavi.sources.vtk_data_source
import mayavi.mlab
In [3]:
ds = netCDF4.Dataset('/home/fedor/Downloads/subgrid_map_fedor.nc')
In [4]:
xc = ds.variables['FlowElemContour_x'][:]
yc = ds.variables['FlowElemContour_y'][:]
xcc = ds.variables['FlowElem_xcc'][:]
ycc = ds.variables['FlowElem_ycc'][:]
# variables at last timestep
ucx = ds.variables['ucx']
ucy = ds.variables['ucy']
s1 = ds.variables['s1']
#ucx += np.random.random(ucx.shape)
#ucy += np.random.random(ucy.shape)
In [5]:
verts = [np.c_[x,y] for x,y in zip(xc, yc)]
In [6]:
polys = matplotlib.collections.PolyCollection(verts, facecolor=(0.1,0.3,0.2,0.3))
In [7]:
fig, ax = plt.subplots(figsize=(20,13))
ax.add_collection(polys)
ax.autoscale()
In [8]:
grid = tvtk.UnstructuredGrid()
ncells = xc.shape[0]
X = xc
Y = yc
Z = np.zeros_like(X)
pts = np.c_[X.ravel(), Y.ravel(), Z.ravel() ]
In [9]:
cell_array = tvtk.CellArray()
cell_types = np.array([tvtk.Quad().cell_type for i in range(ncells)])
cell_idx = np.array([[4] + [i*4 + j for j in range(4)] for i in range(ncells)]).ravel()
cell_offsets = np.array([i*5 for i in range(ncells)])
cell_array.set_cells(ncells, cell_idx)
In [13]:
grid.points = pts
grid.set_cells(cell_types, cell_offsets, cell_array)
grid.cell_data.scalars = s1[1][1:].copy()
grid.cell_data.scalars.name = 's1'
u, v = ucx[-1], ucy[-1]
vectors = np.c_[u, v, np.zeros_like(u)]
grid.cell_data.vectors = vectors
grid.cell_data.vectors.name = 'uc'
grid.print_traits()
In [13]:
In [14]:
#source = mayavi.sources.vtk_data_source.VTKDataSource(data=grid)
xcc.shape
Out[14]:
In [15]:
# Set a tracer in each cell
pd = tvtk.PolyData()
#pd = tvtk.PointData()
pts = np.c_[xcc, ycc, np.zeros_like(xcc)]
pd.points = pts
pd.points.print_traits()
pd.modified()
# Connect the grid
In [16]:
cell2pointgrid = tvtk.CellDataToPointData()
cell2pointgrid.input = grid
cell2pointgrid.update()
In [17]:
# Create a stream tracer
st = tvtk.StreamTracer()
# Set the points in the streamer
st.source = pd
# Set the corner velocities
st.input = cell2pointgrid.output
In [18]:
# Set some options
# Lookup a nice cell size
cells = grid.get_cells().to_array()
points = grid.points.to_array()
In [19]:
def yieldcells(cells, points):
i = 0
while True:
if i >= cells.shape[0]:
break
n = cells[0]
yield points[cells[i+1:(i+n+1)]]
i += (n+1)
In [20]:
cellsizes = np.array([np.abs(x.max(0) - x.min(0)) for x in yieldcells(cells, points)])
cellsize = np.max(np.mean(cellsizes,0))
cellsize
#LENGTH_UNIT = 1,
#CELL_LENGTH_UNIT = 2
Out[20]:
In [21]:
l = 5000.0 # max propagation
n = 100 # max number of steps
st.maximum_propagation = l # 1km max propagation
st.integration_step_unit = 1
st.minimum_integration_step = (l/n)
st.maximum_integration_step = 10*(l/n)
st.maximum_number_of_steps = n
st.maximum_error = 50.0
st.integrator_type = 'runge_kutta45'
st.print_traits()
In [22]:
st.update()
In [23]:
st.output.print_traits()
In [24]:
def st2linelist(st):
start = 0
linelist = []
lines = st.output.lines.data.to_array()
points = st.output.points.to_array()
for i in range(st.output.lines.number_of_cells):
"""loop over all lines"""
n = lines[start]
idx = lines[start+1:start+n]
line = points[idx]
linelist.append(line.tolist())
start += (n + 1)
linelist = [np.c_[[x[0] for x in line], [x[1] for x in line]] for line in linelist]
return linelist
linelist = st2linelist(st)
In [25]:
fig, ax = plt.subplots(figsize=(20,15))
#ax.add_collection(polys)
#ax.autoscale()
#for vert in verts:
# ax.fill(vert[:,0], vert[:,1], 'b', alpha=0.2)
s = np.s_[:4]
#plt.quiver(xcc[s], ycc[s], ucx[s], ucy[s])
fig.set_facecolor('black')
ax.set_axis_off()
ll = matplotlib.collections.LineCollection(linelist, alpha=0.2, linewidth=4, colors=(0.1,0.2,0.8))
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(linelist, alpha=1, linewidth=0.1, colors=(0.9,1.0,0.8))
ax.add_collection(ll)
ax.autoscale()
In [26]:
def num2deg(xtile, ytile, zoom):
n = 2.0 ** zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = np.arctan(np.sinh(np.pi * (1 - 2 * ytile / n)))
lat_deg = np.degrees(lat_rad)
return (lon_deg, lat_deg)
327,329, 11
791, 793, 11
ll, ur = num2deg(1308, 3167, 13), num2deg(1314, 3163, 13)
In [27]:
import requests
import io
# blackwater styles
images = []
for i in range(1308, 1314):
imagerow = []
for j in range(3163, 3167):
url = 'http://c.tile.cloudmade.com/1a1b06b230af4efdbb989ea99e9841af/124531/256/13/{i}/{j}.png'.format(i=i, j=j)
response = requests.get(url)
img = plt.imread(io.BytesIO(response.content))
imagerow.append(img)
images.append(imagerow)
#url = 'http://b.tile.cloudmade.com/8ee2a50541944fb9bcedded5165f09d9/1/256/15/17599/10746.png'
In [28]:
#img = np.concatenate([np.concatenate(row, axis=0) for row in images], axis=1)
#plt.imsave('2sfo_13_1308_1314_3163_3167.png', img)
img = plt.imread('sfo_13_1308_1314_3163_3167.png')
In [29]:
import osgeo.osr
utm_srs = osgeo.osr.SpatialReference()
utm_srs.ImportFromEPSGA(26910)
osm_srs = osgeo.osr.SpatialReference()
osm_srs.ImportFromEPSGA(3857)
wgs_srs = osgeo.osr.SpatialReference()
wgs_srs.ImportFromEPSGA(4326)
utm2wgs = osgeo.osr.CoordinateTransformation(utm_srs, wgs_srs)
utm2osm = osgeo.osr.CoordinateTransformation(utm_srs, osm_srs)
wgs2utm = osgeo.osr.CoordinateTransformation(wgs_srs, utm_srs)
wgs2osm = osgeo.osr.CoordinateTransformation(wgs_srs, osm_srs)
ll_x = wgs2osm.TransformPoint(ll[0], ll[1])
ur_x = wgs2osm.TransformPoint(ur[0], ur[1])
ll_x, ur_x
Out[29]:
In [33]:
x,y,_ = wgs2utm.TransformPoint(-122.4239052, 37.8279125)
# create seeds starting from this location
n_seeds = 100
seeds = np.c_[x + np.random.random(n_seeds) * 100, y + np.random.random(n_seeds) * 100, np.zeros(n_seeds)]
pd.points = seeds
grid.cell_data.vectors = np.c_[ucx[-10], ucy[-10], np.zeros_like(ucx[-10])]
st.update()
linelist_0 = st2linelist(st)
grid.cell_data.vectors = np.c_[ucx[-1], ucy[-1], np.zeros_like(ucx[-1])]
st.update()
linelist_1 = st2linelist(st)
In [34]:
transformed_lines_0 = [np.array(utm2osm.TransformPoints(line))[:,:2] for line in linelist_0]
transformed_lines_1 = [np.array(utm2osm.TransformPoints(line))[:,:2] for line in linelist_1]
In [35]:
fig, ax = plt.subplots(figsize=(20,20))
ax.imshow(img, extent=(ll_x[0], ur_x[0], ll_x[1], ur_x[1]))
col = matplotlib.collections.LineCollection(transformed_lines_0, alpha=0.2, linewidth=1, colors=(0.1,0.2,0.8))
ax.add_collection(col)
col = matplotlib.collections.LineCollection(transformed_lines_0, alpha=0.3, linewidth=0.1, colors='white')
ax.add_collection(col)
col = matplotlib.collections.LineCollection(transformed_lines_1, alpha=0.2, linewidth=1, colors=(0.1,0.8,0.2))
ax.add_collection(col)
col = matplotlib.collections.LineCollection(transformed_lines_1, alpha=0.3, linewidth=0.1, colors='white')
ax.add_collection(col)
plt.plot(x,y, 'rx', alpha=1)
ax.set_xlim(ll_x[0]+6000, ur_x[0]-12000)
ax.set_ylim(ll_x[1]+6000, ur_x[1]-8000)
Out[35]:
In [ ]: