In [2]:
import netCDF4
from tvtk.api import tvtk
import mayavi.sources.vtk_data_source
import mayavi.mlab
import matplotlib.cm
import matplotlib.colors
import osgeo.osr
import datetime
import pandas
import scipy.interpolate
import requests
import io
import collections
import python_subgrid.particles
import logging
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
logging.basicConfig()
logging.root.setLevel(logging.DEBUG)
import pandas
import netCDF4
In [3]:
df = pandas.read_hdf('/Users/baart_f/models/sfo/sfo-3di/parts6.h5', 'particles')
In [4]:
mins = df[['particle', 't']].groupby('particle').min()
In [5]:
import dateutil.rrule
rrule = dateutil.rrule.rrule(dateutil.rrule.HOURLY, dtstart=datetime.datetime(1962, 6, 11, 18, 0))
escapes = rrule.between(datetime.datetime(1962, 6, 11, 20, 0), datetime.datetime(1962, 6, 12, 4, 0), inc=True)
In [6]:
ds = netCDF4.Dataset('/Users/baart_f/models/sfo/sfo-3di/subgrid_map_15min.nc')
times = netCDF4.num2date(ds.variables['time'][:], ds.variables['time'].units)
idx = np.where(np.in1d(times, escapes))[0]
idx
Out[6]:
In [7]:
i0 = 250
swimmers = mins[np.in1d(mins['t'],ds.variables['time'][:][idx-i0])].index
t_swimmers = ds.variables['time'][:][idx-i0]
N = matplotlib.colors.Normalize(0, 90*900)
C = matplotlib.cm.Accent
In [8]:
# generate the spatial reference systems and the transformations
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)
In [9]:
def num2deg(xtile, ytile, zoom):
"""convert x,y zoom number to a lat/long"""
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)
# we're using this area
ll, ur = num2deg(1308, 3167, 13), num2deg(1314, 3163, 13)
# lookup bounding box in different coordinate systems
ll_osm = wgs2osm.TransformPoint(ll[0], ll[1])
ur_osm = wgs2osm.TransformPoint(ur[0], ur[1])
ll_utm = wgs2utm.TransformPoint(ll[0], ll[1])
ur_utm = wgs2utm.TransformPoint(ur[0], ur[1])
ll, ur, ll_osm, ur_osm, ll_utm, ur_utm
Out[9]:
In [10]:
# use this image as a background
# generated with mapbox
img = plt.imread('sfo_13_1308_1314_3163_3167-sepia.png')
2160/5.83
Out[10]:
In [10]:
# print on a4 size
fig, ax = plt.subplots(figsize=(8.27*np.sqrt(2), 5.83*np.sqrt(2)))
# we're shifting time a bit, going from index based time to seconds
# nut sure why I use 1 and not 0 here...
t0 = 1 # 270
t1 = 90
# give the particles a trail of half an our
trail = 1800
# group dataset by particles
grouped = df.groupby('particle')
# Loop over time, 1 frame per minute
for i, t in enumerate(np.arange(t0*900, t1*900, 60)):
# cleanup the image
ax.clear()
# big title
ax.title.set_size(36)
# show background map
ax.imshow(img, extent=(ll_osm[0], ur_osm[0], ll_osm[1], ur_osm[1]))
# no need for axis, we know where we are
ax.axis('off')
# set title to current time
ax.set_title(times[i0] + datetime.timedelta(seconds=t))
# limit axis, we're in google mercator coordinates here.
ax.set_xlim(ll_osm[0]+2000, ur_osm[0]-3000)
ax.set_ylim(ll_osm[1]+2000, ur_osm[1]-1000)
# lookup values with a trail of length trail
selection = df[np.logical_and(df['t'] <= t, df['t'] > t - trail)]
ends = df[df['t'] == t]
if selection.empty:
continue
# generate lines for particles and swimmers, color the swimmers
linelist_osm = []
swimmers_osm = []
swimmer_colors = []
# process all particles
for j, group in selection.groupby('particle'):
# convert to google mercator
x, y, _ = np.array(utm2osm.TransformPoints(np.c_[group['x'], group['y']])).T
coords = np.c_[x,y]
# use two lists, one for swimmers, one for other particles
if j in swimmers:
swimmers_osm.append(coords)
# get color based on water time
swimmer_colors.append(C(N(grouped.get_group(j)['t'].irow(0))))
else:
linelist_osm.append(coords)
if swimmers_osm:
ll = matplotlib.collections.LineCollection(swimmers_osm, alpha=0.25, linewidth=4, colors=swimmer_colors, transOffset=ax.transData)
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(swimmers_osm, alpha=0.3, linewidth=0.25, colors=(0.9,1.0,0.8), transOffset=ax.transData)
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(linelist_osm, alpha=0.25, linewidth=4, colors=(0.1,0.2,0.8), transOffset=ax.transData)
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(linelist_osm, alpha=1, linewidth=0.25, colors=(0.9,1.0,0.8), transOffset=ax.transData)
ax.add_collection(ll)
# add the white dots
pts = np.c_[ends['x'], ends['y']]
x, y, _ = np.array(utm2osm.TransformPoints(pts)).T
ax.plot(x, y, 'w.', alpha=0.5, markeredgecolor='none')
# get rid of border
fig.subplots_adjust(left=0, bottom=0, right=1, top=0.95)
# store (use ffmpeg to generate animation)
fig.savefig('merged_%05d.png' % (i,)) # dpi=370 for a5 == max res at youtube
In [15]:
fig, ax = plt.subplots(figsize=(21.9882, 21.9882/np.sqrt(2)))
t0 = 1 # 270
t1 = 90
trail = 1800
grouped = df.groupby('particle')
ts = np.arange(15, 15+15)*1800
N = matplotlib.colors.Normalize(3600*6, 3600*13)
C = matplotlib.cm.jet
font = {'family' : 'Myriad Pro',
'weight' : 'normal',
'size' : 36,
}
for i, t in enumerate(ts):
ax.clear()
ax.imshow(img, extent=(ll_osm[0], ur_osm[0], ll_osm[1], ur_osm[1]))
ax.axis('off')
ax.set_title(times[i0] + datetime.timedelta(seconds=t), fontdict=font)
ax.set_xlim(ll_osm[0]+2000, ur_osm[0]-3000)
ax.set_ylim(ll_osm[1]+2000, ur_osm[1]-1000)
# lookup values with a trail of length trail
selection = df[np.logical_and(df['t'] <= t, df['t'] > t - trail)]
ends = df[df['t'] == t]
if selection.empty:
continue
linelist_osm = []
swimmers_osm = []
swimmer_colors = []
for j, group in selection.groupby('particle'):
x, y, _ = np.array(utm2osm.TransformPoints(np.c_[group['x'], group['y']])).T
coords = np.c_[x,y]
if j in swimmers:
swimmers_osm.append(coords)
# get color based on water time
swimmer_colors.append(C(N(grouped.get_group(j)['t'].irow(0))))
else:
linelist_osm.append(coords)
if swimmers_osm:
ll = matplotlib.collections.LineCollection(swimmers_osm, alpha=0.25, linewidth=4, colors=swimmer_colors, transOffset=ax.transData)
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(swimmers_osm, alpha=0.3, linewidth=0.25, colors=(0.9,1.0,0.8), transOffset=ax.transData)
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(linelist_osm, alpha=0.25, linewidth=4, colors=(0.1,0.2,0.8), transOffset=ax.transData)
ax.add_collection(ll)
ll = matplotlib.collections.LineCollection(linelist_osm, alpha=1, linewidth=0.25, colors=(0.9,1.0,0.8), transOffset=ax.transData)
ax.add_collection(ll)
pts = np.c_[ends['x'], ends['y']]
x, y, _ = np.array(utm2osm.TransformPoints(pts)).T
ax.plot(x, y, 'w.', alpha=0.5, markeredgecolor='none')
fig.subplots_adjust(left=0, bottom=0, right=1, top=0.95)
fig.savefig('lenticulair_%05d.png' % (i,), dpi=300) # dpi=370 for a5 == max res at youtube
fig.savefig('lenticulair_%05d.pdf' % (i,)) # dpi=370 for a5 == max res at youtube