接上期基于 Python 的自然邻域法空间插值的实现与思考。
上期说到,我们仅仅利用自然邻域法基础原理进行插值,会出现许多空值、异常值,且与ArcGIS相同分辨率、范围下的插值结果对比(对比图如下),结果较差。主要体现在:
- 插值结果范围内有空值,而ArcGIS没有,可能是ArcGIS做了其他的一些处理。
- ArcGIS插值结果仅包含了最外层点组成的面内的数据,显然,边界外的数据插值结果异常值较多。
- 部分区域插值结果较差(例如下图左,左下角),仍有需要改进的地方。
后续优化实践过程中遇到了瓶颈,偶然机会发现可能是泰森多边形构建的问题,主要问题在于:泰森多边形存在无限区域,导致后续加权计算时面积异常。无限区域因为缺少顶点因此无法形成多边形! 解决的方法也比较容易理解,就是有效化泰森多边形无限区域。经过尝试,基本解决上述问题。达到了函数构建的预期效果。主要解决的内容如下:
有效化泰森多边形无限区域
有效化方法如下,有需要可以学习参考。
def VoronoiVertices(Points):
"""
生成泰森多边形标记(区域序号,多边形点端点位置,端点数组,无限区域序号),并使无限区域有限。
"""
VOR = Voronoi(Points, incremental = True)
Vertices = VOR.vertices
NewVertices = Vertices.copy()
# 几何中心和半径
Center = Points.mean(axis=0)
Radius = Points.ptp().max() * 2
# 构建包含给定点的所有脊
AllRidges = {}
RidgePoints = VOR.ridge_points
RidgeVertices = VOR.ridge_vertices
for p1, p2, v1, v2 in np.concatenate([RidgePoints, RidgeVertices], axis = 1):
AllRidges.setdefault(p1, []).append((p2,v1,v2))
AllRidges.setdefault(p2, []).append((p1,v1,v2))
# 重构无限区域
PointRegion = VOR.point_region
Regions = VOR.regions[:]
# 记录无限区域,方便提取边界
InfiniteRegion = []
for p1, pr in enumerate(PointRegion):
Vertice = np.array(Regions[pr])
# 跳过正常区域
if (Vertice != -1).all() and Vertice.size != 0:
continue
Ridges = AllRidges[p1]
NewRegion = Vertice[Vertice != -1]
for p2, v1, v2 in Ridges:
if v2 0:
v1, v2 = v2, v1
if v1 >= 0:
# 有限脊:已在该区域
continue
# 计算无限脊的缺失端点
T = Points[p2] - Points[p1] # 切线
T /= np.linalg.norm(T)
n = np.array([-T[1], T[0]]) # 正常
MidPoint = Points[[p1, p2]].mean(axis=0)
Direction = np.sign(np.dot(MidPoint - Center, n)) * n
FarPoint = Vertices[v2] + Direction * Radius
NewRegion = np.append(NewRegion, len(NewVertices))
NewVertices = np.append(NewVertices, [FarPoint], axis = 0)
# 更新
Regions[pr] = list(NewRegion)
# 记录无限区域点
InfiniteRegion.append(p1)
return PointRegion, Regions, NewVertices, InfiniteRegion
新的泰森多边形
更新后的插值结果对比
总结
插值结果用插值点组成的最大外部多边形裁剪之后与ArcGIS结果几乎一致,满足精度需求。
算法缺点:插值速度较慢。
声明:欢迎转载、转发本号原创内容,可留言区留言或者后台联系小编(微信:gavin7675)进行授权。气象学家公众号转载信息旨在传播交流,其内容由作者负责,不代表本号观点。文中部分图片来源于网络,如涉及作品内容、版权和其他问题,请后台联系小编处理。
欢迎加入气象学家交流群
请备注:
姓名/昵称-单位/学校-研究方向