def draw_client_density():

    m = Basemap(llcrnrlon=-180.,llcrnrlat=-90,urcrnrlon=180.,urcrnrlat=90.,\
                resolution='c',projection='cyl')

    # plot them as filled circles on the map.
    # first, create a figure.
    dpi=100
    dimx=800/dpi
    dimy=400/dpi
    fig=figure(figsize=(dimx,dimy), dpi=dpi, frameon=False, facecolor='blue')
#    ax=fig.add_axes([0.1,0.1,0.7,0.7],axisbg='g')
    ax=fig.add_axes([0.0,0.0,1.0,1.0],axisbg='g')
    canvas = FigureCanvas(fig)
    results = lookup_client_locations()
    X,Y,Z = find_client_density(m,results)
#    s = random.sample(results, 40000)
#    for t in s:
#        lat=t[2]
#        lon=t[3]
#        # draw a red dot at the center.
#        xpt, ypt = m(lon, lat)
#        m.plot([xpt],[ypt],'ro', zorder=10)
    # draw coasts and fill continents.
    m.drawcoastlines(linewidth=0.5)
    m.drawcountries(linewidth=0.5)
    m.drawlsmask([100,100,100,0],[0,0,255,255])
#    m.fillcontinents(color='green')
    palette = cm.YlOrRd
    m.imshow(Z,palette,extent=(m.xmin,m.xmax,m.ymin,m.ymax),interpolation='gaussian',zorder=0)
#    l,b,w,h = ax.get_position()
#    cax = axes([l+w+0.075, b, 0.05, h])
#    colorbar(cax=cax) # draw colorbar

    canvas.print_figure(outdir+'/clientmap.png', dpi=100)

def InvokeMap(coastfile='/media/sda4/map-data/aust-coast-noaa-2000000-1.dat',
		    lllon=80,
		    urlon=166,
		    lllat=-47,
		    urlat=-9,
		    draw_map=True):
    global PYLIB_PATH

    map = Basemap(projection='cyl',
			llcrnrlon=lllon,
			urcrnrlon=urlon,
			llcrnrlat=lllat,
			urcrnrlat=urlat,
			#lat_ts=-35,
			lat_0=-35,
			lon_0=120,
			resolution='l',
			area_thresh=1000.)


    try: 
	coast = p.load(coastfile)
	coast = p.load(coastfile)
	coast_x,coast_y = map(coast[:,0],coast[:,1])
	p.plot(coast_x,coast_y,color='black')    
    except IOError:
	map.drawcoastlines()

    map.drawmapboundary()
    map.drawmeridians(p.arange(0,360,10),labels=[0,0,1,0])
    map.drawparallels(p.arange(-90,0,10),labels=[1,0,0,0])

    return map

def doit():
    map = Basemap(projection='lcc',
		    llcrnrlon=80,
		    urcrnrlon=160,
		    llcrnrlat=-50,
		    urcrnrlat=-8,
		    #lat_ts=-35,
		    lat_0=-35,
		    lon_0=120,
		    resolution='c',
		    area_thresh=1000.)
    p.clf()
    map.drawcoastlines()
    # map.drawcountries()
    
    # map.drawrivers()

    map.drawmeridians(p.arange(0,360,10),labels=[0,0,1,0])
    map.drawparallels(p.arange(-90,0,10),labels=[1,0,0,0])

    traj=p.load('example_traj.dat')
    coast=p.load('/media/sda4/map-data/aust-coast-noaa-2000000-1.dat')

    traj_x,traj_y   = map(traj[:,1],traj[:,0]) 
    # coast_x,coast_y = map(coast[:,0],coast[:,1])
    
    p.plot(traj_x,traj_y)    
    p.plot(coast_x,coast_y,color='black')    

    map.drawmapboundary()
    p.show()
    return map 

# read in jpeg image to rgba array of normalized floats.
pilImage = Image.open('land_shallow_topo_2048.jpg')
rgba = pil_to_array(pilImage)
rgba = rgba.astype(P.Float32)/255. # convert to normalized floats.

# define lat/lon grid that image spans (projection='cyl').
nlons = rgba.shape[1]; nlats = rgba.shape[0]
delta = 360./float(nlons)
lons = P.arange(-180.+0.5*delta,180.,delta)
lats = P.arange(-90.+0.5*delta,90.,delta)

# create new figure
fig=P.figure()
# define cylindrical equidistant projection.
m = Basemap(projection='cyl',llcrnrlon=-180,llcrnrlat=-90,urcrnrlon=180,urcrnrlat=90,resolution='l')
# plot (unwarped) rgba image.
im = m.imshow(rgba)
# draw coastlines.
m.drawcoastlines(linewidth=0.5,color='0.5')
# draw lat/lon grid lines.
m.drawmeridians(P.arange(-180,180,60),labels=[0,0,0,1],color='0.5')
m.drawparallels(P.arange(-90,90,30),labels=[1,0,0,0],color='0.5')
P.title("Blue Marble image - native 'cyl' projection",fontsize=12)
print 'plot cylindrical map (no warping needed) ...'

# create new figure
fig=P.figure()
# define orthographic projection centered on North America.
m = Basemap(projection='ortho',lat_0=40,lon_0=40,resolution='l')
# transform to nx x ny regularly spaced native projection grid

    x = property(get_xdata)
    def get_ydata(self): return self._line.get_ydata()
    y = property(get_ydata)
    


if __name__ == '__main__':
    from matplotlib.toolkits.basemap import Basemap
    from matplotlib.ticker import FuncFormatter
    from pyroms.focus import Focus
    
    m = Basemap(projection='lcc', 
                resolution='i',
                llcrnrlon=5.0,
                llcrnrlat= 52.,
                urcrnrlon=35.5,
                urcrnrlat=68.0,
                lat_0=15.00, 
                lon_0=0.0,
                suppress_ticks=False)
    
    fig=pl.figure()
    # background color will be used for 'wet' areas.
    fig.add_axes([0.1,0.1,0.8,0.8],axisbg='azure')
    m.drawcoastlines(linewidth=0.25)
    m.fillcontinents(color='beige')
    
    x, y, beta = zip(*[(241782.65384467551, 1019981.3539730886,   1.0),
                       (263546.12512432877,  686274.79435173667,  0),
                       (452162.87621465814,  207478.42619936226,  1.0),
                       (1431519.0837990609,  432367.62942244718,  1.0),

lats = 90.-delat*arange(nlats)
lons, lats = meshgrid(lons, lats)
lons = (d2r*lons.flat).tolist()
lats = (d2r*lats.flat).tolist()
# randomly shuffle locations.
random.shuffle(lons)
random.shuffle(lats)
lons = array(lons,'f')
lats = array(lats,'f')

# minimum separation distance in km.
rcrit = 500.

# set up lambert azimuthal map centered on N. Pole.
m = Basemap(llcrnrlon=-150.,llcrnrlat=0.,urcrnrlon=30.,urcrnrlat=0.,
            resolution='l',area_thresh=10000.,projection='stere',\
            lat_0=90.,lon_0=-105.,lat_ts=90.)

print len(lons), ' obs before thinning'

# calculate distance between each ob and all preceding obs in list.
# throw out those that are closer than rcrit.
nob = 0
lats_out = []
lons_out = []
for lon,lat in zip(lons,lats):
   if nob:
      r = (m.rmajor/1000.)*get_dist(lon,lons[0:nob],lat,lats[0:nob])
      if min(r) > rcrit: 
          lats_out.append(lat)
          lons_out.append(lon)

                  figure, title, meshgrid, cm, arange

# examples of filled contour plots on map projections.

# read in data on lat/lon grid.
hgt = load('500hgtdata.gz')
lons = load('500hgtlons.gz')
lats = load('500hgtlats.gz')
# shift data so lons go from -180 to 180 instead of 0 to 360.
hgt,lons = shiftgrid(180.,hgt,lons,start=False)
lons, lats = meshgrid(lons, lats)

# create new figure
fig=figure()
# setup of sinusoidal basemap
m = Basemap(resolution='c',projection='sinu',lon_0=0)
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# make a filled contour plot.
x, y = m(lons, lats)
CS = m.contour(x,y,hgt,15,linewidths=0.5,colors='k')
CS = m.contourf(x,y,hgt,15,cmap=cm.jet)
l,b,w,h=ax.get_position()
cax = axes([l+w+0.075, b, 0.05, h]) # setup colorbar axes
colorbar(drawedges=True, cax=cax) # draw colorbar
axes(ax)  # make the original axes current again
# draw coastlines and political boundaries.
m.drawcoastlines()
m.drawmapboundary()
m.fillcontinents()
# draw parallels and meridians.
parallels = arange(-60.,90,30.)

   """

   def __init__(self, baseMap):
      """baseMap is the Basemap object that will be used to translate between
      lon/lat and native map coordinates.
      """
      self.baseMap = baseMap

   def __call__(self, y, pos=1):
      """Return the label for tick value y at position pos.
      """
      lon, lat = self.baseMap(0.0, y, inverse=True)
      return "%g" % lat

m = Basemap(-180.,-80.,180.,80.,\
            resolution='c',area_thresh=10000.,projection='merc',\
            lat_ts=20.)
xsize = rcParams['figure.figsize'][0]
fig=figure(figsize=(xsize,m.aspect*xsize))
fig.add_axes([0.1,0.1,0.8,0.8])
ax = gca() # get current axis instance
ax.update_datalim(((m.llcrnrx, m.llcrnry),(m.urcrnrx,m.urcrnry)))
ax.set_xlim((m.llcrnrx, m.urcrnrx))
ax.set_ylim((m.llcrnry, m.urcrnry))
m.drawcoastlines(ax)
m.drawcountries(ax)
m.drawstates(ax)
m.fillcontinents(ax)
# draw parallels
delat = 30.
circles = arange(0.,90.,delat).tolist()+\

from matplotlib.toolkits.basemap import Basemap, interp
from pylab import *
import cPickle

# read in data on lat/lon grid.
datadict = cPickle.load(open('500hgt.pickle','rb'))
hgt = datadict['data']; lons = datadict['lons']; lats = datadict['lats']

# set up map projection (lambert azimuthal equal area).
m = Basemap(-150.,-20.,30.,-20.,
             resolution='c',area_thresh=10000.,projection='laea',
             lat_0=90.,lon_0=-105.)
# interpolate to map projection grid.
nx = 101
ny = 101
lonsout, latsout = m.makegrid(nx,ny)
# get rid of negative lons.
lonsout = where(lonsout < 0., lonsout + 360., lonsout)
hgt = interp(hgt,lons,lats,lonsout,latsout)
dx = (m.xmax-m.xmin)/(nx-1)
dy = (m.ymax-m.ymin)/(ny-1)  
x = m.xmin+dx*indices((ny,nx))[1,:,:]
y = m.ymin+dy*indices((ny,nx))[0,:,:]

#m = Basemap(lons[0],lats[0],lons[-1],lats[-1],\
#              resolution='c',area_thresh=10000.,projection='cyl')
#x, y = meshgrid(lons, lats)

fig = figure(figsize=(8,8))

plots = ['contour','pcolor']

from matplotlib.toolkits.basemap import Basemap
from pylab import *
import cPickle

# create Basemap instance. Use 'crude' resolution coastlines.
m = Basemap(llcrnrlon=-11.,llcrnrlat=50.5,urcrnrlon=-5.,urcrnrlat=56.,
            resolution='c',area_thresh=1000.,projection='tmerc',lon_0=-8.,lat_0=0.)
# create figure with same aspect ratio as map.
xsize = rcParams['figure.figsize'][0]
fig=figure(figsize=(xsize,m.aspect*xsize))
fig.add_axes([0.1,0.1,0.8,0.8])
# draw coastlines and fill continents.
m.drawcoastlines()
m.fillcontinents()
# draw political boundaries.
m.drawcountries()
# draw parallels
circles = arange(50,60,1).tolist()
m.drawparallels(circles,labels=[1,1,0,0])
# draw meridians
meridians = arange(-12,0,1)
m.drawmeridians(meridians,labels=[0,0,1,1])
title("Crude Res Boundaries ('c')",y=1.075)
show()

# create Basemap instance. Use 'low' resolution coastlines.
m = Basemap(llcrnrlon=-11.,llcrnrlat=50.5,urcrnrlon=-5.,urcrnrlat=56.,
            resolution='l',area_thresh=1000.,projection='tmerc',lon_0=-8.,lat_0=0.)
# create figure with same aspect ratio as map.
xsize = rcParams['figure.figsize'][0]
fig=figure(figsize=(xsize,m.aspect*xsize))

from matplotlib.toolkits.basemap import Basemap, shiftgrid
from pylab import *

# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
topoin = array(load('etopo20data.gz'),'d')
lons = array(load('etopo20lons.gz'),'d')
lats = array(load('etopo20lats.gz'),'d')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topoin,lons,start=False)

# setup of basemap ('lcc' = lambert conformal conic).
# use major and minor sphere radii from WGS84 ellipsoid.
m = Basemap(llcrnrlon=-145.5,llcrnrlat=1.,urcrnrlon=-2.566,urcrnrlat=46.352,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',area_thresh=1000.,projection='lcc',\
            lat_1=50.,lon_0=-107.)
# transform to nx x ny regularly spaced native projection grid
nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
topodat,x,y = m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True)
# create the figure.
fig=figure(figsize=(8,8))
# add an axes, leaving room for colorbar on the right.
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# plot image over map with imshow.
im = m.imshow(topodat,cm.jet)
# setup colorbar axes instance.
l,b,w,h = ax.get_position()
cax = axes([l+w+0.075, b, 0.05, h])
colorbar(tickfmt='%d', cax=cax) # draw colorbar
axes(ax)  # make the original axes current again

    data.y, 2, oban.barnes_weights, (kappa0, g))
  barnes_analyses.append(field)
  diffs = oban.get_ob_incs(data.x, data.y, data.heights, x, y, field)
  barnes_rms.append(misc.rms(diffs))
  print "Barnes 2-pass gamma=%.1f rms: %f" % (g, barnes_rms[-1])

#Do the barnes 3-pass
barnes_3pass = oban.analyze_grid_multipass(data.heights, x_grid, y_grid,
  data.x, data.y, 3, oban.barnes_weights, (kappa0, 1.0))
diffs = oban.get_ob_incs(data.x, data.y, data.heights, x, y, barnes_3pass)
barnes_3pass_rms = misc.rms(diffs)
print "Barnes 3-pass rms: %f" % barnes_3pass_rms

#Generate a grid for basemap plotting
bm_ps = Basemap(projection='stere',lat_0=90.0,lon_0=-110.0,lat_ts=60.0,
  rsphere=ps.radius,llcrnrlon=-120.5,llcrnrlat=25.1,urcrnrlon=-61.8,
  urcrnrlat=43.4,resolution='l',area_thresh=5000)

#Transform the grid to basemap space for plotting
x_bm,y_bm = ps.to_basemap_xy(bm_ps, x_grid / cm_per_m, y_grid / cm_per_m)

#Check if we want to view or save output
if len(sys.argv) > 1 and sys.argv[1].startswith('silent'):
  save_work = True
  colormap = M.cm.get_cmap('gist_gray')
  fontsize = 10
else:
  save_work = False
  colormap = M.cm.get_cmap('jet')
  fontsize = 14

# make plot of etopo bathymetry/topography data on
# lambert conformal conic map projection, drawing coastlines, state and
# country boundaries, and parallels/meridians.

from matplotlib.toolkits.basemap import Basemap, interp
from pylab import *
import cPickle

# read in topo data from pickle (on a regular lat/lon grid)
# longitudes go from 20 to 380.
topodict = cPickle.load(open('data/etopo20.pickle','rb'))
topoin = topodict['data']; lons = topodict['lons']; lats = topodict['lats']

# setup of basemap ('lcc' = lambert conformal conic).
m = Basemap(-145.5,1.,-2.566,46.352,\
            resolution='c',area_thresh=10000.,projection='lcc',\
            lat_1=50.,lon_0=-107.)
# define grid (nx x ny regularly spaced native projection grid)
nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
lonsout, latsout = m.makegrid(nx,ny)
topodat = interp(topoin,lons,lats,lonsout,latsout)
xsize = rcParams['figure.figsize'][0]
fig=figure(figsize=(xsize,m.aspect*xsize))
ax = fig.add_axes([0.1,0.1,0.75,0.75])
im = imshow(topodat,cm.jet,extent=(m.xmin, m.xmax, m.ymin, m.ymax),origin='lower')

cax = axes([0.875, 0.1, 0.05, 0.75])
colorbar(tickfmt='%d', cax=cax) # draw colorbar

axes(ax)  # make the original axes current again
m.drawcoastlines(ax)

    u1[:,-1] = u1_all[time_idx,:,0]

tmp = copy.copy(lon)
lon = p.zeros((p.size(tmp)+1,))
lon[0:-1] = tmp[:]
lon[-1] = tmp[0]+360
del tmp

f.close()


#--- Mapping information:

map = Basemap( projection='cyl', resolution='l'
             , llcrnrlon=0, urcrnrlon=360
             , llcrnrlat=-76.875, urcrnrlat=76.875
             , lon_0=180, lat_0=0
             )
map.drawmeridians(p.arange(0,361,45), labels=[0,0,0,1])
map.drawparallels(p.arange(-90,90,30), labels=[1,0,0,1])
map.drawcoastlines()


#--- Write out contour map and view using preview:

x, y = p.meshgrid(lon,lat)
CS = map.contourf(x, y, u1, cmap=p.cm.gray)
p.text( 0.5, -0.15, 'Longitude [deg]'
      , horizontalalignment='center'
      , verticalalignment='center'
      , transform = p.gca().transAxes )

	def DrawSiteNetwork(self, g, node_label2pos_counts,pic_area=[-180,-90,180,90], output_fname_prefix=None):
		"""
		2007-07-17
			put ax.plot() right after Basemap() but after m.xxx() so that it'll zoom in
			use 'g' in ax.plot(), otherwise, ax.plot() alternates all colors.
			no parallels, no meridians
		2007-08-29 copied from CreatePopulation.py and renamed from DrawStrainNetwork
		"""
		sys.stderr.write("Drawing Site Network...")
		import pylab
		from matplotlib.toolkits.basemap import Basemap
		pylab.clf()
		fig = pylab.figure()
		fig.add_axes([0.05,0.05,0.9,0.9])	#[left, bottom, width, height]
		m = Basemap(llcrnrlon=pic_area[0],llcrnrlat=pic_area[1],urcrnrlon=pic_area[2],urcrnrlat=pic_area[3],\
		resolution='l',projection='mill')
		
		ax=pylab.gca()
		for e in g.edges():
			lat1, lon1 = node_label2pos_counts[e[0]][0]
			lat2, lon2 = node_label2pos_counts[e[1]][0]
			x1, y1 = m(lon1, lat1)
			x2, y2 = m(lon2, lat2)
			ax.plot([x1,x2],[y1,y2], 'g', alpha=0.5, zorder=12)
		
		#m.drawcoastlines()
		m.drawparallels(pylab.arange(-90,90,30), labels=[1,1,0,1])
		m.drawmeridians(pylab.arange(-180,180,30), labels=[1,1,0,1])
		m.fillcontinents()
		m.drawcountries()
		m.drawstates()

		pylab.title("Network of strains")
		if output_fname_prefix:
			pylab.savefig('%s_site_network.eps'%output_fname_prefix, dpi=300)
			pylab.savefig('%s_site_network.svg'%output_fname_prefix, dpi=300)
			pylab.savefig('%s_site_network.png'%output_fname_prefix, dpi=300)
		del m, pylab, Basemap
		sys.stderr.write("Done.\n")

def draw_graph_on_map(g, node2weight, node2pos, pic_title,  pic_area=[-130,10,140,70], output_fname_prefix=None, need_draw_edge=0):
	"""
	2007-09-13
		identity_pair_ls is a list of pairs of strains (ecotype id as in table ecotype)
	2007-10-08
		correct a bug in 4*diameter_ls, diameter_ls has to be converted to array first.
		sqrt the node weight, 8 times the original weight
	"""
	import os, sys
	sys.stderr.write("Drawing graph on a map ...\n")
	import pylab, math
	from matplotlib.toolkits.basemap import Basemap
	pylab.clf()
	fig = pylab.figure()
	fig.add_axes([0.05,0.05,0.9,0.9])	#[left, bottom, width, height]
	m = Basemap(llcrnrlon=pic_area[0],llcrnrlat=pic_area[1],urcrnrlon=pic_area[2],urcrnrlat=pic_area[3],\
	resolution='l',projection='mill', ax=pylab.gca())
	
	sys.stderr.write("\tDrawing nodes ...")
	euc_coord1_ls = []
	euc_coord2_ls = []
	diameter_ls = []
	for n in g.nodes():
		lat, lon = node2pos[n]
		euc_coord1, euc_coord2 = m(lon, lat)	#longitude first, latitude 2nd
		euc_coord1_ls.append(euc_coord1)
		euc_coord2_ls.append(euc_coord2)
		diameter_ls.append(math.sqrt(node2weight[n]))
	import numpy
	diameter_ls = numpy.array(diameter_ls)
	m.scatter(euc_coord1_ls, euc_coord2_ls, 8*diameter_ls, marker='o', color='r', alpha=0.4, zorder=12, faceted=False)
	sys.stderr.write("Done.\n")
	
	if need_draw_edge:
		sys.stderr.write("\tDrawing edges ...")
		ax=pylab.gca()
		for popid1, popid2, no_of_connections in g.edges():
			lat1, lon1 = node2pos[popid1]
			lat2, lon2 = node2pos[popid2]
			x1, y1 = m(lon1, lat1)
			x2, y2 = m(lon2, lat2)
			ax.plot([x1,x2],[y1,y2], 'g', linewidth=math.log(no_of_connections+1)/2, alpha=0.2, zorder=10)
		sys.stderr.write("Done.\n")
	
	#m.drawcoastlines()
	m.drawparallels(pylab.arange(-90,90,30), labels=[1,1,0,1])
	m.drawmeridians(pylab.arange(-180,180,30), labels=[1,1,0,1])
	m.fillcontinents()
	m.drawcountries()
	m.drawstates()
	
	pylab.title(pic_title)
	if output_fname_prefix:
		pylab.savefig('%s.eps'%output_fname_prefix, dpi=600)
		pylab.savefig('%s.svg'%output_fname_prefix, dpi=600)
		pylab.savefig('%s.png'%output_fname_prefix, dpi=600)
	del fig, m, pylab
	sys.stderr.write("Done.\n")

pylab.clf()

b = 5.
# Plot the region surrounding the coastline rounded to the nearest 5 degree lat/lon.
lonmin = b*math.floor(min(lon)/b)
lonmax = b*math.ceil(max(lon)/b)
latmin = b*math.floor(min(lat)/b)
latmax = b*math.ceil(max(lat)/b)
meridian = (range(lonmin, lonmax+b, b))
parallel = (range(latmin, latmax+b, b))
prj = 'cyl'
res = 'i'
map = Basemap(projection=prj,
              llcrnrlon=lonmin,
              llcrnrlat=latmin,
              urcrnrlon=lonmax,
              urcrnrlat=latmax,
              resolution=res)

map.plot(lon, lat, 'r-')
map.plot(lon, lat, 'k+')

fs = 10

pylab.text(115.9,-31.9, "Perth", fontsize=fs, ha='left', va='top')
pylab.text(114.1,-21.9, "Exmouth", fontsize=fs, ha='left', va='top')
pylab.text(118.5,-20.4, "Port Hedland", fontsize=fs, ha='left', va='top')
pylab.text(130.8,-12.6, "Darwin", fontsize=fs, ha='left', va='top')
pylab.text(145.75,-16.9, "Cairns", fontsize=fs, ha='right', va='top')
pylab.text(149.2,-21.2, "Mackay", fontsize=fs, ha='right', va='top')
pylab.text(153.0,-27.5, "Brisbane", fontsize=fs, ha='right', va='top')

# in which direction to depart for other points on earth and how far
# it will be to reach that destination.
# The specified point shows up as a red dot in the center of the map.

# This example shows how to use the width and height keywords
# to specify the map projection region (instead of specifying
# the lat/lon of the upper right and lower left corners).

# user enters the lon/lat of the point, and it's name
lon_0 = float(raw_input('input reference lon (degrees):'))
lat_0 = float(raw_input('input reference lat (degrees):'))
location = raw_input('name of location:')

# use these values to setup Basemap instance.
width = 28000000
m = Basemap(width=width,height=width,\
            resolution='c',projection='aeqd',\
            lat_0=lat_0,lon_0=lon_0)
# draw coasts and fill continents.
m.drawcoastlines(linewidth=0.5)
m.fillcontinents()
# 20 degree graticule.
m.drawparallels(arange(-80,81,20))
m.drawmeridians(arange(-180,180,20))
# draw a red dot at the center.
xpt, ypt = m(lon_0, lat_0)
m.plot([xpt],[ypt],'ro') 
# draw the title.
title('The World According to Garp in '+location)
show()

import pylab as p
import matplotlib.numerix as nx
from matplotlib.toolkits.basemap import Basemap as Basemap
from matplotlib.collections import LineCollection
from matplotlib.colors import rgb2hex
import random

# requires pyshapelib from Thuban (http://thuban.intevation.org/).
# cd to libraries/pyshapelib in Thuban source distribution, run
# 'python setup.py install'.

# Lambert Conformal map of lower 48 states.
m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
fig=p.figure(figsize=(8,m.aspect*8))
fig.add_axes([0.1,0.1,0.8,0.8])
# draw climate division boundaries.
shp_info = m.readshapefile('divisions','climdivs',drawbounds=True)
print shp_info
# make sure the shapefile has polygons (and not just lines).
if shp_info[1] != 5:
    print 'warning: shapefile does not contain polygons'
# choose a color for each climate division (randomly).
colors={}
divnames=[]
print m.climdivs_info[0].keys()
for shapedict in m.climdivs_info:
    divname = shapedict['ST']+repr(shapedict['DIV'])
    colors[divname] = (random.uniform(0,1),random.uniform(0,1),random.uniform(0,1))
    divnames.append(divname)
# cycle through climate divnames, color each one.

for line in file.readlines():
   l = line.replace('\n','').split()
   ul.append(float(l[0]))
   vl.append(float(l[1]))
   pl.append(float(l[2]))
u = reshape(array(ul,Float32),(nlats,nlons))
v = reshape(array(vl,Float32),(nlats,nlons))
p = reshape(array(pl,Float32),(nlats,nlons))
lats1 = -90.+dellat*arange(nlats)
lons1 = -180.+dellon*arange(nlons)
lons, lats = meshgrid(lons1, lats1)

# plot vectors in geographical (lat/lon) coordinates.

# north polar projection.
m = Basemap(lon_0=-135,boundinglat=25,
            resolution='c',area_thresh=10000.,projection='npstere')
# create a figure, add an axes.
fig=figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# rotate wind vectors to map projection coordinates.
# (also compute native map projections coordinates of lat/lon grid)
# only do Northern Hemisphere.
urot,vrot,x,y = m.rotate_vector(u[36:,:],v[36:,:],lons[36:,:],lats[36:,:],returnxy=True)
# plot filled contours over map.
cs = m.contourf(x,y,p[36:,:],15,cmap=cm.jet)
# plot wind vectors over map.
Q = m.quiver(x,y,urot,vrot) #or specify, e.g., width=0.003, scale=400)
qk = quiverkey(Q, 0.95, 1.05, 25, '25 m/s', labelpos='W')
cax = axes([0.875, 0.1, 0.05, 0.7]) # setup colorbar axes.
colorbar(cax=cax) # draw colorbar
axes(ax)  # make the original axes current again