Roberto De Almeida roberto@marinexplore.com
This is notebook inspired by the article Straight as an arrow: humpback whales swim constant course tracks during long-distance migration. Here we adapt Dijkstra's algorithm to calculate the shortest (in time) path between two points given a non-homogeneous, non-stationary field of superficial ocean currents. A similar application of the algorithm can be seen in the article Minimum-Risk Planning for AUVs using Ocean Current Predictions.
import heapq
import numpy as np
import scipy.interpolate
from mpl_toolkits.basemap import Basemap
from IPython.display import HTML
from pydap.client import open_url
import coards
Whale tracks are shown on the Fig. 1 from Horton et al. 2011. Since satellite derived ocean currents used in this notebook start in October 2009, we compare it to the circles in Fig 1b., which represent the whales also tracked in 2009.
HTML("""
<a href="http://rsbl.royalsocietypublishing.org/content/7/5/674/F1.large.jpg" target="_blank"><img src="http://rsbl.royalsocietypublishing.org/content/7/5/674/F1.medium.gif"></a>
""")
This dataset at Marinexplore contains Ekman currents from 2009-10-03 to 2010-02-12. It can be accessed directly using OPeNDAP:
dataset1 = open_url("http://podaac-opendap.jpl.nasa.gov/opendap/allData/oscar/preview/L4/oscar_third_deg/oscar_vel2009.nc.gz")
dataset2 = open_url("http://podaac-opendap.jpl.nasa.gov/opendap/allData/oscar/preview/L4/oscar_third_deg/oscar_vel2010.nc.gz")
print dataset1
print dataset2
DatasetType with children ['time', 'year', 'depth', 'latitude', 'longitude', 'u', 'v', 'um', 'vm'] DatasetType with children ['time', 'year', 'depth', 'latitude', 'longitude', 'u', 'v', 'um', 'vm']
i0, i1 = ((310 <= dataset1.longitude) & (dataset1.longitude <= 355)).nonzero()[0][[0,-1]]
j0, j1 = ((-65 <= dataset1.latitude) & (dataset1.latitude <= 0)).nonzero()[0][[0,-1]]
U = np.concatenate((dataset1.u.u[:,0,j0:j1,i0:i1], dataset2.u.u[:,0,j0:j1,i0:i1]), axis=0)
U = np.ma.masked_where(isnan(U), U)
U = U[:,0] # remove vertical dimension
V = np.concatenate((dataset1.v.v[:,0,j0:j1,i0:i1], dataset2.v.v[:,0,j0:j1,i0:i1]), axis=0)
V = np.ma.masked_where(isnan(V), V)
V = V[:,0]
x = dataset1.longitude[i0:i1]
y = dataset1.latitude[j0:j1]
X, Y = np.meshgrid(x, y)
fig = figure(figsize=(15,9))
ax = fig.gca()
ax.set_title('OSCAR currents (%s)' % coards.parse(dataset1.time[0][0], dataset1.time.units))
m = Basemap(ax=ax, llcrnrlon=310, llcrnrlat=-65, urcrnrlon=355, urcrnrlat=0)
m.drawcoastlines(linewidth=0.25)
m.fillcontinents(color='#000000', lake_color='#666666')
m.drawmapboundary(fill_color='#666666')
colors = m.contourf(X, Y, np.sqrt(U[0]**2 + V[0]**2))
m.quiver(X[::5,::5], Y[::5,::5], U[0,::5,::5], V[0,::5,::5], color='white')
cb = pylab.colorbar(colors)
cb.set_label('m/s')
Let's increase the resolution, interpolating linearly:
xi = linspace(310, 355, 500)
yi = linspace(0, -65, 500)
Xi, Yi = np.meshgrid(xi, yi)
Ui = np.ma.masked_all((U.shape[0], len(yi), len(xi)), U.dtype)
Vi = np.ma.masked_all((V.shape[0], len(yi), len(xi)), V.dtype)
for l in range(U.shape[0]):
valid = ~U[l].mask & ~V[l].mask
points = X[valid], Y[valid]
Ui[l] = scipy.interpolate.griddata(points, U[l][valid], (Xi, Yi)).reshape(Xi.shape)
Vi[l] = scipy.interpolate.griddata(points, V[l][valid], (Xi, Yi)).reshape(Xi.shape)
We can replace missing values with zero, and allow whales to navigate there on their own:
Ui[isnan(Ui)] = 0
Vi[isnan(Vi)] = 0
fig = figure(figsize=(15,9))
ax = fig.gca()
ax.set_title('OSCAR currents (%s)' % coards.parse(dataset1.time[0][0], dataset1.time.units))
m = Basemap(ax=ax, llcrnrlon=310, llcrnrlat=-65, urcrnrlon=355, urcrnrlat=0)
m.drawcoastlines(linewidth=0.25)
m.fillcontinents(color='#000000', lake_color='#666666')
m.drawmapboundary(fill_color='#666666')
colors = m.contourf(Xi, Yi, np.sqrt(Ui[0]**2 + Vi[0]**2))
m.quiver(Xi[::10,::10], Yi[::10,::10], Ui[0,::10,::10], Vi[0,::10,::10], color='white')
cb = pylab.colorbar(colors)
cb.set_label('m/s')
Now that we have the current field we can compute the best trajectory for a whale travelling at 2.5m/s. Let's try to simulate the eastern most trajectory from Fig 1b, starting from 12S 37W and going to 60S 25W:
#2009-10-03 to 2010-02-12
fig = figure(figsize=(15,9))
ax = fig.gca()
ax.set_title('OSCAR currents (%s)' % coards.parse(dataset1.time[0][0], dataset1.time.units))
#m = Basemap(ax=ax, llcrnrlon=-50, llcrnrlat=-65, urcrnrlon=-5, urcrnrlat=0)
width = 10000000
height = 10000000
m = Basemap(width=width, height=height, projection='aeqd', lat_0=-36, lon_0=323)
m.drawcoastlines(linewidth=0.25)
m.fillcontinents(color='#000000', lake_color='#666666')
m.drawmapboundary(fill_color='#666666')
xx, yy = m(Xi[::10,::10], Yi[::10,::10])
m.quiver(xx, yy, Ui[0,::10,::10], Vi[0,::10,::10], color='white')
xx, yy = m([-37, -25], [-12, -60])
m.plot(xx, yy, 'r-', linewidth=2)
[<matplotlib.lines.Line2D at 0xc006d10>]
def shortest_path1(X, Y, U, V, T, G, start, end, s0=0, t0=0):
q = [(t0, start, ())]
visited = {}
while 1:
cost, v1, path = heapq.heappop(q)
if v1 not in visited or visited[v1] > cost:
path = path + (v1,)
visited[v1] = cost
if v1 == end:
return list(path), cost
for v2 in G[v1]:
cost2 = calculate_cost(X, Y, U, V, T, cost, v1, v2, s0)
if (v2 not in visited or visited[v2] > cost2) and cost+cost2 <= T[-1]:
heapq.heappush(q, (cost + cost2, v2, path))
def shortest_path2(X, Y, U, V, T, G, start, end, s0=0, t0=0):
q = [(t0, start, ())]
visited = set()
while 1:
cost, v1, path = heapq.heappop(q)
if v1 not in visited:
path = path + (v1,)
visited.add(v1)
if v1 == end:
return list(path), cost
for v2 in G[v1]:
cost2 = calculate_cost(X, Y, U, V, T, cost, v1, v2, s0)
if v2 not in visited and cost+cost2 <= T[-1]:
heapq.heappush(q, (cost + cost2, v2, path))
shortest_path = shortest_path1
def distance(lat1, lon1, lat2, lon2):
# http://www.movable-type.co.uk/scripts/latlong.html
R = 6.371e6
lat1 *= pi/180.
lon1 *= pi/180.
lat2 *= pi/180.
lon2 *= pi/180.
return R*np.arccos(
np.sin(lat1)*np.sin(lat2) +
np.cos(lat1)*np.cos(lat2)*np.cos(lon2-lon1))
def bearing(lat1, lon1, lat2, lon2):
# http://www.movable-type.co.uk/scripts/latlong.html
lat1 *= pi/180.
lon1 *= pi/180.
lat2 *= pi/180.
lon2 *= pi/180.
y = np.sin(lon2-lon1)*np.cos(lat2)
x = np.cos(lat1)*np.sin(lat2) - np.sin(lat1)*np.cos(lat2)*np.cos(lon2-lon1)
return np.arctan2(y, x)
def calculate_cost(X, Y, U, V, T, t, v1, v2, s0):
l = (T <= t).nonzero()[0].max()
j1, i1 = v1
j2, i2 = v2
u = (U[l,j1,i1] + U[l,j2,i2])/2.
v = (V[l,j1,i1] + V[l,j2,i2])/2.
ds = distance(Y[v1], X[v1], Y[v2], X[v2])
a = bearing(Y[v1], X[v1], Y[v2], X[v2])
# velocity along track
s = s0 + u*np.cos(a) + v*np.sin(a)
if s < 0:
return np.inf
else:
return ds/s
G = {}
for i in range(len(xi)):
for j in range(len(yi)):
G[j, i] = []
if i > 0:
if j > 0:
G[j, i].append((j-1, i-1))
if j+1 < len(yi):
G[j, i].append((j+1, i-1))
G[j, i].append((j, i-1))
if i+1 < len(xi):
if j > 0:
G[j, i].append((j-1, i+1))
if j+1 < len(yi):
G[j, i].append((j+1, i+1))
G[j, i].append((j, i+1))
if j > 0:
G[j, i].append((j-1, i))
if j+1 < len(yi):
G[j, i].append((j+1, i))
s0 = 1.0 # m/s
# find start 12S 37W
i = (xi <= 323).nonzero()[0].max()
j = (yi <= -12).nonzero()[0].min()
#j = (yi <= -30).nonzero()[0].min()
start = j, i
# find end 60S 25W
i = (xi <= 335).nonzero()[0].max()
j = (yi <= -60).nonzero()[0].min()
end = j, i
T = np.concatenate((dataset1.time[:], dataset2.time[:]), axis=0) * 24 * 60 * 60 # convert from days to seconds
path, end = shortest_path(Xi, Yi, Ui, Vi, T, G, start, end, s0, T[51])
fig = figure(figsize=(15,9))
ax = fig.gca()
ax.set_title('OSCAR currents (%s)' % coards.parse(dataset1.time[0][0], dataset1.time.units))
m = Basemap(ax=ax, llcrnrlon=310, llcrnrlat=-65, urcrnrlon=355, urcrnrlat=0)
width = 8000000
height = 8000000
#m = Basemap(width=width, height=height, projection='aeqd', lat_0=-36, lon_0=-37)
m.drawcoastlines(linewidth=0.25)
m.fillcontinents(color='#000000', lake_color='#666666')
m.drawmapboundary(fill_color='#666666')
xx, yy = m(Xi[::10,::10], Yi[::10,::10])
m.quiver(xx, yy, Ui[0,::10,::10], Vi[0,::10,::10], color='white')
xx = [Xi[s] for s in path]
yy = [Yi[s] for s in path]
xx, yy = m(xx, yy)
m.plot(xx, yy, 'ro')
xx, yy = m([323, 335], [-12, -60])
m.plot(xx, yy, 'r-', linewidth=2)
fig.savefig('traj.png')
print coards.parse(T[51]/(24*60*60), dataset1.time.units)
print coards.parse(end/(24*60*60), dataset1.time.units)
2009-09-16 00:00:00 2009-11-18 17:33:24.974011
fig = figure(figsize=(15,9))
ax = fig.gca()
ax.set_title('OSCAR currents (%s)' % coards.parse(dataset1.time[0][0], dataset1.time.units))
#m = Basemap(ax=ax, llcrnrlon=310, llcrnrlat=-65, urcrnrlon=355, urcrnrlat=0)
width = 8000000
height = 8000000
m = Basemap(width=width, height=height, projection='aeqd', lat_0=-36, lon_0=-37)
m.drawcoastlines(linewidth=0.25)
m.fillcontinents(color='#000000', lake_color='#666666')
m.drawmapboundary(fill_color='#666666')
xx, yy = m(Xi[::10,::10], Yi[::10,::10])
m.quiver(xx, yy, Ui[0,::10,::10], Vi[0,::10,::10], color='white')
xx = [Xi[s] for s in path]
yy = [Yi[s] for s in path]
xx, yy = m(xx, yy)
m.plot(xx, yy, 'ro')
xx, yy = m([323, 335], [-12, -60])
m.plot(xx, yy, 'r-', linewidth=2)
fig.savefig('traj.png')