python - 获取重叠纬度经度域的包络线(或边界)

Question

我有许多重叠的经度、纬度网格域。现在，我需要围绕这些域绘制“外线”，但我找不到合适的解决方案。所有这些都需要在 python 中完成，我正在其中进行数据处理和字段绘图(来自不同网格域的值的组合)。我知道这可以通过 postgresql+postgis 通过某种方式与域边界相交来实现，但我不想仅仅为了这个问题而涉及数据库。

凸包不是我想要的，因为这不会覆盖“问题的凹部分”。在我看来，Alpha 形状也不是正确的方法，因为它们没有提供明确定义的解决方案。显然这一切都取决于所选的 alpha 值。发现的例子例如这里:http://blog.thehumangeo.com/2014/05/12/drawing-boundaries-in-python/

我认为这个问题必须有一个明确定义的单一解决方案，因为我只是在寻找通过经/纬度点定义的一些线的交点。

现在，为了更清楚起见，这里有一个人工示例，我只创建两个域，将它们连接起来以获得所有(两个)域的经/纬度点数组，并尝试使用 shapely 进行一些操作(如下所示)似乎最有希望):

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import MultiPoint, Polygon

#two sample domains of longitudes and latitudes
dom1 = np.array([np.linspace(12.0,19.0,20), np.linspace(41.0,42.0,20)])
dom1 = np.append(dom1, [np.linspace(19.0,19.5,20), np.linspace(42.0,32.0,20)],axis=1)
dom1 = np.append(dom1, [np.linspace(19.5,11.5,20), np.linspace(32.0,31.0,20)],axis=1)
dom1 = np.append(dom1, [np.linspace(11.5,12.0,20), np.linspace(31.0,41.0,20)],axis=1)

dom2 = np.array([np.linspace(14.0,21.0,20), np.linspace(44.0,44.5,20)])
dom2 = np.append(dom2, [np.linspace(21.0,19.5,20), np.linspace(44.5,35.0,20)],axis=1)
dom2 = np.append(dom2, [np.linspace(19.5,14.,20), np.linspace(35.0,35.2,20)],axis=1)
dom2 = np.append(dom2, [np.linspace(14.,14.,20), np.linspace(35.2,44.0,20)],axis=1)

plt.plot(dom1[0,:],dom1[1,:],'k')
plt.plot(dom2[0,:],dom2[1,:],'r')         
plt.show()

下图显示了两个示例域重叠。我现在对构建两个域的“外边界”的一组点感兴趣(即，如果一个人沿着域边界行走(顺时针)并始终在交叉点左转)。

#join the domains in one array
doms = np.hstack([dom1,dom2])

#convex hull
points = MultiPoint(zip(doms[0,:],doms[1,:]))
x,y = points.convex_hull.exterior.xy

plt.plot(doms[0,:],doms[1,:],'k.')
plt.plot(x,y,'ro-')
plt.show()

凸包缺少“内部”部分，不是所需的解决方案:

#polygon boundary
poly = Polygon(zip(doms[0,:],doms[1,:]))
x1,y1 = poly.boundary.xy
plt.plot(doms[0,:],doms[1,:],'k.')
plt.plot(x1,y1,'go-')
plt.show()

对多边形边界进行实验也不会产生任何有用的结果:

可以用shapely来解决这个问题吗？或者还有其他与 gis 相关的软件包可以在这里提供帮助吗？

Answer 1

在我看来，您所追求的是两个多边形并集的外部边界:

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import MultiPoint, Polygon

#two sample domains of longitudes and latitudes
dom1 = np.array([np.linspace(12.0,19.0,20), np.linspace(41.0,42.0,20)])
dom1 = np.append(dom1, [np.linspace(19.0,19.5,20), np.linspace(42.0,32.0,20)],axis=1)
dom1 = np.append(dom1, [np.linspace(19.5,11.5,20), np.linspace(32.0,31.0,20)],axis=1)
dom1 = np.append(dom1, [np.linspace(11.5,12.0,20), np.linspace(31.0,41.0,20)],axis=1)

dom2 = np.array([np.linspace(14.0,21.0,20), np.linspace(44.0,44.5,20)])
dom2 = np.append(dom2, [np.linspace(21.0,19.5,20), np.linspace(44.5,35.0,20)],axis=1)
dom2 = np.append(dom2, [np.linspace(19.5,14.,20), np.linspace(35.0,35.2,20)],axis=1)
dom2 = np.append(dom2, [np.linspace(14.,14.,20), np.linspace(35.2,44.0,20)],axis=1)

P = Polygon(dom1.T)
Q = Polygon(dom2.T)
P = P.union(Q)

dom3 = np.asarray(P.exterior.coords).T

plt.plot(dom1[0,:],dom1[1,:],'k')
plt.plot(dom2[0,:],dom2[1,:],'r')         
plt.plot(dom3[0,:],dom3[1,:],'go-')         
plt.show()