Fix up shapely polygon object when discontinuous a

Question

This demo program (intended to be run in an IPython notebook; you need matplotlib, mpl_toolkits.basemap, pyproj, and shapely) 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)

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 wrong in that the polygon is invalid (.is_valid is False), and much more importantly, the nonzero-winding-number interior of the polygon does not enclose the correct region of the map.

I believe this is happening because shapely.ops.transform is blind to the coordinate singularity at +180==-180 longitude. The question 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 south pole, the fixup lines would need to go south instead. And if the circle had gone over both 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.

Downloadable notebook: https://gist.github.com/zackw/e48cb1580ff37acfee4d0a7b1d43a037

## 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)
    )

## 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()

## 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()

(The sunlit region shown at http://www.die.net/earth/rectangular.html is an example of what a circle that crosses a pole should look like when projected onto an equirectangular map, as long as it's not an equinox today.)

Answer 1

Manually fixing up the projected polygon turns out not to be that bad. There are two steps: first, find all segments of the polygon that cross the coordinate singularity at longitude ±180, and replace them with excursions to either the north or south pole, whichever is nearest; second, if the resulting polygon doesn't contain the origin point, invert it. Note that both steps must be carried out whether or not shapely thinks the projected polygon is "invalid"; depending on where the starting point is, it may cross one or both poles without being invalid.

This probably isn't the most efficient way to do it, but it works.

import pyproj
from shapely.geometry import Point, Polygon, box as Box
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):
    """Generate a shapely.Polygon object representing a disk on the
    surface of the Earth, containing all points within RADIUS meters
    of latitude/longitude LAT/LON."""

    aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84',
                       lat_0=lat, lon_0=lon)
    disk = sh_transform(
        partial(pyproj.transform, aeqd, wgs84_globe),
        Point(0, 0).buffer(radius)
    )

    # Fix up segments that cross the coordinate singularity at longitude ±180.
    # We do this unconditionally because it may or may not create a non-simple
    # polygon, depending on where the initial point was.
    boundary = np.array(disk.boundary)
    i = 0
    while i < boundary.shape[0] - 1:
        if abs(boundary[i+1,0] - boundary[i,0]) > 180:
            assert (boundary[i,1] > 0) == (boundary[i,1] > 0)
            vsign = -1 if boundary[i,1] < 0 else 1
            hsign = -1 if boundary[i,0] < 0 else 1
            boundary = np.insert(boundary, i+1, [
                [hsign*179, boundary[i,1]],
                [hsign*179, vsign*89],
                [-hsign*179, vsign*89],
                [-hsign*179, boundary[i+1,1]]
            ], axis=0)
            i += 5
        else:
            i += 1
    disk = Polygon(boundary)

    # If the fixed-up polygon doesn't contain the origin point, invert it.
    if not disk.contains(Point(lon, lat)):
        disk = Box(-180, -90, 180, 90).difference(disk)

    assert disk.is_valid
    assert disk.boundary.is_simple
    assert disk.contains(Point(lon, lat))
    return disk

The other problem -- mpl_toolkits.basemap.Basemap.plot producing garbage -- is not corrected by fixing up the polygon as above. However, if you manually project the polygon into map coordinates and then draw it using a descartes.PolygonPatch, that works, as long as the projection has a rectangular boundary, and that's enough of a workaround for me. (I think it would work for any projection if one added a lot of extra points along all straight lines at the map boundary.)

%matplotlib inline
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
from descartes import PolygonPatch

plt.figure(figsize=(14, 12))
map_ = Basemap(projection='cea', resolution='c')
map_.drawcoastlines(linewidth=0.25)

for rad in range(3,19,2):
    plt.gca().add_patch(PolygonPatch(
        sh_transform(map_,
            disk_on_globe(40.439, -79.976, rad * 1000 * 1000)),
        alpha=0.1))    
plt.show()