为数据分析匹配两个网格,有什么好的算法吗?

2 投票
5 回答
3489 浏览
提问于 2025-04-16 19:19

我想在Python中比较两个间隔不同的数据集。我总是想找到最接近的匹配(最近邻),并重新调整数据,看看这个例子:

数据集A:

ALTITUDE[m]   VALUE
1.            a1
2.            a2
3.            a3
4.            a4

数据集B:

ALTITUDE[m]   VALUE
0.7           b1
0.9           b2
1.7           b3
2.            b4
2.4           b5
2.9           b6
3.1           b7
3.2           b8
3.9           b9
4.1           b10

aibi包含双精度数字,但也有nan字段。

我想把数据集B转换到数据集A的高度网格上,但由于数据集A的高度级别比数据集B少,我想对它们进行平均。

ALTITUDE[m]   VALUE
1.            median(b1,b2)
2.            median(b3,b4,b5)
3.            median(b6,b7,b8)
4.            median(b9,b10)

也就是说,找到了最接近的高度级别并进行了平均。

反过来,如果我想把数据集A匹配到数据集B的网格上,数据集A应该看起来像这样(最近邻):

ALTITUDE[m]   VALUE
0.7           a1
0.9           a1
1.7           a2
2.            a2
2.4           a2
2.9           a3
3.1           a3
3.2           a3
3.9           a4
4.1           a4

也许这个问题甚至有个名字(我想这应该是个常见问题),但我不知道,所以无法搜索。我相信有一种高效的方法来解决这个问题,除了自己编写明显的解决方案(但我担心这样效率不高,还可能引入很多bug)。

最好使用numpy。

编辑:感谢四位贡献者的意见。我学到了一些东西,对我没有问得很清楚表示歉意。我自己也在理解这个问题的过程中。你们的回答让我注意到了interp1d的用法,而这个回答让我能够利用它。我会很快发布结果。我只能接受一个答案,但任何答案都可以。

5 个回答

1

这里有一种方法:

import numpy as np

def regrid_op(x, y, xi, op=np.median):
    x, y, xi = np.atleast_1d(x, y, xi)
    if (x.ndim, y.ndim, xi.ndim) != (1, 1, 1):
        raise ValueError("works only for 1D data")

    yi = np.zeros(xi.shape, dtype=y.dtype)
    yi.fill(np.nan)

    # sort data
    j = np.argsort(x)
    x = x[j]
    y = y[j]

    # group items by nearest neighbour
    x0s = np.r_[xi, np.inf]
    xc = .5*(x0s[:-1] + x0s[1:])

    j0 = 0
    for i, j1 in enumerate(np.searchsorted(x, xc)):
        print "x =", xi[i], ", y =", y[j0:j1] # print some debug info
        yi[i] = op(y[j0:j1])
        j0 = j1

    return yi

xi = np.array([1, 2, 3, 4])
x = np.array([0.7, 0.9, 1.7, 2.0, 2.4, 2.9, 3.1, 3.2, 3.9, 4.1])
y = np.array([1,   2,   3,   4,   5,   6,   7,   8,   9,   10.])

print regrid_op(x, y, xi)

我没有找到办法可以让对xi数组中每个项目的循环变得更快,所以只要网格A中的点数不是太多,这种方法应该是高效的。

补充说明:这也假设xi中的点是排好序的。

2

有两个前提假设:

1:你不是在找最近的邻居,而是想找某个范围内的所有高度。比如说,对于a1,你想找所有在a1上下0.5范围内的bn(根据你的例子,这样会得到b1和b2)。我会把“最近的邻居”定义得更严格一些。

2:你在计算中位数时不算nan(缺失值),因为numpy会把它们当作无穷大,这在某些IEEE标准下是这样,但我觉得这有点奇怪。所以我们可以用scipy.stats里的nanmedian来处理这个问题。

我会这样做:

from numpy import *
from pylab import *

A_Alt = array([1,2,3,4])
A_Val = array([.33, .5, .6, .8])
B_Alt = array([.7, 0.9, 1.7, 2., 2.4, 2.9, 3.1, 3.2, 3.9, 4.1])
B_Val = array([.3, NaN, .8, .6, .7, .4, .3, NaN, .99, 1.3])

range = .5

B_Agrid = [nanmedian(B_Val[abs(B_Alt - k)<range]).item() for k in A_Alt]
A_Bgrid = [nanmedian(A_Val[abs(A_Alt - k)<range]).item() for k in B_Alt]    

我们找到所有B_Alt到A_Alt中k的距离小于指定范围的索引。然后我们取这些B_Val的中位数。对于A_Bgrid也是一样,结果会如你所要求的那样。

==编辑==

关于你提到的最近邻居的假设有所不同:我们把最近邻居定义为绝对高度差最小的条目(如果有并列的情况就取多个),同时这些值不能是nan。注意,这样的结果和你的例子不一致,因为b1并不是a1的最近邻居,因为b2更近。

在这个假设下,以下代码应该可以工作:

from numpy import *
from pylab import *
from scipy.stats import nanmedian

A_Alt = array([1,2,3,4])
A_Val = array([.33, .5, .6, .8])
B_Alt = array([.7, 0.9, 1.7, 2., 2.4, 2.9, 3.1, 3.2, 3.9, 4.1])
B_Val = array([.3, NaN, .8, .6, .7, .4, .3, NaN, .99, 1.3])

def ReGridMedian(AltIn, ValIn, AltOut):
    part = isfinite(ValIn)
    q = [abs(AltIn[part]-k) for k in AltOut]
    q = [nonzero(abs(k - min(k))<3*finfo(k.dtype).eps) for k in q]
    q = [ValIn[part][k] for k in q]
    return [median(k) for k in q]

B_Agrid = ReGridMedian(B_Alt, B_Val, A_Alt)    
A_Bgrid = ReGridMedian(A_Alt, A_Val, B_Alt)

我拼凑了一个检查两个值在计算机精度范围内是否相同的东西,但我觉得可能有更好的方法。无论如何,我们首先过滤掉所有不是nan的值,然后找到最接近的匹配,再检查是否有重复的最小值,最后得到这些值的中位数。

====

这样回答你的问题了吗?还是我的假设有误?

1

可以看看 numpy.interp 这个函数:

http://docs.scipy.org/doc/numpy/reference/generated/numpy.interp.html

(补充说明: numpy.interp 只提供线性插值,这显然不是提问者想要的。可以使用 scipy 的方法,比如 interp1d,并设置 kind='nearest' 来实现。)

http://docs.scipy.org/doc/scipy/reference/interpolate.html

看起来你想做的是用一个数据集的高度点来推算另一个数据集的值。这可以通过 numpy 的方法或者 scipy 的插值方法很简单地实现。

撰写回答