In [26]:

    

import io
import collections
import logging
import datetime

import numpy as np
import pandas
import scipy.interpolate
import requests

import matplotlib.cm
import matplotlib.colors
import matplotlib.pyplot as plt
import matplotlib.style
matplotlib.style.use('grayscale')
from tvtk.api import  tvtk
import mayavi.sources.vtk_data_source 
import mayavi.mlab

import osgeo.osr
import netCDF4

import python_subgrid.particles

logging.basicConfig()
logging.root.setLevel(logging.DEBUG)

%matplotlib inline


    

In [27]:

    

# 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 [28]:

    

# define the location of the Alcatraz escape location
alcatraz = {}
alcatraz['lon'] = -122.4239052
alcatraz['lat'] = 37.8279125
alcatraz['x_osm'], alcatraz['y_osm'], _ = wgs2osm.TransformPoint(alcatraz['lon'], alcatraz['lat'])
alcatraz['x_utm'], alcatraz['y_utm'], _ = wgs2utm.TransformPoint(alcatraz['lon'], alcatraz['lat'])


    

In [29]:

    

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[29]:





((-122.51953125, 37.75334401310657),
 (-122.255859375, 37.892195547244356),
 (-13638811.830980573, 4544639.953723437, 0.0),
 (-13609460.012119064, 4564207.832964446, 0.0),
 (542324.8953110904, 4178557.256203641, 0.0),
 (565429.5404843381, 4194114.96778207, 0.0))

In [30]:

    

# I don't know a good way to use a TMS service as a background
# Download blackwater styles

images = []
for i in range(1308, 1314):
    imagerow = []
    for j in range(3163, 3167):
        # don't download now, remove this if you want to redownload images
        continue
        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)
# This concatenates all images and saves it as a file    
#img = np.concatenate([np.concatenate(row, axis=0) for row in images], axis=1)
#plt.imsave('2sfo_13_1308_1314_3163_3167.png', img)


    

In [31]:

    

# we'll use the version with a sepia color (for contrast with the blue lines)
img = plt.imread('sfo_13_1308_1314_3163_3167-sepia.png')
plt.imshow(img, extent=(ll_osm[0], ur_osm[0], ll_osm[1], ur_osm[1]))


    

    
Out[31]:





<matplotlib.image.AxesImage at 0x11901ca90>






    

In [33]:

    

#ds2 = netCDF4.Dataset('/home/fedor/Checkouts/models/sanfrancisco/subgrid_map.nc')
ds = netCDF4.Dataset('/home/fedor/Downloads/subgrid_map_15min.nc')


    

    




---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-33-cb74a31c46d5> in <module>()
      1 #ds2 = netCDF4.Dataset('/home/fedor/Checkouts/models/sanfrancisco/subgrid_map.nc')
----> 2 ds = netCDF4.Dataset('/home/fedor/Downloads/subgrid_map_15min.nc')

/Users/baart_f/.virtualenvs/main/lib/python2.7/site-packages/netCDF4.so in netCDF4.Dataset.__init__ (netCDF4.c:19738)()

RuntimeError: No such file or directory

In [8]:

    

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']


    

In [9]:

    

times = netCDF4.num2date(ds.variables['time'][:], ds.variables['time'].units)
datetime.timedelta(hours=16)


    

    
Out[9]:





datetime.timedelta(0, 57600)

In [10]:

    

import bisect
bisect.bisect(times, datetime.datetime(1962,6,11, 22,0)), bisect.bisect(times, datetime.datetime(1962,6,12, 8,0))
#times[472]
#times[512]


    

    
Out[10]:





(281, 321)

In [11]:

    

verts = [np.array(utm2osm.TransformPoints(np.c_[x,y]))[:,:2] for x,y in zip(xc, yc)]


    

In [12]:

    

fig, ax = plt.subplots(figsize=(20,13))
polys = matplotlib.collections.PolyCollection(verts, facecolor=(0.1,0.3,0.2,0.1), linewidth=0.1, transOffset=ax.transData)
ax.imshow(img, extent=(ll_osm[0], ur_osm[0], ll_osm[1], ur_osm[1]))
ax.add_collection(polys)
ax.set_xlim(ll_osm[0]-50000, ur_osm[0]+50000)
ax.set_ylim(ll_osm[1]-70000, ur_osm[1]+70000)
ax.set_xlim(ll_osm[0], ur_osm[0])
ax.set_ylim(ll_osm[1], ur_osm[1])

ax.plot(alcatraz['x_osm'], alcatraz['y_osm'], 'ro')


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0xde5ebd0>]






    

In [13]:

    

reload(python_subgrid.particles)
system = python_subgrid.particles.ParticleSystem(ds)


    

In [14]:

    

n_new = 500
x_src = np.random.uniform(alcatraz['x_utm'], alcatraz['x_utm']+100, size=n_new)
y_src = np.random.uniform(alcatraz['y_utm'], alcatraz['y_utm']+100, size=n_new)
pts = np.c_[x_src, y_src, np.zeros_like(x_src)]
system.reseed(pts=pts)
#system.tracer.terminal_speed = 1e-2


    

    




