Невозможно обрезать озера с помощью Cartopy

Question

Я пытаюсь рисовать озера на карте, используя 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() 

Я пытался исправить ошибку в течение некоторого времени, но безрезультатно. I думаю это связано с неправильным допущением, что type(polygon) is sgeom.Polygonhere. Фактически, переменная иногда имеет тип sgeom.MultiPolygon или sgeom.GeometryCollection.

В то время как мы в этом, мне кажется, что линия 544 из crs.py может использовать 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?

Answer 1

Whoa! Я, наконец, понял это. Все работает, когда я изменить этот фрагмент:

 if not self._boundary.is_ccw: 
     self._boundary.coords = list(self._boundary.coords)[::-1] 

к

 if self._boundary.is_ccw: 
     self._boundary.coords = list(self._boundary.coords)[::-1] 

, что означает, что границы проекций должны быть по часовой стрелке. Оглядываясь назад, я мог бы предположить это с lines 123–149 of crs.py.