title: Fix up shapely polygon object when discontinuous after map projection
body: <p>This demo program (intended to be run in an IPython notebook; you need <code>matplotlib</code>, <code>mpl_toolkits.basemap</code>, <code>pyproj</code>, and <code>shapely</code>) is supposed to plot increasingly large circles on the surface of the Earth.  It works correctly as long as the circle does not cross over one of the poles.  If that happens, the result is complete nonsense when plotted on a map (see below cell 2)</p>
  
  <p>If I plot them "in a void" instead of on a map (see below cell 3) the results are correct in the sense that, if you removed the horizontal line going from +180 to -180 longitude, the rest of the curve would indeed delimit the boundary between the interior and exterior of the desired circle.  However, they are <em>wrong</em> in that the polygon is invalid (<code>.is_valid</code> is False), and much more importantly, the nonzero-winding-number interior of the polygon does <em>not</em> enclose the correct region of the map.</p>
  
  <p>I believe this is happening because <code>shapely.ops.transform</code> is blind to the coordinate singularity at +180==-180 longitude.  The <em>question</em> is, how do I detect the problem and repair the polygon, so that it does enclose the correct region of the map?  In this case, an appropriate fixup would be to replace the horizontal segment from (X,+180) -- (X,-180) with three lines, (X,+180) -- (+90,+180) -- (+90,-180) -- (X,-180); but note that if the circle had gone over the <em>south</em> pole, the fixup lines would need to go south instead.  And if the circle had gone over <em>both</em> poles, we'd have a valid polygon again but its interior would be the complement of what it should be.  I need to detect all of these cases and handle them correctly.  Also, I do not know how to "edit" a shapely geometry object.</p>
  
  <p>Downloadable notebook: <a href="https://gist.github.com/zackw/e48cb1580ff37acfee4d0a7b1d43a037" rel="nofollow noreferrer">https://gist.github.com/zackw/e48cb1580ff37acfee4d0a7b1d43a037</a></p>
  
  <pre><code>## cell 1
  %matplotlib inline
  
  import numpy as np
  import matplotlib.pyplot as plt
  from mpl_toolkits.basemap import Basemap
  
  import pyproj
  from shapely.geometry import Point, Polygon, MultiPolygon
  from shapely.ops import transform as sh_transform
  from functools import partial
  
  wgs84_globe = pyproj.Proj(proj='latlong', ellps='WGS84')
  
  def disk_on_globe(lat, lon, radius):
      aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84',
                         lat_0=lat, lon_0=lon)
      return sh_transform(
          partial(pyproj.transform, aeqd, wgs84_globe),
          Point(0, 0).buffer(radius)
      )
  </code></pre>
  
  
  
  <pre><code>## cell 2
  def plot_poly_on_map(map_, pol):
      if isinstance(pol, Polygon):
          map_.plot(*(pol.exterior.xy), '-', latlon=True)
      else:
          assert isinstance(pol, MultiPolygon)
          for p in pol:
              map_.plot(*(p.exterior.xy), '-', latlon=True)
  
  plt.figure(figsize=(14, 12))
  map_ = Basemap(projection='cyl', resolution='c')
  map_.drawcoastlines(linewidth=0.25)
  
  for rad in range(1,10):
      plot_poly_on_map(
          map_,
          disk_on_globe(40.439, -79.976, rad * 1000 * 1000)
  )
  plt.show()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/VZmsn.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/VZmsn.png" alt="output of cell 2"></a></p>
  
  <pre><code>## cell 3
  def plot_poly_in_void(pol):
      if isinstance(pol, Polygon):
          plt.plot(*(pol.exterior.xy), '-')
      else:
          assert isinstance(pol, MultiPolygon)
          for p in pol:
              plt.plot(*(p.exterior.xy), '-', latlon=True)
  
  plt.figure()
  for rad in range(1,10):
      plot_poly_in_void(
          disk_on_globe(40.439, -79.976, rad * 1000 * 1000)
  )
  plt.show()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/9GBgk.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9GBgk.png" alt="output of cell 3"></a></p>
  
  <p>(The sunlit region shown at <a href="http://www.die.net/earth/rectangular.html" rel="nofollow noreferrer">http://www.die.net/earth/rectangular.html</a> is an  example of what a circle that crosses a pole <em>should</em> look like when projected onto an equirectangular map, as long as it's not an equinox today.)</p>