INFO:python_subgrid.particles:No output lines found
INFO:python_subgrid.particles:got (0,) particles at t_stop
INFO:python_subgrid.particles:got (0,) particles still alive
INFO:python_subgrid.particles:(0, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (500, 3) points float64

In [15]:

    

system.update_grid(i=270)
system.tracer.update()
system.save_grid()


    

In [15]:

    

x_src = np.random.uniform(ll_utm[0], ur_utm[0], size=n_new)
y_src = np.random.uniform(ll_utm[1], ur_utm[1], size=n_new)
pts = np.c_[x_src, y_src, np.zeros_like(x_src)]


    

In [16]:

    

dfs = []
t0 = 304
t1 = 306
for i in range(t0, 321):
    print times[i]
    system.update_grid(i=i)
    system.tracer.update()
    df = system.get_particles()
    system.reseed(pts)
    dfs.append(df)


    

    




INFO:python_subgrid.particles:got (500,) particles at t_stop
INFO:python_subgrid.particles:got (500,) particles still alive
INFO:python_subgrid.particles:(500, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (1000, 3) points float64






    




1962-06-12 04:00:00
1962-06-12 04:15:00





    




INFO:python_subgrid.particles:got (823,) particles at t_stop
INFO:python_subgrid.particles:got (708,) particles still alive
INFO:python_subgrid.particles:(708, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (1208, 3) points float64






    




1962-06-12 04:30:00





    




INFO:python_subgrid.particles:got (1016,) particles at t_stop
INFO:python_subgrid.particles:got (895,) particles still alive
INFO:python_subgrid.particles:(895, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (1395, 3) points float64






    




1962-06-12 04:45:00





    




INFO:python_subgrid.particles:got (1199,) particles at t_stop
INFO:python_subgrid.particles:got (1084,) particles still alive
INFO:python_subgrid.particles:(1084, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (1584, 3) points float64






    




1962-06-12 05:00:00





    




INFO:python_subgrid.particles:got (1382,) particles at t_stop
INFO:python_subgrid.particles:got (1235,) particles still alive
INFO:python_subgrid.particles:(1235, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (1735, 3) points float64






    




1962-06-12 05:15:00





    




INFO:python_subgrid.particles:got (1541,) particles at t_stop
INFO:python_subgrid.particles:got (1407,) particles still alive
INFO:python_subgrid.particles:(1407, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (1907, 3) points float64






    




1962-06-12 05:30:00





    




INFO:python_subgrid.particles:got (1716,) particles at t_stop
INFO:python_subgrid.particles:got (1560,) particles still alive
INFO:python_subgrid.particles:(1560, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2060, 3) points float64






    




1962-06-12 05:45:00





    




INFO:python_subgrid.particles:got (1866,) particles at t_stop
INFO:python_subgrid.particles:got (1713,) particles still alive
INFO:python_subgrid.particles:(1713, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2213, 3) points float64






    




1962-06-12 06:00:00





    




INFO:python_subgrid.particles:got (2017,) particles at t_stop
INFO:python_subgrid.particles:got (1871,) particles still alive
INFO:python_subgrid.particles:(1871, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2371, 3) points float64






    




1962-06-12 06:15:00





    




INFO:python_subgrid.particles:got (2165,) particles at t_stop
INFO:python_subgrid.particles:got (2016,) particles still alive
INFO:python_subgrid.particles:(2016, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2516, 3) points float64






    




1962-06-12 06:30:00





    




INFO:python_subgrid.particles:got (2317,) particles at t_stop
INFO:python_subgrid.particles:got (2113,) particles still alive
INFO:python_subgrid.particles:(2113, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2613, 3) points float64






    




1962-06-12 06:45:00





    




INFO:python_subgrid.particles:got (2412,) particles at t_stop
INFO:python_subgrid.particles:got (2219,) particles still alive
INFO:python_subgrid.particles:(2219, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2719, 3) points float64






    




1962-06-12 07:00:00





    




INFO:python_subgrid.particles:got (2516,) particles at t_stop
INFO:python_subgrid.particles:got (2334,) particles still alive
INFO:python_subgrid.particles:(2334, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2834, 3) points float64






    




1962-06-12 07:15:00





    




INFO:python_subgrid.particles:got (2619,) particles at t_stop
INFO:python_subgrid.particles:got (2462,) particles still alive
INFO:python_subgrid.particles:(2462, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (2962, 3) points float64






    




1962-06-12 07:30:00





    




INFO:python_subgrid.particles:got (2741,) particles at t_stop
INFO:python_subgrid.particles:got (2536,) particles still alive
INFO:python_subgrid.particles:(2536, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (3036, 3) points float64






    




1962-06-12 07:45:00





    




INFO:python_subgrid.particles:got (2820,) particles at t_stop
INFO:python_subgrid.particles:got (2588,) particles still alive
INFO:python_subgrid.particles:(2588, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (3088, 3) points float64






    




1962-06-12 08:00:00





    




INFO:python_subgrid.particles:got (2881,) particles at t_stop
INFO:python_subgrid.particles:got (2734,) particles still alive
INFO:python_subgrid.particles:(2734, 3) particles used from previous run (float64)
INFO:python_subgrid.particles:Setting particles to (3234, 3) points float64






    

In [19]:

    

df = pandas.concat(dfs)
df.drop_duplicates(cols=['particle', 't'], inplace=True)
df.sort(columns=['particle', 't'], inplace=True)


    

In [21]:

    

fig, ax = plt.subplots(figsize=(20,15))
ax.imshow(img, extent=(ll_osm[0], ur_osm[0], ll_osm[1], ur_osm[1]))


linelist_osm = []
for j, group in df.groupby('particle'):
    x, y, _ = np.array(utm2osm.TransformPoints(np.c_[group['x'], group['y']])).T
    coords = np.c_[x,y]
    linelist_osm.append(coords)
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)


ax.set_xlim(ll_osm[0]+2000, ur_osm[0]-3000)
ax.set_ylim(ll_osm[1]+2000, ur_osm[1]-1000)


    

    
Out[21]:





(4546639.953723437, 4563207.832964446)






    

In [ ]: