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.

Hyperbaric

Trying to figure out the same thing myself. Any luck? – Augustina

This will be a Christmas project for me. So if you can wait until 2015... – Hyperbaric

You actually want a surface, not a (potentially very dense) set of points? Do you mean the type of function that interp2d outputs? – Augustina

Yes, a dense set of points would be fine at this stage. I don't seem to be able to find the time to do this. I should delete the question. I will set myself a deadline for April 1st for this. – Hyperbaric

Check out this example: webrian.ch/2012/11/heat-maps-with-python.html – Augustina

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

Undamped

@Hyperbaric does this give you what you want, or is there something missing that here? Feel free to ping me if there's something that doesn't add up for you. – Undamped

I will give you the bounty. If it expires I will issue it again. I just want to provide an example based on your answer with more data points using bicubic interpolation. Your example fails with 16 data points or more. Sorry about the delay. I'll get it done this weekend. Thanks for your excellent answer, which I have already up voted. – Hyperbaric

Cool - no problem. If you've got 16 data points (or can point me where to get them), I'm happy to have a look at it with that dataset. I can't work out what the data are derived from, but I'm no expert on African cities. I'll revisit even if you award this bounty - it's an interesting problem and I don't want you to waste more of your rep than you need. – Undamped

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à:

Downtown

This solution is not without issues. First of all it is sloooow, then the small landmasses on the edge of the map are not properly rendered and finally it is not tested for robustness, and it not work for other projection. I will try to improve those aspects and update the answer, but feel free to come up with suggestions. – Downtown

What is descartes? I don't seem to have that. – Hyperbaric

Lagos seems to be cyan instead of red on your map. – Hyperbaric

Well, the interpolation raster may be different from yours. My focus was to mask the oceans. As for the bi-cubic interpolation, I do not think 4 points are enough data to allow you to do that. You will need a larger set or a linear interpolation. Does that answer your comment? – Downtown