Соответствует ли две сетки для анализа данных, есть ли хороший алгоритм для моей проблемы?

Question

Я бы хотел сравнить два разных набора данных в 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

ai и bi содержат двойные числа, а также поля 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

Может быть, у этого даже есть имя (я думаю, что это обычная проблема), но я не знаю его и поэтому не могу его искать. Я считаю, что есть эффективный способ сделать это, кроме очевидного решения, кодирующего его самостоятельно (но я боюсь, что оно не будет эффективным, и я бы привел много ошибок).

Предпочтительно с использованием NumPy.

РЕДАКТИРОВАТЬ: Спасибо за ваш вклад для всех четырех участников. Я немного узнал, и я прошу прощения за то, что не спросил очень четко. Я был в процессе понимания проблемы. Ваши ответы указывали на использование interp1d, а этот ответ позволял мне злоупотреблять им. Я опубликую результат в ближайшее время. Я могу принять только один ответ, но любой может сделать.

Answer 1

Два предположения: 1: Вы ищете не ближайшего соседа, а все высоты в пределах некоторого диапазона.Итак, допустим, что для a1 вы хотите, чтобы все bn были в пределах 0.5 от a1 (давая вам b1 и b2 в соответствии с вашим примером).Я бы определил «ближайшего соседа» как нечто иное.

2: Вы не учитываете nan в своих медианах (numpy считает их бесконечностью согласно некоторому соглашению IEEE, но мне это кажется странным).Таким образом, согласно вашему предложению мы используем nanmedian из scipy.stats.

Я бы сделал следующее:

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 до k в A_Alt меньшеуказанный диапазон.Затем мы берем медиану тех B_Val.То же самое работает для A_Bgrid с результатами в соответствии с запросом.

== Редактировать ==

Другое предположение относительно ваших ближайших соседей: Давайте возьмем запись ближайшего соседа (или записи в случае, еслисвязи) с наименьшим абсолютным перепадом высот, не имея значения nan в качестве значения 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, затем находим наиболее близкое соответствие, затем проверяем наличие дублирующих минимумов, а затем получаем медиану этих значений.

====

Это охватывает ваш вопрос, или мои предположения неверны?

Answer 2

Вот один из способов:

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 отсортированы.

Answer 3

Это не совсем тот ответ, который вы искали, но это мой ответ 50c ...

A = {1:'a1',2:'a2',3:'a3',4:'a4'}
B = {0.7:'b1',0.9:'b2',1.7:'b3',2:'b4', 2.4:'b5'}

C = {} # result

# find altitude in A that is the closest to altitude in B
def findAltitude( altB,A):
    toto = [ ((alt-altB)**2,alt) for alt in A.keys() ]
    toto.sort()
    return toto[0][1]

#iter on each altitude of B
for altB,valueB in B.iteritems():
    altC = findAltitude( altB,A)
    if altC in C:
        C[altC].append(valueB)
    else:
        C[altC] = [valueB,]

# then do the median operation
#for altC,valueC in C.iteritems():
#   C[altC] = map( median, valueC ) # where median is your median function

print C

Это НЕ лучшее решение вообще (особенно если у вас много значений), но писать быстрее всего ...

На самом деле, это зависит от того, как хранятся ваши данные.Словарь не лучший выбор.

Более интересно / умно использовать тот факт, что ваши высоты отсортированы.Вы должны предоставить более подробную информацию о том, как хранятся ваши данные (массив с NumPy?)попробуйте что-нибудь более «умное», основываясь на том факте, что ваши высоты отсортированы.

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

# add val into C at the end of C or in the last place (depending if alt already exists in C or not)
def addto(C,val,alt):
    if C and C[-1][0]==alt:
        C[-1][1].append(valB)
    else:
        C.append( (alt,[valB,] ))



# values
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])

#intermediate list of tuple (altitude, list_of_values)
C= []

#iterator on A
Aa = iter(A_Alt)
ainf = Aa.next()
asup = Aa.next()  # two first values of A_Alt

#iterator on B
Ba = iter(B_Alt)
Bv = iter(B_Val)

# regrid
try:
    while True:
        altB = Ba.next()
        valB = Bv.next()

        # find ainf and asup in A_Alt such that ainf < altB < asup
        while asup<altB:
            try:
                ainf,asup = asup, Aa.next()
            except StopIteration:
                break

        # find closest
        if abs(ainf-altB)<=abs(asup-altB):
            addto(C, valB, ainf)
        else:
            addto(C, valB, asup)

except StopIteration:
    pass

# do the median
res = [ nanmedian(k[1]) for k in C ] 

print res

Идея состоит в том, чтобы перебрать два вектора / списка высот и найти для каждой высоты Bдве высоты А, которые его окружают.Тогда легко найти ближайший ...

Это менее читабельно, чем решение Даана, но оно должно быть более эффективным (линейным по размеру ваших данных).

Вы простонужно изменить, если ваши данные не хранятся таким образом.

Answer 4

Посмотрите на numpy.interp:

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

( EDIT : numpy.interp обеспечивает только линейную интерполяцию, которая, очевидно, не является тем, что ищет OP. Вместо этого используйте методы scipy, такие как interp1d, используя kind='nearest')

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

То, что вы хотите сделать, это использовать высотные точки одного набора данных для интерполяции значений другого. Это можно сделать довольно легко либо с помощью метода numpy, либо с помощью одного из методов интерполяции scipy.

Answer 5

Один из способов охватить второй случай (таблица B - A, то есть от нескольких высот до нескольких высот) заключается в следующем:

Функция экстраполяции (от здесь )

from scipy.interpolate import interp1d

def extrap1d(interpolator):
    xs = interpolator.x
    ys = interpolator.y

    def pointwise(x):
        if x < xs[0]:
            return ys[0]
        elif x > xs[-1]:
            return ys[-1]
        else:
            return interpolator(x)

    def ufunclike(xs):
        return array(map(pointwise, array(xs)))

    return ufunclike

Значения

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])

Фактическое сопоставление:

f_i = interp1d(A_Alt, A_Val, kind='nearest')
f_x = extrap1d(f_i)

f_x(B_Alt)

Вывод:

array([ 0.33,  0.33,  0.5 ,  0.5 ,  0.5 ,  0.6 ,  0.6 ,  0.6 ,  0.8 ,  0.8 ])