In [1]:

    

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
logging.basicConfig()
logging.root.setLevel(logging.DEBUG)


    

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

    

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

    

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





((-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 [5]:

    

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

    

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





<matplotlib.image.AxesImage at 0x8fc3a90>






    

In [7]:

    

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


    

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