How can I get the contour portion superimposed on the base map

Question

This is a question that I asked a few months ago, and I'm still trying to come to a solution. My code gives me a basemap and a contour plot side by side (but printing to a file gives only a contour plot), but I want them to overlap. A better solution would be here https://gist.github.com/oblakeobjet/7546272 , but it doesn’t show how to enter data, and it is difficult when you study for scratches on the Internet. Without the tedious very kind people, I hope the solution is as 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 graphs 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 

Answer 1

You are almost there, but Basemap can be temperamental, and you need to control the z-order of the details / data of the map. In addition, you must convert your lon / lat coordinates to match projection coordinates before creating them using the base map.

Here is a complete solution that gives the following result. I changed some colors and line widths to make it all more legible, YMMV. I also scaled the size of the scatter points by the normalized "average" value ( data['Z'] ) - you can simply delete it and replace it, for example. 50 if you prefer a constant size (it will look like the largest marker).

Please also indicate the units of precipitation and the duration of the measurement in which the funds were received, 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 the matplotlib griddata method is convenient, but it can also be slow. Alternatively, you can use scipy griddata methods, which are 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 the scipy griddata method, you will also need to determine which of the methods ( nearest , linear , cubic ) gives the best result.

I must add that the interpolation methods demonstrated and discussed above are the simplest and not necessarily valid for interpolating rain data. This article provides a good overview of valid approaches and considerations for use in hydrology and hydrological modeling. The implementation of these (possibly using Scipy) is left as an exercise & c.

Answer 2

I do not have everything installed here to run your code, but you should try building based on the base map m , 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()) 

(tell me if this works)