How can I get my contour plot superimposed on a basemap

Question

This is a question I asked several months ago and am still struggling to come to a solution. My code gives me a basemap and a contour plot side by side (but printing to file only gives the contour plot), but I want them superimposed. The best solution would the one here https://gist.github.com/oblakeobjet/7546272 but this does not show how to introduce the data, and it is difficult when you are learning from scratch online. Without tiring very kind people, I hope the solution is easy as changing a line of code and that someone can help. My code

#!/usr/bin/python
# vim: set fileencoding=UTF8

import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

#fig = plt.figure(figsize=(10,8))  #when uncommented draws map with colorbar but no contours

#prepare a basemap

m = Basemap(projection = 'merc',llcrnrlon = 21, llcrnrlat = -18, urcrnrlon = 34, urcrnrlat = -8, resolution='h')

# draw country outlines.

m.drawcountries(linewidth=0.5, linestyle='solid', color='k', antialiased=1, ax=None, zorder=None)
m.drawmapboundary(fill_color = 'white')
m.fillcontinents(color='coral',lake_color='blue')
parallels = np.arange(-18, -8, 2.)
m.drawparallels(np.arange(-18, -8, 2.), color = 'black', linewidth = 0.5)
m.drawparallels(parallels,labels=[True,False,False,False])
meridians = np.arange(22,34, 2.)
m.drawmeridians(np.arange(21,36, 2.), color = '0.25', linewidth = 0.5)
m.drawmeridians(meridians,labels=[False,False,False,True])

fig = plt.figure(figsize=(10,8))       # At this position or commented draws teo figures side by side

#-- Read the data.

data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True)

#-- Now gridding data.  First making a regular grid to interpolate onto

numcols, numrows = 300, 300
xi = np.linspace(data.Lon.min(), data.Lon.max(), numcols)
yi = np.linspace(data.Lat.min(), data.Lat.max(), numrows)
xi, yi = np.meshgrid(xi, yi)

#-- Interpolating at the points in xi, yi

x, y, z = data.Lon.values, data.Lat.values, data.Z.values
zi = griddata(x, y, z, xi, yi)

#-- Display and write the results

m = plt.contourf(xi, yi, zi)
plt.scatter(data.Lon, data.Lat, c=data.Z, s=100,
       vmin=zi.min(), vmax=zi.max())
fig.colorbar(m)
plt.savefig("rainfall.jpg", format="jpg")

The plots I get look like this and

and my data

32.6  -13.6   41
27.1  -16.9   43
32.7  -10.2   46
24.2  -13.6   33
28.5  -14.4   43
28.1  -12.6   33
27.9  -15.8   46
24.8  -14.8   44
31.1  -10.2   35
25.9  -13.5   24
29.1   -9.8   10
25.8  -17.8   39
33.2  -12.3   44
28.3  -15.4   46
27.6  -16.1   47
28.9  -11.1   31
31.3   -8.9   39
31.9  -13.3   45
23.1  -15.3   31
31.4  -11.9   39
27.1  -15.0   42
24.4  -11.8   15
28.6  -13.0   39
31.3  -14.3   44
23.3  -16.1   39
30.2  -13.2   38
24.3  -17.5   32
26.4  -12.2   23
23.1  -13.5   27

Earldom

Apparently you get two figures because m is the context from which you call every map drawing command, and there after you draw in current active figure using plt based plotting commands. – Bissextile

Why the R tag? You seem pretty well set in Python. – Cataldo

Can you share the meansr.txt to run you code? – Movement

Answer 1

You're almost there, but Basemap can be temperamental, and you have to manage the z-order of plots / map details. Also, you have to transform your lon / lat coordinates to map projection coordinates before you plot them using basemap.

Here's a complete solution, which gives the following output. I've changed some colours and linewidths in order to make the whole thing more legible, YMMV. I've also scaled the size of the scatter points by the normalised 'mean' value (data['Z']) – you can simply remove it and substitute e.g. 50 if you'd prefer a constant size (it'll look like the largest marker).

Please also specify the units of rainfall, and the duration of the measurement which generated the means, if possible:

import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
%matplotlib inline

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

# grab data
data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True)
norm = Normalize()

# 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')

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

# 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)

# 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 plot
con = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')
# scatter plot
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 and title
# add colour bar, title, and scale
cbar = plt.colorbar(conf, orientation='horizontal', fraction=.057, pad=0.05)
cbar.set_label("Mean Rainfall - mm")

m.drawmapscale(
    24., -9., 28., -13,
    100,
    units='km', fontsize=10,
    yoffset=None,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#000000',
    fontcolor='#000000',
    zorder=5)

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

Using matplotlib's griddata method is convenient, but it can also be slow. As an alternative, you can use scipy's griddata methods, which are both faster, and more flexible:

from scipy.interpolate import griddata as gd

zi = gd(
    (data[['projected_lon', 'projected_lat']]),
    data.Z.values,
    (xi, yi),
    method='linear')

If you use scipy's griddata method, you'll also have to determine which of the methods (nearest, linear, cubic) give the best resulting plot.

I should add that the interpolation methods demonstrated and discussed above are the simplest possible, and aren't necessarily valid for the interpolation of rainfall data. This article gives a good overview of valid approaches and considerations for use in hydrology and hydrological modelling. The implementation of these (probably using Scipy) is left as an exercise &c.

Mastoiditis

Thank you very much for your precious time @urschrei. This will serve me well as learning material as well. I have been working on this for close to a year. I have noted all the suggestions you have made. BTW what does YMMV mean? I started with R. I was able to get contours on a map in r but I understood that one could not extrapolate (to fill the domain) in r.Though am still not extrapolating I think I can be happy living with python. Thanks so much for all kind and helpful people. – Earldom

@ZiloreMumba YMMV means 'your mileage may vary' – that you might feel differently about the colour map and line width choices I made, and that you should feel free to experiment with these. Good luck! – Mastoiditis

@AndreAraujo It's just the original data, saved as a txt file, with a header row: Lon Lat Values. The file is whitespace-delimited. – Mastoiditis

Answer 2

I have not everything installed here to run your code, but you should try plotting to the basemap m you created, like this:

# fig = plt.figure(figsize=(10,8)) # omit this at line 28

(...)

m.contourf(xi, yi, zi)
m.scatter(data.Lon, data.Lat, c=data.Z, s=100,
   vmin=zi.min(), vmax=zi.max())

(please tell if this doesn't work)

Bissextile

replacing my plotting functions withe snippet you sent plots the map without contours nor scatter plot. – Earldom

Thanks all those who spent their time reading my code (above), and those who helped me make this code to work (six months ago), particularly @urschrei. Since then I have been wondering if there is any extrapolation in python. I know extrapolation can give weird results, but I have not seen anything in documentation. Is it possible to extrapolate to at least fill the map? – Earldom

I am surprised there is no answer to this "simple" question: Is there a command for extrapolation in python? – Earldom

@ZiloreMumba you could use, for example, scypy.rbf module to create an RBF interpolator for a set of data. You could then interpolate further away from the original dataset (that is, extrapolate), but check results, because the farther away, the more error you get. And also, check out other interpolators, not only RBF (there is also kriging, but not natively in scipy I think. – Bissextile