matplotlibのベースマップ投影の楕円を描く

Question

ベースマップ投影に楕円を描画しようとしています。多角形のような円を描画するには、次の例のようにTissot's indicatrix'を描画するために使用されるtissot関数があります。

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

tissotと同様の効果を持つベースマップの投影に省略記号をプロットする方法はあります：一方で、私は次のサンプルコードを使用することができ、通常のデカルト格子上の省略記号を描画するには？

Answer 1

ellipsesのいくつかのプロパティと多くのデバッグを学んで、ベースマップのtissot関数のソースコードを分析した後、私は私の問題を解決しました。

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

次のような結果を有する：私は

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 

この新機能は、この例のよう速やかに使用することができ、次のようにellipseと呼ばれる新しい機能を備えたベースマップクラスを拡張しました