Dibujo de elipses en las proyecciones del mapa base matplotlib

Question

Estoy tratando de dibujar elipsis en una proyección de mapa base. Para dibujar un círculo como un polígono, se usa la función tissot para dibujar Tissot's indicatrix' como se ilustra en el siguiente ejemplo.

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) 

Sin embargo, estoy interesado en el trazado de una elipsis en forma (x-x0)**2/a**2 + (y-y0)**2/2 = 1. Por otro lado, para dibujar una elipse en una cuadrícula cartesiana regulares puedo utilizar el siguiente código de ejemplo:

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

¿Hay una manera de trazar una elipse en una proyección de mapa base con el un efecto similar al tissot?

Answer 1

Después de horas analizando el código fuente de la función tissot del mapa base, aprendiendo algunas propiedades de ellipses y mucha depuración, vine con una solución a mi problema. He ampliado la clase mapa base con una nueva función llamada ellipse de la siguiente manera,

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 

Esta nueva función se puede utilizar de inmediato como en este ejemplo:

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

que tiene el siguiente resultado: