In [107]:
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
%matplotlib inline
In [108]:
# Create the coastline shape
import netCDF4
#ds = netCDF4.Dataset('http://opendap.deltares.nl/thredds/dodsC/opendap/noaa/gshhs/gshhs_f.nc')
ds = netCDF4.Dataset('/Users/baart_f/Downloads/gshhs_f.nc')
lat = ds.variables['lat'][:]
lon = ds.variables['lon'][:]
In [109]:
# Lookup the correct continent, ignore lakes and stuff
idx = np.logical_and.reduce([lon >= -122.8, lon <= -122, lat >= 37.4 , lat < 38])
plt.plot(lon[idx], lat[idx], 'k-')
nanidx = np.where(np.isnan(lon))[0]
#min, max = np.where(idx)[0].min(), np.where(idx)[0].max()
min, max = np.where(idx)[0][0], np.where(idx)[0][0]
minmin = nanidx[nanidx < min].max()
maxmax = nanidx[nanidx > max ].min()
minmin, maxmax, min, max
min, max = 1455261, 2674864
# usa
min, max
plt.plot(lon[min:max], lat[min:max], 'g-')
Out[109]:
In [110]:
# 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 [111]:
import shapely.geometry
USA = shapely.geometry.Polygon(shell=np.c_[lon[min:max], lat[min:max]][1:])
box = shapely.geometry.box(-122.8, 37.4 ,-122, 38)
bay = box.intersection(USA)
In [198]:
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_buffered = bay_utm.buffer(400)
bay_utm_buffered_prepped = shapely.prepared.prep(bay_utm_buffered)
bay_buffereds = []
for poly in bay_utm_buffered:
xy = np.array(poly.exterior)
xy = np.array(utm2wgs.TransformPoints(xy))[:,:2]
bay_buffereds.append(shapely.geometry.Polygon(xy))
bay_buffered = shapely.geometry.MultiPolygon(polygons=bay_buffereds)
bay_buffered_prepped = shapely.prepared.prep(bay_buffered)
fig, (ax1, ax2) = plt.subplots(1,2)
for poly in bay_utm:
xy = np.array(poly.exterior)
ax1.plot(xy[:,0], xy[:,1], 'k-')
for poly in bay_utm_buffered:
xy = np.array(poly.exterior)
ax1.plot(xy[:,0], xy[:,1], 'g-')
for poly in bay:
xy = np.array(poly.exterior)
ax2.plot(xy[:,0], xy[:,1], 'k-')
for poly in bay_buffered:
xy = np.array(poly.exterior)
ax2.plot(xy[:,0], xy[:,1], 'g-')
In [199]:
df = pandas.read_csv('/Users/baart_f/Downloads/parts4a.csv')
df = pandas.read_hdf('/Users/baart_f/src/siggyf-python-subgrid/notebooks/parts5.h5', 'particles')
df.head()
len(df)
Out[199]:
In [200]:
import dateutil.rrule
rrule = dateutil.rrule.rrule(dateutil.rrule.HOURLY, dtstart=datetime.datetime(1962, 6, 11, 19, 0)) #, interval=30)
escapes = rrule.between(datetime.datetime(1962, 6, 11, 20, 0), datetime.datetime(1962, 6, 12, 4, 0), inc=True)
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[200]:
In [201]:
i0 = 250
mins = df[['particle', 't']].groupby('particle').min()
swimmers = mins[np.in1d(mins['t'],ds.variables['time'][:][idx-i0])].index
t_swimmers = ds.variables['time'][:][idx-i0]
In [202]:
swimming = df[np.in1d(df['particle'], swimmers)]
swimming['date'] = np.fromiter([times[i0] + datetime.timedelta(seconds=x) for x in swimming['t']], dtype='datetime64[s]')
In [203]:
latlonz = utm2wgs.TransformPoints(np.c_[swimming['x'], swimming['y']])
swimming['lon'], swimming['lat'], _ = np.array(latlonz).T
In [204]:
mins['start'] = np.fromiter([times[i0] + datetime.timedelta(seconds=x) for x in mins['t']], dtype='datetime64[s]')
results = pandas.merge(swimming, mins['start'].reset_index(), left_on='particle', right_on='particle')
results['hour'] = results['start'].apply(lambda x:x.strftime('%H:%M'))
results['dt'] = np.subtract(results['t'], [x.total_seconds() for x in results['start'] - times[i0]])
results = results[results['dt'] <= 60*60*5]
In [205]:
pts = shapely.geometry.MultiPoint(points=np.array(results[['lon', 'lat']]))
on_land = [bay_buffered_prepped.contains(pt) for pt in pts]
results['on_land'] = on_land
any(on_land)
survive_dict = results[['on_land','particle']].groupby('particle').any()
In [206]:
results2 = results.copy()
results2['reason'].ix[np.logical_and(results['dt'] == 60*60*5, results['reason']==4)] = 11
a= results2.groupby(['particle'])['dt'].max()
records = []
for particle, tstop in a.iteritems():
record = results2[np.logical_and(results2['particle'] == particle, results2['dt']== tstop)]
survive = survive_dict.ix[particle]
if survive.on_land:
record['reason'] = 12
records.append(record.irow(0))
reasons = pandas.DataFrame.from_records(records)[['reason', 'hour', 'particle']]
survival = reasons.pivot_table(rows=['reason'], cols=['hour'], aggfunc='count')
print(survival)
In [207]:
plt.subplots(figsize=(13,8))
for poly in bay:
xy = np.array(poly.exterior)
plt.plot(xy[:,0], xy[:,1], 'k-')
for poly in bay_buffered:
xy = np.array(poly.exterior)
plt.plot(xy[:,0], xy[:,1], 'g-', linewidth=0.1)
plt.plot(results['lon'], results['lat'], 'b.', linewidth=0.1, alpha=0.5, markersize=0.1)
for i in reasons[reasons['reason'] == 3].particle:
a = results[results['particle'] == i]
plt.plot(a['lon'], a['lat'], 'k.', linewidth=0.1, markersize=0.1)
plt.savefig('a.pdf')
In [100]:
In [ ]: