该演示程序(旨在在 IPython 笔记本中运行;您需要matplotlib
, mpl_toolkits.basemap
, pyproj
, and shapely
)应该在地球表面绘制越来越大的圆圈。只要圆不越过其中一根极点,它就能正常工作。如果发生这种情况,绘制在地图上时的结果完全是无稽之谈(见下面的单元格 2)
如果我将它们绘制在“虚空”而不是地图上(参见下面的单元格 3),则结果是正确的,因为如果您删除了从 +180 经度到 -180 经度的水平线,则曲线的其余部分将实际上划定了所需圆的内部和外部之间的边界。然而,他们是wrong因为多边形是无效的(.is_valid
是 False),更重要的是,多边形的非零绕数内部确实not圈出地图的正确区域。
我相信这一切正在发生,因为shapely.ops.transform
对经度 +180==-180 处的坐标奇点视而不见。这question是,如何检测问题并修复多边形,以便它确实包含地图的正确区域?在这种情况下,适当的修复方法是将 (X,+180) -- (X,-180) 中的水平线段替换为三行: (X,+180) -- (+90,+180) -- (+90,-180) -- (X,-180);但请注意,如果圆圈超出了south杆子上,固定线需要向南走。如果圆圈已经过去了both极点,我们将再次获得一个有效的多边形,但其内部将是其应有的补充。我需要检测所有这些情况并正确处理它们。另外,我不知道如何“编辑”形状优美的几何对象。
可下载的笔记本:https://gist.github.com/zackw/e48cb1580ff37acfee4d0a7b1d43a037 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/rectangle.html http://www.die.net/earth/rectangular.html这是一个穿过极点的圆的例子should看起来就像投影到等距柱状图上一样,只要今天不是春分点。)