Contours with map overlay on an irregular grid in python

Question

Here are my details:

Lon Lat Z Z2 pos 32.6 -13.6 41 9 CHIP 27.1 -16.9 43 12 CHOM 32.7 -10.2 46 14 ISOK 24.2 -13.6 33 13 KABO 28.5 -14.4 43 11 KABW 28.1 -12.6 33 16 KAFI 27.9 -15.8 46 13 KAFU 24.8 -14.8 44 9 KAOM 31.1 -10.2 35 14 KASA 25.9 -13.5 24 8 KASE 29.1 -9.8 10 13 KAWA 25.8 -17.8 39 11 LIVI 33.2 -12.3 44 8 LUND 28.3 -15.4 46 12 LUSA 27.6 -16.1 47 9 MAGO 28.9 -11.1 31 15 MANS 31.3 -8.9 39 9 MBAL 31.9 -13.3 45 9 MFUW 23.1 -15.3 31 9 MONG 31.4 -11.9 39 9 MPIK 27.1 -15.0 42 12 MUMB 24.4 -11.8 15 9 MWIN 28.6 -13.0 39 9 NDOL 31.3 -14.3 44 12 PETA 23.3 -16.1 39 5 SENA 30.2 -13.2 38 11 SERE 24.3 -17.5 32 10 SESH 26.4 -12.2 23 12 SOLW 23.1 -13.5 27 14 ZAMB 

And here is the code I tried. The part that causes me the problem is the grid of data and loops. The code below is mostly related to the various posts in stackoverflow, but it is not easy to find anything regarding my specific question.

 #!/usr/bin/python from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import griddata #prepare a basemap m = Basemap(projection = 'merc',llcrnrlon = 21, llcrnrlat = -18, urcrnrlon = 34, urcrnrlat = -8) m.drawmapboundary(fill_color = 'white') m.fillcontinents(color = '0.85') # draw country outlines. m.drawcountries(linewidth=0.5, linestyle='solid', color='k', antialiased=1, ax=None, zorder=None) m.drawparallels(np.arange(-18, -8, 2.), color = '0.25', linewidth = 0.5) m.drawparallels(parallels,labels=[True,False,False,False]) m.drawmeridians(np.arange(21,36, 2.), color = '0.25', linewidth = 0.5) m.drawmeridians(meridians,labels=[False,False,False,True]) #plt.show() ##Code works up to here, drawing basemap ####################################################### #prepare data for contouring #declare empty arrays for xyz xarray = [] yarray = [] zarray = [] #get data and pack in the arrays with open("meansr.txt", "r") as f: for line in f: (x, y, z, z2, pos) = line.split() xarray.append(x) yarray.append(y) zarray.append(z) f.close() #prepare and grid the data extent = (min(xarray), max(xarray), min(yarray), max(yarray)) xs,ys = np.mgrid[extent[0]:extent[1]:30j, extent[2]:extent[3]:30j] z = griddata((xarray, yarray), zarray, (xs, ys)) x,y = m(xarray,yarray) #make contour plot cs = m.contourf(x, y, z) cs2 = m.contour(x,y,z, levels = range(10,50,5),colors = 'blue') plt.clabel(cs, fmt = '%.0f', inline = True) plt.show() 

Here is the error I get on the contour part:

 Traceback (most recent call last): File "./contour12.py", line 44, in <module> xs,ys = np.mgrid[extent[0]:extent[1]:29j, extent[2]:extent[3]:29j] File "/home/zmumba/anaconda/lib/python2.7/site-packages/numpy/lib/index_tricks.py", line 185, in __getitem__ step = key.step AttributeError: 'tuple' object has no attribute 'step' 

I don't know what 30j should do in the above code, although my data has 29 points

Answer 1

To begin with, ignore some of the things based on the map, and simply treat your lat, long coordinates as a Cartesian coordinate system.

 import numpy as np import pandas as pd from matplotlib.mlab import griddata import matplotlib.pyplot as plt #-- Read the data. # I'm going to use `pandas` to read in and work with your data, mostly due to # the text site names. Using pandas is optional, however. data = pd.read_csv('your_data.txt', delim_whitespace=True) #-- Now let grid your data. # First we'll make a regular grid to interpolate onto. This is equivalent to # your call to `mgrid`, but it broken down a bit to make it easier to # understand. The "30j" in mgrid refers to 30 rows or columns. numcols, numrows = 30, 30 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) #-- Interpolate at the points in xi, yi # "griddata" expects "raw" numpy arrays, so we'll pass in # data.x.values instead of just the pandas series data.x x, y, z = data.Lon.values, data.Lat.values, data.Z.values zi = griddata(x, y, z, xi, yi) #-- Display the results fig, ax = plt.subplots() im = ax.contourf(xi, yi, zi) ax.scatter(data.Lon, data.Lat, c=data.Z, s=100, vmin=zi.min(), vmax=zi.max()) fig.colorbar(im) plt.show() 

The “block” border is due to the coarse (30x30) resolution of the grid. griddata uses the triangulation method, so nothing outside the convex hull of your data points is interpolated. To see this more clearly, raise the number and number to 300x300:

You can also use several other interpolation methods (especially if you want to extend the interpolation beyond the convex hull of the data).