在地图投影后不连续时修复形状多边形对象
Fix up shapely polygon object when discontinuous after map projection
此演示程序（打算在 IPython 笔记本中 运行；您需要 matplotlib、mpl_toolkits.basemap、pyproj 和 shapely) 应该在地球表面绘制越来越大的圆圈。只要圆圈不越过其中一个极点，它就可以正常工作。如果发生这种情况，在地图上绘制时结果完全是胡说八道（见下面的单元格 2）
如果我绘制它们 "in a void" 而不是在地图上（见下面的单元格 3），结果是正确的，如果你删除了从 +180 到 -180 经度的水平线，曲线的其余部分确实会划定所需圆的内部和外部之间的边界。然而，它们是 错误的 因为多边形是无效的（.is_valid 是错误的），更重要的是，多边形的非零缠绕数内部确实 不包围地图的正确区域。
我相信这是因为 shapely.ops.transform 对 +180==-180 经度处的坐标奇点视而不见。 问题 是，如何检测问题并修复多边形，以便它包围地图的正确区域？在这种情况下，适当的修正是用三行替换 (X,+180) -- (X,-180) 的水平线段，(X,+180) -- (+90,+180) -- (+90,-180) -- (X,-180);但请注意，如果圆圈已经越过 南 极，则修正线需要向南移动。如果圆已经越过了 两个 极点，我们将再次得到一个有效的多边形，但它的内部将是它应该是的补充。我需要检测所有这些情况并正确处理它们。另外，我不知道如何 "edit" 形状优美的几何对象。
可下载笔记本：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()
（在 http://www.die.net/earth/rectangular.html 处显示的阳光照射区域是一个例子，当投影到等距柱状地图上时，穿过极点的圆 应该 看起来像，只要今天不是春分。)
手动修复投影多边形结果 并没有那么糟糕。
有两个步骤：首先，找到在经度 ±180 处穿过 coordinate singularity 的所有多边形线段，并将它们替换为到北极或南极的短途旅行，以最近的为准；其次，如果生成的多边形不包含原点，则将其反转。请注意，无论是否 shapely 认为投影的多边形是 "invalid"，这两个步骤都必须进行 ；根据起点的位置，它可能会越过一个或两个极点 而不会 无效。
这可能不是最有效的方法，但它确实有效。
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
另一个问题——mpl_toolkits.basemap.Basemap.plot 产生垃圾——没有通过修复上面的多边形来纠正。但是，如果您手动将多边形投影到地图坐标中，然后使用 descartes.PolygonPatch 绘制它，只要投影具有矩形边界，这就可以了，这对我来说已经足够了。 （我 认为 如果沿着地图边界的所有直线添加很多额外的点，它对任何投影都有效。）
%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()
此演示程序（打算在 IPython 笔记本中 运行；您需要 matplotlib、mpl_toolkits.basemap、pyproj 和 shapely) 应该在地球表面绘制越来越大的圆圈。只要圆圈不越过其中一个极点，它就可以正常工作。如果发生这种情况，在地图上绘制时结果完全是胡说八道（见下面的单元格 2）
如果我绘制它们 "in a void" 而不是在地图上（见下面的单元格 3），结果是正确的，如果你删除了从 +180 到 -180 经度的水平线，曲线的其余部分确实会划定所需圆的内部和外部之间的边界。然而，它们是 错误的 因为多边形是无效的（.is_valid 是错误的），更重要的是，多边形的非零缠绕数内部确实 不包围地图的正确区域。
我相信这是因为 shapely.ops.transform 对 +180==-180 经度处的坐标奇点视而不见。 问题 是，如何检测问题并修复多边形，以便它包围地图的正确区域？在这种情况下，适当的修正是用三行替换 (X,+180) -- (X,-180) 的水平线段，(X,+180) -- (+90,+180) -- (+90,-180) -- (X,-180);但请注意，如果圆圈已经越过 南 极，则修正线需要向南移动。如果圆已经越过了 两个 极点，我们将再次得到一个有效的多边形，但它的内部将是它应该是的补充。我需要检测所有这些情况并正确处理它们。另外，我不知道如何 "edit" 形状优美的几何对象。
可下载笔记本：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()
（在 http://www.die.net/earth/rectangular.html 处显示的阳光照射区域是一个例子，当投影到等距柱状地图上时，穿过极点的圆 应该 看起来像，只要今天不是春分。)
手动修复投影多边形结果 并没有那么糟糕。 有两个步骤：首先，找到在经度 ±180 处穿过 coordinate singularity 的所有多边形线段，并将它们替换为到北极或南极的短途旅行，以最近的为准；其次，如果生成的多边形不包含原点，则将其反转。请注意，无论是否 shapely 认为投影的多边形是 "invalid"，这两个步骤都必须进行 ；根据起点的位置，它可能会越过一个或两个极点 而不会 无效。
这可能不是最有效的方法，但它确实有效。
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
另一个问题——mpl_toolkits.basemap.Basemap.plot 产生垃圾——没有通过修复上面的多边形来纠正。但是，如果您手动将多边形投影到地图坐标中，然后使用 descartes.PolygonPatch 绘制它，只要投影具有矩形边界，这就可以了，这对我来说已经足够了。 （我 认为 如果沿着地图边界的所有直线添加很多额外的点，它对任何投影都有效。）
%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()