Map interpolated values using a contour plot, and a scatter plot

contour.ipynb

contour.py

# coding: utf-8

# In[144]:

import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from mpl_toolkits.basemap import Basemap, interp
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
mpl.rcParams['figure.figsize'] = (16, 12)
get_ipython().magic(u'matplotlib inline')


# In[96]:

dd = {"Lon": {0: 32.600000000000001,
  1: 27.100000000000001,
  2: 32.700000000000003,
  3: 24.199999999999999,
  4: 28.5,
  5: 28.100000000000001,
  6: 27.899999999999999,
  7: 24.800000000000001,
  8: 31.100000000000001,
  9: 25.899999999999999,
  10: 29.100000000000001,
  11: 25.800000000000001,
  12: 33.200000000000003,
  13: 28.300000000000001,
  14: 27.600000000000001,
  15: 28.899999999999999,
  16: 31.300000000000001,
  17: 31.899999999999999,
  18: 23.100000000000001,
  19: 31.399999999999999,
  20: 27.100000000000001,
  21: 24.399999999999999,
  22: 28.600000000000001,
  23: 31.300000000000001,
  24: 23.300000000000001,
  25: 30.199999999999999,
  26: 24.300000000000001,
  27: 26.399999999999999,
  28: 23.100000000000001},
 "Lat": {0: -13.6,
  1: -16.899999999999999,
  2: -10.199999999999999,
  3: -13.6,
  4: -14.4,
  5: -12.6,
  6: -15.800000000000001,
  7: -14.800000000000001,
  8: -10.199999999999999,
  9: -13.5,
  10: -9.8000000000000007,
  11: -17.800000000000001,
  12: -12.300000000000001,
  13: -15.4,
  14: -16.100000000000001,
  15: -11.1,
  16: -8.9000000000000004,
  17: -13.300000000000001,
  18: -15.300000000000001,
  19: -11.9,
  20: -15.0,
  21: -11.800000000000001,
  22: -13.0,
  23: -14.300000000000001,
  24: -16.100000000000001,
  25: -13.199999999999999,
  26: -17.5,
  27: -12.199999999999999,
  28: -13.5},
 "Z": {0: 41,
  1: 43,
  2: 46,
  3: 33,
  4: 43,
  5: 33,
  6: 46,
  7: 44,
  8: 35,
  9: 24,
  10: 10,
  11: 39,
  12: 44,
  13: 46,
  14: 47,
  15: 31,
  16: 39,
  17: 45,
  18: 31,
  19: 39,
  20: 42,
  21: 15,
  22: 39,
  23: 44,
  24: 39,
  25: 38,
  26: 32,
  27: 23,
  28: 27}}


# In[159]:

# uncomment to get from CSV
# data = pd.read_csv(
#     'means.tsv',
#     delim_whitespace=True, header=None,
#     names=["Lon", "Lat", "Z"])
data = pd.DataFrame(dd)


# In[160]:

# define map extent
lllon = 21
lllat = -18
urlon = 34
urlat = -8

# set up Basemap instance
m = Basemap(
    projection = 'merc',
    llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat,
    resolution='h')


# In[187]:

# transform lon / lat coordinates to map projection
data['projected_lon'], data['projected_lat'] = m(*(data["Lon"].values, data["Lat"].values))
norm = Normalize()

# grid data
numcols, numrows = 1000, 1000
xi = np.linspace(data['projected_lon'].min(), data['projected_lon'].max(), numcols)
yi = np.linspace(data['projected_lat'].min(), data['projected_lat'].max(), numrows)
xi, yi = np.meshgrid(xi, yi)

# interpolate
x, y, z = data['projected_lon'].values, data['projected_lat'].values, data['Z'].values
zi = griddata(x, y, z, xi, yi)


# In[185]:

# set up plot
plt.clf()
fig = plt.figure(figsize=(6.4, 5.12))
ax = fig.add_subplot(111, axisbg='w', frame_on=False)

# draw map details
m.drawmapboundary(fill_color = 'white')
m.fillcontinents(color='#C0C0C0', lake_color='#7093DB')
m.drawcountries(
    linewidth=.75, linestyle='solid', color='#000073',
    antialiased=True,
    ax=ax, zorder=3)

m.drawparallels(
    np.arange(lllat, urlat, 2.),
    color = 'black', linewidth = 0.5,
    labels=[True, False, False, False])
m.drawmeridians(
    np.arange(lllon, urlon, 2.),
    color = '0.25', linewidth = 0.5,
    labels=[False, False, False, True])

# contour plots
con = m.contour(xi, yi, zi, 15, zorder=4, linewidths=.25, linestyles='dashed', colors='k', alpha=0.6)
conf = m.contourf(xi, yi, zi, 15, zorder=4, alpha=0.6, cmap='RdPu')
# scatter plot - vmin/max for colormap compat
m.scatter(
    data['projected_lon'],
    data['projected_lat'],
    color='#545454',
    edgecolor='#ffffff',
    alpha=.75,
    s=50 * norm(data['Z']),
    cmap='RdPu',
    ax=ax,
    vmin=zi.min(), vmax=zi.max(), zorder=4)

# add colour bar, title, and scale
cbar = plt.colorbar(orientation='horizontal', fraction=.057, pad=0.05)
plt.title("Mean Rainfall")
m.drawmapscale(
    24., -9., 28., -13,
    100,
    units='km', fontsize=7,
    yoffset=None,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#000000',
    fontcolor='#000000',
    zorder=5)

plt.savefig("rainfall.png", format="png", transparent=True, dpi=300)
plt.show()


# Also look at http://earthpy.org/interpolation_between_grids_with_basemap.html for grid interpolation