color matplotlib map using bicubic interpolation

Question

I know that matplotlib and scipy can do bicubic interpolation: http://matplotlib.org/examples/pylab_examples/image_interp.html http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

I also know that it is possible to draw a map of the world with matplotlib: http://matplotlib.org/basemap/users/geography.html http://matplotlib.org/basemap/users/examples.html http://matplotlib.org/basemap/api/basemap_api.html

But can I do a bicubic interpolation based on 4 data points and only color the land mass?

For example using these for 4 data points (longitude and latitude) and colors:

Lagos: 6.453056, 3.395833; red HSV 0 100 100 (or z = 0)
Cairo: 30.05, 31.233333; green HSV 90 100 100 (or z = 90)
Johannesburg: -26.204444, 28.045556; cyan HSV 180 100 100 (or z = 180)
Mogadishu: 2.033333, 45.35; purple HSV 270 100 100 (or z = 270)

I am thinking that it must be possible to do the bicubic interpolation across the range of latitudes and longitudes and then add oceans, lakes and rivers on top of that layer? I can do this with drawmapboundary. Actually there is an option maskoceans for this: http://matplotlib.org/basemap/api/basemap_api.html#mpl_toolkits.basemap.maskoceans

I can interpolate the data like this:

xnew, ynew = np.mgrid[-1:1:70j, -1:1:70j]
tck = interpolate.bisplrep(x, y, z, s=0)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

Or with scipy.interpolate.interp2d: http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

Here it is explained how to convert to map projection coordinates: http://matplotlib.org/basemap/users/mapcoords.html

But I need to figure out how to do this for a calculated surface instead of individual points. Actually there is an example of such a topographic map using external data, which I should be able to replicate: http://matplotlib.org/basemap/users/examples.html

P.S. I am not looking for a complete solution. I would much prefer to solve this myself. Rather I am looking for suggestions and hints. I have been using gnuplot for more than 10 years and only switched to matplotlib within the past few weeks, so please don't assume I know even the simplest things about matplotlib.

Answer 1

I think this is what you are looking for (roughly). Note the crucial things are masking the data array before you plot the pcolor and passing in the hsv colormap (Docs: cmap parameter for pcolormesh and available colormaps).

I've kept the code for plotting the maps quite close to the examples so it should be easy to follow. I've kept your interpolation code for the same reason. Note that the interpolation is linear rather than cubic - kx=ky=1 - because you don't give enough points to do cubic interpolation (you'd need at least 16 - scipy will complain with less saying that "m must be >= (kx+1)(ky+1)", although the constraint is not mentioned in the documentation).

I've also extended the range of your meshgrid and kept in lat / lon for x and y throughout.

Code

from mpl_toolkits.basemap import Basemap,maskoceans
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
# set up orthographic map projection with
# perspective of satellite looking down at 0N, 20W (Africa in main focus)
# use low resolution coastlines.
map = Basemap(projection='ortho',lat_0=0,lon_0=20,resolution='l')

# draw coastlines, country boundaries
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
# Optionally (commented line below) give the map a fill colour - e.g. a blue sea
#map.drawmapboundary(fill_color='aqua')
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))


data = {'Lagos': (6.453056, 3.395833,0),
        'Cairo': (30.05, 31.233333,90),
        'Johannesburg': (-26.204444, 28.045556,180),
        'Mogadishu': (2.033333, 45.35, 270)}

x,y,z = zip(*data.values())

xnew, ynew = np.mgrid[-30:60:0.1, -50:50:0.1]

tck = interpolate.bisplrep(x, y, z, s=0,kx=1,ky=1)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)
znew = maskoceans(xnew, ynew, znew)

col_plot = map.pcolormesh(xnew, ynew, znew, latlon=True, cmap='hsv')
plt.show()

Output

Answer 2

Observe that doing the opposite, that is putting a raster on the sea and lay a mask over the continents, is easy as pie. Simply use map.fillcontinents(). So the basic idea of this solution is to modify the fillcontinents function so that it lays polygons over the oceans.

The steps are:

1. Create a large circle-like polygon that covers the entire globe.
2. Create a polygon for each shape in the map.coastpolygons array.
3. Cut the shape of the landmass polygon away from the circle using shapely and its difference method.
4. Add the remaining polygons, which have the shape of the oceans, on the top, with a high zorder.

The code:

from mpl_toolkits.basemap import Basemap
import numpy as np
from scipy import interpolate
from shapely.geometry import Polygon
from descartes.patch import PolygonPatch

def my_circle_polygon( (x0, y0), r, resolution = 50 ):
    circle = []
    for theta in np.linspace(0,2*np.pi, resolution):
        x = r * np.cos(theta) + x0
        y = r * np.sin(theta) + y0    
        circle.append( (x,y) )

    return Polygon( circle[:-1] )

def filloceans(the_map, color='0.8', ax=None):
    # get current axes instance (if none specified).
    if not ax:     
        ax = the_map._check_ax()

    # creates a circle that covers the world
    r =  0.5*(map.xmax - map.xmin) # + 50000 # adds a little bit of margin
    x0 = 0.5*(map.xmax + map.xmin)
    y0 = 0.5*(map.ymax + map.ymin)
    oceans = my_circle_polygon( (x0, y0) , r, resolution = 100 )

    # for each coastline polygon, gouge it out of the circle
    for x,y in the_map.coastpolygons:
        xa = np.array(x,np.float32)
        ya = np.array(y,np.float32)

        xy = np.array(zip(xa.tolist(),ya.tolist()))
        continent = Polygon(xy)

        ##  catches error when difference with lakes 
        try:
            oceans = oceans.difference(continent)
        except: 
            patch = PolygonPatch(continent, color="white", zorder =150)
            ax.add_patch( patch ) 

    for ocean in oceans:
        sea_patch = PolygonPatch(ocean, color="blue", zorder =100)
        ax.add_patch( sea_patch )

###########  DATA
x = [3.395833, 31.233333, 28.045556, 45.35   ]
y = [6.453056, 30.05,    -26.204444, 2.033333]
z = [0, 90, 180, 270]

# set up orthographic map projection
map = Basemap(projection='ortho', lat_0=0, lon_0=20, resolution='l')

## Plot the cities on the map
map.plot(x,y,".", latlon=1)

# create a interpolated mesh and set it on the map
interpol_func = interpolate.interp2d(x, y, z, kind='linear')
newx = np.linspace( min(x), max(x) )
newy = np.linspace( min(y), max(y) )
X,Y = np.meshgrid(newx, newy)
Z = interpol_func(newx, newy)
map.pcolormesh( X, Y, Z, latlon=1, zorder=3)

filloceans(map, color="blue")

Voilà: