这个演示程序(打算在IPython笔记本上运行;您需要matplotlib
、mpl_toolkits.basemap
、pyproj
、和{
如果我把它们画成“空白”而不是在地图上(见下面的单元格3),结果是正确的,因为如果你删除了从+180到-180经度的水平线,曲线的其余部分确实会划定所需圆的内部和外部之间的边界。然而,它们是错误的,因为多边形是无效的(.is_valid
是假的),更重要的是,多边形内部的非零缠绕数并没有包围地图的正确区域。在
我相信这是因为shapely.ops.transform
对经度+180==-180的坐标奇点是视而不见的。问题是,如何检测问题并修复多边形,使其包含地图的正确区域?在这种情况下,一个合适的固定方法是用三条线(X,+180)--(X,-180)替换水平段,(X,+180)--(+90,-180)--(X,-180);但请注意,如果圆越过南极,则修正线需要向南。如果圆经过两个极点,我们就会有一个有效的多边形,但它的内部将是它应该是什么的补充。我需要发现所有这些案件并正确处理。此外,我不知道如何“编辑”一个形状几何体对象。在
可下载笔记本: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)
)
^{pr2}$
## 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是否认为投影多边形“无效”,都必须执行这两个步骤;根据起点的位置,它可以跨越一个或两个极点,而不是无效。在
这可能不是最有效的方法,但它是有效的。在
另一个产生垃圾的问题是通过修正多边形来纠正的。但是,如果您手动将多边形投影到地图坐标系中,然后使用
^{pr2}$descartes.PolygonPatch
绘制它,只要投影有一个矩形边界,这对我来说就足够了。(我认为,如果在地图边界的所有直线上添加许多额外的点,那么对于任何投影都是有效的。)相关问题 更多 >
编程相关推荐