python - 在 matplotlib basemap 投影上绘制椭圆

Question

我正在尝试在 basemap 投影上绘制椭圆。要绘制像多边形一样的圆，可以使用 tissot 函数来绘制 Tissot's indicatrix'如以下示例所示。

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)

但是，我对以 (x-x0)**2/a**2 + (y-y0)**2/2 = 1 形式绘制省略号很感兴趣。另一方面，要在规则的笛卡尔网格上绘制省略号，我可以使用以下示例代码:

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])

有没有办法在 basemap 投影上绘制一个省略号，效果类似于tissot？

Answer 1

经过数小时分析 basemap 的tissot 函数的源代码，了解ellipses 的一些属性以及大量的调试，我找到了解决问题的方法。我用一个名为 ellipse 的新函数扩展了 basemap 类，如下所示，

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

这个新功能可以像这个例子一样迅速使用:

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")

结果如下: