Today's Question:  What are you most afraid of as a programmer?

# Build route graph of Hurricane Sandy

sonic0002      2012-11-17 07:58:13      12,189    0    0

Hurricane Sandy swept US east side and landed in New York, it killed 113 persons and incurred 50 billion US dollars economic loss. Sandy is also considered as the most expensive hurricane. We will now use matplotlib and basemap libraries in Python to build a route graph of Sandy.

Below is the animated GIF.

Shadow is added in the graph to show the time at night, we can see from the graph that Sandy stayed a while in Panama after emerging, then it went through Cuba and became Hurricane-2. Later Sandy moved along US east coast and then turned to west suddenly, headed directly to New York and landed in south New York on 29th Oct.

The unit of wind in data file is knot and unit of pressure is mb and time is UTC time.

Source code of drawing the route graph:

```# Written by Vamei

from datetime import datetime, timedelta
import re
import numpy as np

import matplotlib
matplotlib.rcParams['text.color'] = 'white'
matplotlib.rcParams['xtick.color'] = 'white'
matplotlib.rcParams['ytick.color'] = 'white'

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

font = 'monospace'

# This function is to plot the base map
def plotBase(fig, dt=None):
m = Basemap(projection='merc',
lon_0=0,lat_0=0,lat_ts=0,
llcrnrlat=0,urcrnrlat=50,
llcrnrlon=-100,urcrnrlon=-50,
resolution='l')
m.drawcountries(linewidth=1, color='k')
m.drawmapscale(-90, 5, -90, 5, 1000, barstyle='fancy')
m.bluemarble(scale=1)

# Get Position of NYC, longitude -74.0064, latitude 40.7142
x,y    = m(-74.0064, 40.7142)
# Plot NYC
m.scatter(x, y, s=100,  marker='*', color='0.5', alpha=1)
plt.text(x,y,'NYC', fontsize='15')

if dt is not None: m.nightshade(dt, alpha = 0.3)
return m

# Hurricane category colors
color_dict = {'TROPICAL DEPRESSION':'#AEF100', 'TROPICAL STORM':'#FFD600', 'HURRICANE-1':'#FF6440', 'HURRICANE-2':'#8506A9'}

# Read data file, unzip from track.zip to get track.dat
fn  = 'track.dat'
rec = {'lat':[],'lon':[],'wind':[],'press':[],'dt':[],'cat':[]}
for i,line in enumerate(file(fn)):
if i == 0: continue  # Jump over the first line
# replace multiple whitespaces with a single whitespace
line   = re.sub(r"\s+", ' ', line)
pieces = line.split(" ")
# retrieve information
rec['lat'].append(float(pieces[0]))
rec['lon'].append(float(pieces[1]))
rec['wind'].append(float(pieces[3]))
rec['press'].append(float(pieces[4]))
rec['cat'].append((" ".join(pieces[5:])).strip())
time   = pieces[2]
time   = "2012/" + time
rec['dt'].append(datetime.strptime(time,"%Y/%m/%d/%HZ"))

# Plot the track and the else
N = len(rec['lat'])
for idx in range(N):
dt     = rec['dt'][idx]
# Adjust time zone according to NYC
lt     = dt - timedelta(hours=5)
lon    = rec['lon'][idx]
lat    = rec['lat'][idx]
wind   = rec['wind'][idx]
press  = rec['press'][idx]
cat    = rec['cat'][idx]
fig    = plt.figure()
m      = plotBase(fig, dt)
# From lon,lat to pixels
x,y    = m(lon, lat)
# Plot track
for i in range(idx):
a0,b0 = m(rec['lon'][i], rec['lat'][i])
a1,b1 = m(rec['lon'][i+1], rec['lat'][i+1])
m.plot((a0,a1),(b0,b1), linewidth=2.5,
color=color_dict[rec['cat'][i+1]])
# Plot Sandy's current position
m.scatter(x, y, s=100, c=color_dict[cat], alpha=0.8)
# Annotate current position
plt.annotate(
cat,
xy = (x, y), xytext = (-5, -30),
textcoords = 'offset points', ha = 'right', va = 'bottom',
bbox = dict(boxstyle = 'round,pad=0.5', fc=color_dict[cat], alpha = 0.8),
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))

tx, ty = m(-98, 40)
plt.text(tx,ty,
lt.strftime("Hurricane Sandy\n\nNYC LT:\n%Y-%m-%d %H:00:00\nData Source: NOAA\nBy Vamei"),
family=font,ha='left')
# add a small axes to show pressure
a.set_ylim((950,1000))
a.set_xlim((-10,70))
a.set_yticks([900,1050])
a.set_xticks([0,60])
a.set_title("Center Pressure (mb)", fontsize=10)
a.plot(rec['press'])
a.axvline(x=idx, color='r')
# add a small axes to show wind
a.set_ylim((0,100))
a.set_xlim((-10,70))
a.set_yticks([0,100])
a.set_xticks([0,60])
a.set_title("Max Wind (knots)", fontsize=10)
a.plot(rec['wind'])
a.axvline(x=idx, color='r')

fig.savefig(('%04d.png' % idx))
plt.close()```

Finally, hope those people who lost their lives in this natural disaster rest in peace.

Original author : Vamei Source : http://www.cnblogs.com/vamei/archive/2012/11/07/2758006.html