Drawing ellipses on matplotlib basemap projections

Question

I am trying to draw ellipses on a basemap projection. To draw a circle like polygon there is the tissot function used to draw Tissot's indicatrix' as illustrates the following example.

from mpl_toolkits.basemap import Basemap

x0, y0 = 35, -50
R = 5

m = Basemap(width=8000000,height=7000000, resolution='l',projection='aea', 
    lat_1=-40.,lat_2=-60,lon_0=35,lat_0=-50)
m.drawcoastlines()
m.tissot(x0, y0, R, 100, facecolor='g', alpha=0.5)

However, I am interested in plotting an ellipsis in the form (x-x0)**2/a**2 + (y-y0)**2/2 = 1. On the other hand, to draw an ellipsis on a regular Cartesian grid I can use the following sample code:

import pylab
from matplotlib.patches import Ellipse

fig = pylab.figure()
ax = pylab.subplot(1, 1, 1, aspect='equal')

x0, y0 = 35, -50
w, h = 10, 5
e = Ellipse(xy=(x0, y0), width=w, height=h, linewidth=2.0, color='g')
ax.add_artist(e)
e.set_clip_box(ax.bbox)
e.set_alpha(0.7)
pylab.xlim([20, 50])
pylab.ylim([-65, -35])

Is there a way to plot an ellipsis on a basemap projection with the an effect similar to tissot?

Answer 1

After hours analyzing the source code of basemap's tissot function, learning some properties of ellipses and lot's of debugging, I came with a solution to my problem. I've extended the basemap class with a new function called ellipse as follows,

from __future__ import division
import pylab
import numpy

from matplotlib.patches import Polygon
from mpl_toolkits.basemap import pyproj
from mpl_toolkits.basemap import Basemap

class Basemap(Basemap):
    def ellipse(self, x0, y0, a, b, n, ax=None, **kwargs):
        """
        Draws a polygon centered at ``x0, y0``. The polygon approximates an
        ellipse on the surface of the Earth with semi-major-axis ``a`` and 
        semi-minor axis ``b`` degrees longitude and latitude, made up of 
        ``n`` vertices.

        For a description of the properties of ellipsis, please refer to [1].

        The polygon is based upon code written do plot Tissot's indicatrix
        found on the matplotlib mailing list at [2].

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.patches.Polygon

        RETURNS
            poly : a maptplotlib.patches.Polygon object.

        REFERENCES
            [1] : http://en.wikipedia.org/wiki/Ellipse


        """
        ax = kwargs.pop('ax', None) or self._check_ax()
        g = pyproj.Geod(a=self.rmajor, b=self.rminor)
        # Gets forward and back azimuths, plus distances between initial
        # points (x0, y0)
        azf, azb, dist = g.inv([x0, x0], [y0, y0], [x0+a, x0], [y0, y0+b])
        tsid = dist[0] * dist[1] # a * b

        # Initializes list of segments, calculates \del azimuth, and goes on 
        # for every vertex
        seg = [self(x0+a, y0)]
        AZ = numpy.linspace(azf[0], 360. + azf[0], n)
        for i, az in enumerate(AZ):
            # Skips segments along equator (Geod can't handle equatorial arcs).
            if numpy.allclose(0., y0) and (numpy.allclose(90., az) or
                numpy.allclose(270., az)):
                continue

            # In polar coordinates, with the origin at the center of the 
            # ellipse and with the angular coordinate ``az`` measured from the
            # major axis, the ellipse's equation  is [1]:
            #
            #                           a * b
            # r(az) = ------------------------------------------
            #         ((b * cos(az))**2 + (a * sin(az))**2)**0.5
            #
            # Azymuth angle in radial coordinates and corrected for reference
            # angle.
            azr = 2. * numpy.pi / 360. * (az + 90.)
            A = dist[0] * numpy.sin(azr)
            B = dist[1] * numpy.cos(azr)
            r = tsid / (B**2. + A**2.)**0.5
            lon, lat, azb = g.fwd(x0, y0, az, r)
            x, y = self(lon, lat)

            # Add segment if it is in the map projection region.
            if x < 1e20 and y < 1e20:
                seg.append((x, y))

        poly = Polygon(seg, **kwargs)
        ax.add_patch(poly)

        # Set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

        return poly

This new function can be used promptly like in this example:

pylab.close('all')
pylab.ion()

m = Basemap(width=12000000, height=8000000, resolution='l', projection='stere',
            lat_ts=50, lat_0=50, lon_0=-107.)
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(numpy.arange(-80.,81.,20.))
m.drawmeridians(numpy.arange(-180.,181.,20.))
m.drawmapboundary(fill_color='aqua') 
# draw ellipses
ax = pylab.gca()
for y in numpy.linspace(m.ymax/20, 19*m.ymax/20, 9):
    for x in numpy.linspace(m.xmax/20, 19*m.xmax/20, 12):
        lon, lat = m(x, y, inverse=True)
        poly = m.ellipse(lon, lat, 3, 1.5, 100, facecolor='green', zorder=10,
            alpha=0.5)
pylab.title("Ellipses on stereographic projection")

Which has the following outcome: