我正在尝试使用 Cartopy 0.14 和 Shapely 1.5.12 在地图上绘制湖泊。使用我的自定义投影，保存或显示图形有时会失败，堆栈跟踪以结尾
File "/usr/local/lib/python2.7/dist-packages/Cartopy-0.14.dev0-py2.7-linux-x86_64.egg/cartopy/crs.py", line 291, in _project_multipolygon
r = self._project_polygon(geom, src_crs)
File "/usr/local/lib/python2.7/dist-packages/Cartopy-0.14.dev0-py2.7-linux-x86_64.egg/cartopy/crs.py", line 330, in _project_polygon
return self._rings_to_multi_polygon(rings, is_ccw)
File "/usr/local/lib/python2.7/dist-packages/Cartopy-0.14.dev0-py2.7-linux-x86_64.egg/cartopy/crs.py", line 589, in _rings_to_multi_polygon
multi_poly = sgeom.MultiPolygon(polygon_bits)
File "/usr/local/lib/python2.7/dist-packages/Shapely-1.5.12-py2.7-linux-x86_64.egg/shapely/geometry/multipolygon.py", line 62, in __init__
self._geom, self._ndim = geos_multipolygon_from_polygons(polygons)
File "/usr/local/lib/python2.7/dist-packages/Shapely-1.5.12-py2.7-linux-x86_64.egg/shapely/geometry/multipolygon.py", line 178, in geos_multipolygon_from_polygons
geom, ndims = polygon.geos_polygon_from_py(shell, holes)
File "/usr/local/lib/python2.7/dist-packages/Shapely-1.5.12-py2.7-linux-x86_64.egg/shapely/geometry/polygon.py", line 503, in geos_polygon_from_py
geos_shell, ndim = geos_linearring_from_py(shell)
File "shapely/speedups/_speedups.pyx", line 214, in shapely.speedups._speedups.geos_linearring_from_py (shapely/speedups/_speedups.c:3679)
ValueError: A LinearRing must have at least 3 coordinate tuples
当湖泊的边界与投影的边界相交时，就会发生这种情况。我无法使用内置的 Cartopy 投影重现该行为。这是我能想到的最小测试用例：
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from matplotlib import pyplot as plt
import numpy as np
from shapely import geometry as sgeom
class Polyconic(ccrs.Projection):
NUM_BOUNDARY_SEGMENTS = 30
def __init__(self, central_longitude, globe=None):
proj4_params = [
('proj', 'poly'),
('lon_0', central_longitude)]
super(Polyconic, self).__init__(proj4_params, globe=globe)
bounds = self.ToPolygon(self.GetLimits(central_longitude)).bounds
self._x_limits = bounds[0], bounds[2]
self._y_limits = bounds[1], bounds[3]
self._boundary = self.ToPolygon(self.GetDomain(central_longitude)).exterior
if not self._boundary.is_ccw:
self._boundary.coords = list(self._boundary.coords)[::-1]
@staticmethod
def GetDomain(central_longitude):
lats = np.linspace(0, +90, Polyconic.NUM_BOUNDARY_SEGMENTS + 1)
lons = np.linspace(
central_longitude - 15., central_longitude + 15.,
Polyconic.NUM_BOUNDARY_SEGMENTS + 1)
domain = []
for lat in lats:
domain.append((central_longitude - 15., lat))
for lat in reversed(lats):
domain.append((central_longitude + 15., lat))
return domain
@staticmethod
def GetLimits(central_longitude):
return [
(central_longitude - 15., 0.),
(central_longitude + 15., 0.),
(central_longitude + 15., +90.),
(central_longitude - 15., +90.)]
def ToPolygon(self, polygon):
return sgeom.Polygon(self.transform_points(
ccrs.PlateCarree(),
np.array([p[0] for p in polygon]),
np.array([p[1] for p in polygon])))
@property
def threshold(self):
return 1e3
@property
def boundary(self):
return self._boundary
@property
def x_limits(self):
return self._x_limits
@property
def y_limits(self):
return self._y_limits
plt.figure()
# ax = plt.axes(projection=Polyconic(180)) works.
ax = plt.axes(projection=Polyconic(0))
lakes = cfeature.NaturalEarthFeature('physical', 'lakes', '50m')
ax.add_feature(lakes)
plt.show()
我试图修复错误一段时间，但无济于事。我认为type(polygon)它源于sgeom.Polygon 这里的错误假设。事实上，变量有时是sgeom.MultiPolygon或类型sgeom.GeometryCollection。
在我看来，crs.py 的第 544 行可能会用到prep_polygon，第 562-577 行可以简化如下：
y4 += by
box = sgeom.box(x3, y3, x4, y4)
for ring in interior_rings:
polygon = sgeom.Polygon(ring)
if polygon.is_valid:
# Invert the polygon
polygon = box.difference(polygon)
我的问题是：我的代码或 Cartopy 中的错误是什么？