问题: 对于点位置,用Python计算逆距离加权(IDW)插值的最佳方法是什么?
背景知识:
目前,我正在使用RPy2与R及其gstat模块进行接口。不幸的是,gstat模块与arcgisscripting冲突,我通过在单独的过程中运行基于RPy2的分析来解决这个问题。即使在最新的/将来的版本中解决了此问题,并且可以提高效率,我仍然希望消除对安装R的依赖。
gstat网站确实提供了一个独立的可执行文件,该文件更容易与我的python脚本一起打包,但是我仍然希望有一个Python解决方案,该解决方案不需要多次写入磁盘并启动外部进程。在执行的处理中,分别由点和值组成的对插值函数的调用次数可能接近20,000。
我特别需要对点进行插值,因此就性能而言,使用ArcGIS中的IDW函数生成栅格声音比使用R还要差,除非有一种方法可以有效地仅掩盖我需要的点。即使进行了此修改,我也不希望性能如此出色。我将研究此选项作为另一种选择。更新:这里的问题是您与所使用的像元大小有关。如果减小像元大小以获得更好的精度,则处理将花费很长时间。如果需要特定点的值,还需要通过按点提取.....来进行后续处理。
我看过scipy文档,但是看起来没有直接计算IDW的方法。
我正在考虑滚动自己的实现,可能会使用某些scipy功能来定位最接近的点并计算距离。
我是否缺少明显的东西?有没有我没有看到过的python模块,正是我想要的吗?在scipy的帮助下创建自己的实现是明智的选择吗?
10月20日更改:此类Invdisttree结合了反距离权重和
scipy.spatial.KDTree。
忘记原始的蛮力答案;这是分散数据插值的首选方法。
""" invdisttree.py: inverse-distance-weighted interpolation using KDTree
fast, solid, local
"""
from __future__ import division
import numpy as np
from scipy.spatial import cKDTree as KDTree
# http://docs.scipy.org/doc/scipy/reference/spatial.html
__date__ = "2010-11-09 Nov" # weights, doc
#...............................................................................
class Invdisttree:
""" inverse-distance-weighted interpolation using KDTree:
invdisttree = Invdisttree( X, z ) -- data points, values
interpol = invdisttree( q, nnear=3, eps=0, p=1, weights=None, stat=0 )
interpolates z from the 3 points nearest each query point q;
For example, interpol[ a query point q ]
finds the 3 data points nearest q, at distances d1 d2 d3
and returns the IDW average of the values z1 z2 z3
(z1/d1 + z2/d2 + z3/d3)
/ (1/d1 + 1/d2 + 1/d3)
= .55 z1 + .27 z2 + .18 z3 for distances 1 2 3
q may be one point, or a batch of points.
eps: approximate nearest, dist <= (1 + eps) * true nearest
p: use 1 / distance**p
weights: optional multipliers for 1 / distance**p, of the same shape as q
stat: accumulate wsum, wn for average weights
How many nearest neighbors should one take ?
a) start with 8 11 14 .. 28 in 2d 3d 4d .. 10d; see Wendel's formula
b) make 3 runs with nnear= e.g. 6 8 10, and look at the results --
|interpol 6 - interpol 8| etc., or |f - interpol*| if you have f(q).
I find that runtimes don't increase much at all with nnear -- ymmv.
p=1, p=2 ?
p=2 weights nearer points more, farther points less.
In 2d, the circles around query points have areas ~ distance**2,
so p=2 is inverse-area weighting. For example,
(z1/area1 + z2/area2 + z3/area3)
/ (1/area1 + 1/area2 + 1/area3)
= .74 z1 + .18 z2 + .08 z3 for distances 1 2 3
Similarly, in 3d, p=3 is inverse-volume weighting.
Scaling:
if different X coordinates measure different things, Euclidean distance
can be way off. For example, if X0 is in the range 0 to 1
but X1 0 to 1000, the X1 distances will swamp X0;
rescale the data, i.e. make X0.std() ~= X1.std() .
A nice property of IDW is that it's scale-free around query points:
if I have values z1 z2 z3 from 3 points at distances d1 d2 d3,
the IDW average
(z1/d1 + z2/d2 + z3/d3)
/ (1/d1 + 1/d2 + 1/d3)
is the same for distances 1 2 3, or 10 20 30 -- only the ratios matter.
In contrast, the commonly-used Gaussian kernel exp( - (distance/h)**2 )
is exceedingly sensitive to distance and to h.
"""
# anykernel( dj / av dj ) is also scale-free
# error analysis, |f(x) - idw(x)| ? todo: regular grid, nnear ndim+1, 2*ndim
def __init__( self, X, z, leafsize=10, stat=0 ):
assert len(X) == len(z), "len(X) %d != len(z) %d" % (len(X), len(z))
self.tree = KDTree( X, leafsize=leafsize ) # build the tree
self.z = z
self.stat = stat
self.wn = 0
self.wsum = None;
def __call__( self, q, nnear=6, eps=0, p=1, weights=None ):
# nnear nearest neighbours of each query point --
q = np.asarray(q)
qdim = q.ndim
if qdim == 1:
q = np.array([q])
if self.wsum is None:
self.wsum = np.zeros(nnear)
self.distances, self.ix = self.tree.query( q, k=nnear, eps=eps )
interpol = np.zeros( (len(self.distances),) + np.shape(self.z[0]) )
jinterpol = 0
for dist, ix in zip( self.distances, self.ix ):
if nnear == 1:
wz = self.z[ix]
elif dist[0] < 1e-10:
wz = self.z[ix[0]]
else: # weight z s by 1/dist --
w = 1 / dist**p
if weights is not None:
w *= weights[ix] # >= 0
w /= np.sum(w)
wz = np.dot( w, self.z[ix] )
if self.stat:
self.wn += 1
self.wsum += w
interpol[jinterpol] = wz
jinterpol += 1
return interpol if qdim > 1 else interpol[0]
#...............................................................................
if __name__ == "__main__":
import sys
N = 10000
Ndim = 2
Nask = N # N Nask 1e5: 24 sec 2d, 27 sec 3d on mac g4 ppc
Nnear = 8 # 8 2d, 11 3d => 5 % chance one-sided -- Wendel, mathoverflow.com
leafsize = 10
eps = .1 # approximate nearest, dist <= (1 + eps) * true nearest
p = 1 # weights ~ 1 / distance**p
cycle = .25
seed = 1
exec "\n".join( sys.argv[1:] ) # python this.py N= ...
np.random.seed(seed )
np.set_printoptions( 3, threshold=100, suppress=True ) # .3f
print "\nInvdisttree: N %d Ndim %d Nask %d Nnear %d leafsize %d eps %.2g p %.2g" % (
N, Ndim, Nask, Nnear, leafsize, eps, p)
def terrain(x):
""" ~ rolling hills """
return np.sin( (2*np.pi / cycle) * np.mean( x, axis=-1 ))
known = np.random.uniform( size=(N,Ndim) ) ** .5 # 1/(p+1): density x^p
z = terrain( known )
ask = np.random.uniform( size=(Nask,Ndim) )
#...............................................................................
invdisttree = Invdisttree( known, z, leafsize=leafsize, stat=1 )
interpol = invdisttree( ask, nnear=Nnear, eps=eps, p=p )
print "average distances to nearest points: %s" % \
np.mean( invdisttree.distances, axis=0 )
print "average weights: %s" % (invdisttree.wsum / invdisttree.wn)
# see Wikipedia Zipf's law
err = np.abs( terrain(ask) - interpol )
print "average |terrain() - interpolated|: %.2g" % np.mean(err)
# print "interpolate a single point: %.2g" % \
# invdisttree( known[0], nnear=Nnear, eps=eps )
问题内容: 我需要有效地计算给定数组中每个点到另一个数组中每个其他点的欧几里德 加权 距离。这是我拥有的代码,可以正常工作: 最终列表包含的子列表数量与第一个数组中的点数量相同,每个元素中的元素数量与第二个数组中的点数量相同(加权距离)。 我想知道是否有一种方法可以提高此代码的效率,也许使用cdist函数,因为当数组中有很多元素(在我的情况下,它们具有)并且必须检查时,此过程会变得非常昂贵许多数组
问题内容: 在Python + Sqlite中是否有可用的字符串相似性度量,例如与模块有关? 用例示例: 此查询应匹配ID为1的行,但不匹配ID为2的行: 如何在Sqlite + Python中做到这一点? 关于我到目前为止发现的注释: 该Levenshtein距离,即单字符编辑(插入,删除或替换)的最小数量需要改变一个字到另一个,可能是有用的,但我不知道是否SQLite中存在的正式实施(我看到一
问题内容: 我正在寻找Python中的地球移动器的距离(或快速EMD)实现。关于在哪里找到它的任何线索,我在网上已经看够了。我想在我正在做的图像检索项目中使用它。谢谢。 编辑:我发现了使用纸浆库的一个非常好的解决方案。该页面还包含设置所需的说明。 问题答案: OpenCv for Python中有一个出色的实现。该函数的名称为CalcEMD2,用于比较两个图像直方图的简单代码如下所示: 我使用Py
我需要计算汽车行驶的距离!不是距离,不是距离到否。如果我们通过谷歌提供的API计算,距离可以完全不同。谷歌可以提供从一个点到另一个点的1公里距离,但汽车可以按照骑手想要的方式行驶800米。使用加速计没有帮助。它适用于步行,但绝不适用于更快的速度。 我尝试过使用Google的位置API:距离到或距离之间根本不是一个选项。它可以给出与IN REAL截然不同的结果。在真实的汽车中,可以通过非常短的地方并
本文向大家介绍编辑距离相关面试题,主要包含被问及编辑距离时的应答技巧和注意事项,需要的朋友参考一下 参考回答: 概念 编辑距离的作用主要是用来比较两个字符串的相似度的 编辑距离,又称Levenshtein距离(莱文斯坦距离也叫做Edit Distance),是指两个字串之间,由一个转成另一个所需的最少编辑操作次数,如果它们的距离越大,说明它们越是不同。许可的编辑操作包括将一个字符替换成另一个字符,
问题内容: 我有一个一维数字数组,想计算所有成对的欧几里得距离。我有一种方法(感谢SO)在广播中执行此操作,但是它效率低下,因为它两次计算每个距离。而且它的伸缩性不好。 这是一个示例,它为我提供了1000个数字的数组。 我可以使用scipy / numpy / scikit-learn中最快的实现来执行此操作,因为它必须扩展到一维数组具有> 10k值的情况。 注意:矩阵是对称的,所以我猜想通过解决