获取重叠的纬度经度域的包络(或边界)
Get envelope (or boundary) of overlapping latitude longitude domains
我有多个重叠的经度、纬度网格域。现在,我需要围绕这些域绘制 "outer line",但我找不到合适的解决方案。所有这些都需要在 python 中完成,我在那里进行数据处理和绘制字段(来自不同网格域的值的组合)。
我知道这可以通过 postgresql+postgis 以某种方式与域边界相交来实现,但我不想仅仅为了这个问题而涉及数据库。
凸包不是我要找的,因为这不会涵盖 "the concave part of the problem"。
在我看来,Alpha 形状也不是正确的方法,因为它们没有提供明确定义的解决方案。显然这完全取决于所选的 alpha 值。找到的例子,例如这里:http://blog.thehumangeo.com/2014/05/12/drawing-boundaries-in-python/
我认为这个问题必须有一个明确定义的单一解决方案,因为我只是在寻找通过 lon/lat 点定义的一些线的交点。
现在,为了更清楚,这是一个人工示例,我只创建两个域,加入它们以获得所有(两个)域的 lon/lat 点数组,并尝试一些事情匀称(因为这看起来最有希望):
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()
下图显示了重叠的两个示例域。我现在对两个域的一组点构建 "the outer boundary" 感兴趣(即,如果一个人沿着域边界行走(顺时针)并且总是在交叉路口左转)。
#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()
凸包缺少 "inner" 部分,不是所需的解决方案:
#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 相关的包可以在这里提供帮助?
在我看来,你所追求的是两个多边形的并集的外边界:
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()
我有多个重叠的经度、纬度网格域。现在,我需要围绕这些域绘制 "outer line",但我找不到合适的解决方案。所有这些都需要在 python 中完成,我在那里进行数据处理和绘制字段(来自不同网格域的值的组合)。 我知道这可以通过 postgresql+postgis 以某种方式与域边界相交来实现,但我不想仅仅为了这个问题而涉及数据库。
凸包不是我要找的,因为这不会涵盖 "the concave part of the problem"。 在我看来,Alpha 形状也不是正确的方法,因为它们没有提供明确定义的解决方案。显然这完全取决于所选的 alpha 值。找到的例子,例如这里:http://blog.thehumangeo.com/2014/05/12/drawing-boundaries-in-python/
我认为这个问题必须有一个明确定义的单一解决方案,因为我只是在寻找通过 lon/lat 点定义的一些线的交点。
现在,为了更清楚,这是一个人工示例,我只创建两个域,加入它们以获得所有(两个)域的 lon/lat 点数组,并尝试一些事情匀称(因为这看起来最有希望):
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()
下图显示了重叠的两个示例域。我现在对两个域的一组点构建 "the outer boundary" 感兴趣(即,如果一个人沿着域边界行走(顺时针)并且总是在交叉路口左转)。
#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()
凸包缺少 "inner" 部分,不是所需的解决方案:
#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 相关的包可以在这里提供帮助?
在我看来,你所追求的是两个多边形的并集的外边界:
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()