不能用 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()

我曾尝试修复该错误一段时间，但无济于事。我认为它源于 type(polygon) 是 sgeom.Polygon here 的错误假设。事实上，变量有时是 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 中？

Answer 1

哇哦！我终于想通了。当我更改此片段时一切正常：

      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.

中推断出来

不能用 Cartopy 裁剪湖泊

Cannot crop lakes with Cartopy

cartopy