In [1]:
import pandas
import datetime
import numpy as np
import netCDF4
import osgeo.osr
import tables
import shapely.geometry
import shapely.prepared
import matplotlib.pyplot as plt
import dateutil.rrule
%matplotlib inline
In [2]:
# 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 [3]:
import simplejson
# load the USA topography
sfo = simplejson.load(open('sfo.geojson'))
polygons = [shapely.geometry.asShape(feature['geometry']) for feature in sfo['features']]
polygons = [polygon[0] if len(polygon) == 1 else polygon for polygon in polygons]
sfo = shapely.geometry.MultiPolygon(polygons)
In [4]:
import shapely.geometry
#USA = shapely.geometry.Polygon(shell=np.c_[lon[min:max], lat[min:max]][1:])
USA = sfo
box = shapely.geometry.box(-122.8, 37.4 ,-122, 38)
bay = box.intersection(USA)
In [5]:
bay_utms = []
for poly in bay:
xy = np.array(poly.exterior)
xy = np.array(wgs2utm.TransformPoints(xy))[:,:2]
bay_utms.append(shapely.geometry.Polygon(xy))
bay_utm = shapely.geometry.MultiPolygon(polygons=bay_utms)
bay_utm_without_alcatraz = shapely.geometry.MultiPolygon(
polygons=[poly for i, poly in enumerate(bay_utms) if i != 5]
)
bay_utm_buffered = bay_utm_without_alcatraz.buffer(500)
bay_utm_buffered_prepped = shapely.prepared.prep(bay_utm_buffered)
for i, poly in enumerate(bay_utm):
xy = np.array(poly.exterior)
plt.fill(xy[:,0], xy[:,1], facecolor=(0,0,0,0.3), edgecolor='none')
plt.text(poly.centroid.x, poly.centroid.y, str(i))
for poly in bay_utm_buffered:
xy = np.array(poly.exterior)
plt.fill(xy[:,0], xy[:,1], facecolor='none', edgecolor='black')
plt.xlim(535000, 557000)
plt.ylim(4180000, 4200000)
labels = [
(bay_utms[2].centroid, 'Angel Island'),
(bay_utms[5].centroid, 'Alcatraz')
]
In [6]:
filename = '/Users/baart_f/src/python-subgrid/notebooks/parts.h5'
name = 'particles_19'
df = pandas.read_hdf(filename, name)
df['particle'] = df['particle'].astype('int64')
df['t'] = df['t'].astype('int64')
In [19]:
# inspect the particle file
f = tables.open_file(filename)
# lookup the table/node with name
node = f.get_node('/' + name)
attrs = node._v_attrs
# show the details
print node, {name:getattr(attrs,name) for name in attrs._v_attrnames if name in ('u', 'v', 'velocity', 'speed')}
# show all the nodes
root = f.root
for name_ in f.root._v_children:
node = f.get_node('/' + name_)
attrs = node._v_attrs
print node, {name_: getattr(attrs, name_) for name_ in attrs._v_attrnames if name_ in ('u', 'v', 'velocity', 'speed')}
f.close()
del f
In [8]:
# define the time rule when swimmers started
#rrule = dateutil.rrule.rrule(dateutil.rrule.HOURLY, dtstart=datetime.datetime(1962, 6, 11, 19, 0), interval=1)
rrule = dateutil.rrule.rrule(dateutil.rrule.MINUTELY, dtstart=datetime.datetime(1962, 6, 11, 19, 0), interval=30)
# given the rule we can define the escape times
escapes = rrule.between(datetime.datetime(1962, 6, 11, 20, 0), datetime.datetime(1962, 6, 12, 4, 0), inc=True)
# lookup the indices in the 2d result file
ds = netCDF4.Dataset('/Users/baart_f/models/sfo/sfo-3di/subgrid_map_15min.nc')
times = netCDF4.num2date(ds.variables['time'][:], ds.variables['time'].units)
# convert to indices
idx = np.where(np.in1d(times, escapes))[0]
idx
Out[8]:
In [9]:
# we use t_idx=250 as a reference time (that's where the particle simulation starts)
i0 = 250
# lookup the start time per particle
mins = df[['particle', 't']].groupby('particle').min()
# if the start time is a swimmer time, we have a swimmer
swimmers = mins[np.in1d(mins['t'],ds.variables['time'][:][idx-i0])].index
# define the real times where these started
t_swimmers = ds.variables['time'][:][idx-i0]
In [10]:
# extract the dataset with swimmers
swimming = df[np.in1d(df['particle'], swimmers)].copy()
# convert the times
swimming['date'] = np.fromiter([times[i0] + datetime.timedelta(seconds=x) for x in swimming['t']], dtype='datetime64[s]')
In [11]:
# convert to lat/lon
latlonz = utm2wgs.TransformPoints(np.c_[swimming['x'], swimming['y']])
swimming['lon'], swimming['lat'], _ = np.array(latlonz).T
In [12]:
# Lookup the start time
mins['start'] = np.fromiter([times[i0] + datetime.timedelta(seconds=x) for x in mins['t']], dtype='datetime64[s]')
# merge the start time per particle
results = pandas.merge(swimming, mins['start'].reset_index(), left_on='particle', right_on='particle')
# add the start hour (as a label)
results['hour'] = results['start'].apply(lambda x:x.strftime('%H:%M'))
# compute the time since escape
results['dt'] = np.subtract(results['t'], [x.total_seconds() for x in results['start'] - times[i0]])
In [13]:
# store for cartodb
results[['lat', 'lon', 'date', 'hour', 'reason']].to_csv('/Users/baart_f/models/sfo/sfo-3di/parts8a.csv')
In [14]:
results.head()
Out[14]:
In [15]:
# add death condition
import itertools
max_alive = 60*60*5.5
df_death = results.copy()
# die after 5 hours (reason 11)
reasons = []
sunset = datetime.datetime(1962, 6, 12, 5, 47)
counter = itertools.count()
for id_, particle in df_death.groupby(['particle']):
if counter.next() > 10:
pass
# assert sorted in t
if not list(particle['t']) == list(sorted(particle['t'])):
raise ValueError(str(particle))
reason = {'fate': 'unknown', 'row': None, 'particle': id_}
for i, row in particle.iterrows():
if row.reason == 1:
reason['fate'] = 'ocean/bay'
reason['row'] = i
break
elif row['dt'] >= 60*60*5:
reason['fate'] = 'exhausted'
reason['row'] = i
break
elif row['date'] > sunset:
reason['fate'] = 'caught'
reason['row'] = i
break
elif bay_utm_buffered_prepped.contains(shapely.geometry.Point(row.x, row.y)):
reason['fate'] = 'escaped'
reason['row'] = i
break
reasons.append(reason)
#df_death.loc[np.logical_and(results['dt'] == 60*60*5, results['reason']==4),'reason'] = 11
#df_death.loc[np.logical_and(results['dt'] == 60*60*5, results['reason']==4),'reason'] = 11
#df_death.loc[np.logical_and(results['date'] > datetime.datetime(1962, 6, 12, 5,0), results['reason']==4),'reason'] = 11
# VTK termination reasons
OUT_OF_DOMAIN = 1
NOT_INITIALIZED = 2
UNEXPECTED_VALUE = 3
OUT_OF_TIME = 4
OUT_OF_STEPS = 5
STAGNATION = 6
# this is extra for living particles
DEATH = 11
ESCAPE = 12
CAUGHT = 13
reasons_df = pandas.DataFrame(data=reasons)
In [16]:
import matplotlib.colors
_data = (
(0.49803921568627452, 0.78823529411764703, 0.49803921568627452),
(0.74509803921568629, 0.68235294117647061, 0.83137254901960789),
(0.99215686274509807, 0.75294117647058822, 0.52549019607843139),
(0.2196078431372549, 0.42352941176470588, 0.69019607843137254),
(0.94117647058823528, 0.00784313725490196, 0.49803921568627452),
(0.74901960784313726, 0.35686274509803922, 0.09019607843137254),
(1.0, 1.0, 0.6 ),
(0.4, 0.4, 0.4 ),
)
accent = matplotlib.colors.ListedColormap(_data[:4], 'accent')
In [17]:
#plt.plot(df_death['x'],df_death['y'], 'k.')
endings = pandas.merge(df_death, reasons_df, left_index=True, right_on='row')
fates = ['escaped', 'caught', 'exhausted', 'ocean/bay']
colors = [(fates.index(x)) for x in endings['fate']]
accent.N = len(fates)
fig, ax = plt.subplots(figsize=(13,8))
label = 'San Francisco'
for poly in bay_utm:
xy = np.array(poly.exterior)
ax.fill(xy[:,0], xy[:,1], facecolor=(0,0,0,0.3), edgecolor='none', label=label)
label = None
label='escape buffer (500m)'
for poly in bay_utm_buffered:
xy = np.array(poly.exterior)
ax.fill(xy[:,0], xy[:,1], facecolor='none', edgecolor='black', label=label)
# only one label
label = None
sc = ax.scatter(endings['x'], endings['y'], c=colors, edgecolor='none', s=30, cmap=accent, alpha=0.5)
ax.plot(endings[endings['fate']=='escaped']['x'],
endings[endings['fate']=='escaped']['y'],
'gx', markersize=0)
cbar = plt.colorbar(sc, ax=ax)
cbar.ax.get_yaxis().set_ticks([])
for j, lab in enumerate(fates):
cbar.ax.text(1.5, (2 * j + 1) / 8.0, lab, va='center')
ax.set_xlabel('x coordinate [m] (NAD83/UTM10N)')
ax.set_ylabel('y coordinate [m] (NAD83/UTM10N)')
ax.axis('equal')
plt.legend()
ax.set_xlim(535000, 557000)
ax.set_ylim(4180000, 4200000)
fig.savefig('plot_%s.pdf' % (name,))
In [18]:
endings
cross = endings[['fate', 'hour', 'particle_x']].pivot_table(columns=['fate'], index=['hour'], aggfunc='count')
cross.to_latex(open('table_%s.tex' % (name,), 'w'))
In [18]:
In [18]: