#!/usr/bin/env python

import matplotlib.pyplot as plt           # pyplot module import
from mpl_toolkits.basemap import Basemap  # basemap import
import numpy as np                        # Numpy import

# Cities names and coordinates
cities = ['London', 'New York', 'Madrid', 'Cairo', 'Moscow', 'Delhi', 'Dakar']
lat = [51.507778, 40.716667, 40.4, 30.058, 55.751667, 28.61, 14.692778]
lon = [-0.128056, -74, -3.683333, 31.229, 37.617778, 77.23, -17.446667]

# orthogonal projection of the Earth
m = Basemap(projection='ortho', lat_0=45, lon_0=10, resolution='l',area_thresh=1000)

m.bluemarble()

## draw the borders of the map
#m.drawmapboundary()
## draw the coasts borders and fill the continents
#m.drawcoastlines()
#m.fillcontinents(color='coral')

# map city coordinates to map coordinates
x, y = m(lon, lat)

# draw a red dot at cities coordinates
plt.plot(x, y, 'ro')

# for each city,
for city, xc, yc in zip(cities, x, y):
    # draw the city name in a yellow (shaded) box
    plt.text(xc+250000, yc-150000, city)

# make up some data on a regular lat/lon grid.
nlats = 73; nlons = 145; delta = 2.*np.pi/(nlons-1)
lats = (0.5*np.pi-delta*np.indices((nlats,nlons))[0,:,:])
lons = (delta*np.indices((nlats,nlons))[1,:,:])
wave = 0.75*(np.sin(2.*lats)**8*np.cos(4.*lons))
mean = 0.5*np.cos(2.*lats)*((np.sin(2.*lats))**2 + 2.)
# compute native map projection coordinates of lat/lon grid.
x, y = m(lons*180./np.pi, lats*180./np.pi)
# contour data over the map.
CS = m.contour(x,y,wave+mean,15,linewidths=1.5)

plt.savefig('fig_dataOverMap2.png')
plt.show()