Fix up formschön Polygon-Objekt, wenn diskontinuierlich nach Kartenprojektion

Question

Dieses Demo-Programm (bestimmt in einem IPython Notebook ausgeführt werden, Sie müssen matplotlib, mpl_toolkits.basemap, pyproj und shapely) soll immer größere Kreise auf der Oberfläche der Erde plotten . Es funktioniert korrekt, solange der Kreis einen der Pole nicht überquert. Wenn das passiert, ist das Ergebnis völliger Unsinn, wenn es auf einer Karte geplottet wird (siehe unten Zelle 2)

Wenn ich sie "in einer Lücke" anstatt auf einer Karte (siehe unten Zelle 3) plotten, sind die Ergebnisse in der Wenn Sie die horizontale Linie von +180 auf -180 Länge entfernen, würde der Rest der Kurve tatsächlich die Grenze zwischen dem Inneren und Äußeren des gewünschten Kreises begrenzen. Allerdings sind sie falsch in dem das Polygon ist ungültig (.is_valid ist False), und viel wichtiger, die nicht-Null-Windungszahl Innenraum des Polygons tut nicht umschließen die richtige Region der Karte.

Ich glaube, das passiert, weil shapely.ops.transform ist blind für die Koordinaten Singularität bei +180 == - 180 Länge. Die Frage ist, wie erkenne ich das Problem und das Polygon reparieren, so dass es den richtigen Bereich der Karte umschließt? In diesem Fall wäre eine angemessene Korrektur das Ersetzen des horizontalen Segments von (X, + 180) - (X, -180) durch drei Zeilen (X, + 180) - (+ 90, + 180) - (+ 90, -180) - (X, -180); aber beachte, dass, wenn der Kreis über die Süd Stange gegangen wäre, die Reparaturlinien stattdessen nach Süden gehen müssten. Und wenn der Kreis über beide Pole gegangen wäre, hätten wir wieder ein gültiges Polygon, aber sein Inneres wäre die Ergänzung dessen, was es sein sollte. Ich muss alle diese Fälle erkennen und sie richtig behandeln. Ich weiß auch nicht, wie man ein formschönes Geometrieobjekt "bearbeitet".

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

(Die sonnigen Region beigezeigtist ein Beispiel dafür, was ein Kreis, der eine Stange kreuzt sollte aussehen, wenn auf eine equirectangular Karte projiziert wird, solange es heute kein Equinox.)

Answer 1

manuell das projizierte Polygon bis zur Festsetzung erweist sich als nicht zu sein dass schlecht. Es gibt zwei Schritte: Zuerst finden Sie alle Segmente des Polygons, die die coordinate singularity bei Länge ± 180 überqueren, und ersetzen Sie sie durch Exkursionen entweder zum Nord- oder zum Südpol, was auch immer am nächsten ist; Zweitens, wenn das resultierende Polygon den Ursprungspunkt nicht enthält, invertieren Sie es. Beachten Sie, dass beide Schritte ausgeführt werden müssen ob formschön denkt das projizierte Polygon ist "ungültig"; Je nachdem, wo der Startpunkt liegt, kann er einen oder beide Pole ohne überschreiten.

Dies ist wahrscheinlich nicht der effizienteste Weg, es zu tun, aber es funktioniert.

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 

Das andere Problem - mpl_toolkits.basemap.Basemap.plot Herstellung Müll - ist nicht korrigiert durch das Polygon, wie oben Fixierung auf. Wenn Sie das Polygon jedoch manuell in Kartenkoordinaten projizieren und es dann mit einer descartes.PolygonPatch zeichnen, funktioniert das, solange die Projektion eine rechteckige Grenze hat, und das reicht für mich aus. (I denken es wäre für jede Projektion arbeiten, wenn man eine Menge Extra-Punkte entlang aller geraden Linien auf der Karte Grenze hinzugefügt.)

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