经过数小时分析basemap的tissot
函数源代码,学习了椭圆的一些属性以及进行了大量调试后,我找到了解决问题的方法。我扩展了basemap类并添加了一个名为ellipse
的新函数,如下所示:
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)
azf, azb, dist = g.inv([x0, x0], [y0, y0], [x0+a, x0], [y0, y0+b])
tsid = dist[0] * dist[1]
seg = [self(x0+a, y0)]
AZ = numpy.linspace(azf[0], 360. + azf[0], n)
for i, az in enumerate(AZ):
if numpy.allclose(0., y0) and (numpy.allclose(90., az) or
numpy.allclose(270., az)):
continue
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)
if x < 1e20 and y < 1e20:
seg.append((x, y))
poly = Polygon(seg, **kwargs)
ax.add_patch(poly)
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')
m.drawparallels(numpy.arange(-80.,81.,20.))
m.drawmeridians(numpy.arange(-180.,181.,20.))
m.drawmapboundary(fill_color='aqua')
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")
以下是结果:
这将产生以下结果:![示例图片](https://istack.dev59.com/AOFCy.webp)