обратная интерполяция xarray (по координатам, а не по данным)

Question

У меня есть следующий массив данных

arr = xr.DataArray([[0.33, 0.25],[0.55, 0.60],[0.85, 0.71],[0.92,0.85],[1.50,0.96],[2.5,1.1]],[('x',[0.25,0.5,0.75,1.0,1.25,1.5]),('y',[1,2])])

Это дает следующий вывод

<xarray.DataArray (x: 6, y: 2)>
array([[0.33, 0.25],
       [0.55, 0.6 ],
       [0.85, 0.71],
       [0.92, 0.85],
       [1.5 , 0.96],
       [2.5 , 1.1 ]])
Coordinates:
  * x        (x) float64 0.25 0.5 0.75 1.0 1.25 1.5
  * y        (y) int32 1 2

или отсортирован ниже с x и выводом (z) рядом для удобства.

x         z (y=1)   z(y=2)
0.25      0.33      0.25
0.50      0.55      0.60
0.75      0.85      0.71
1.00      0.92      0.85
1.25      1.50      0.96
1.50      2.50      1.10

Данные, которые я имею, являются результатом нескольких входных значений. Одним из них является значение х. Есть несколько других измерений (таких как y) для других входных значений. Я хочу знать, когда мое выходное значение (z) становится больше, чем 1,00, сохраняя другие размеры фиксированными и меняя значение x. В приведенном выше 2-мерном примере я хотел бы получить ответ [1,03 1,32]. Потому что значение 1,03 для x даст мне 1,00 для z, когда y = 1, и значение 1,32 для x даст мне 1,00 для z, когда y = 2.

edit: Поскольку выходной сигнал z будет расти с увеличением x, есть только одна точка, в которой z будет иметь значение 1.0.

Есть ли эффективный способ добиться этого с помощью xarray? Моя фактическая таблица намного больше и имеет 4 входа (измерения).

Спасибо за любую помощь!

Answer 1

xarray имеет очень удобную функцию для этого: xr.interp, которая будет выполнять кусочно-линейную интерполяцию xarray.

В вашем случае вы можете использовать ее для получения кусочно интерполяция точек (x, y1) и (x, y1). Как только это будет сделано, единственное, что остается сделать, - это получить значение вашего интерполированного массива x, связанного со значением закрытия вашего интерполированного массива y1/y2/.., до целевого числа (1.00 в вашем примере).

Вот как это может выглядеть:

y_dims = [0, 1,] 
target_value = 1.0
# create a 'high resolution` version of your data array:
arr_itp = arr.interp(x=np.linspace(arr.x.min(), arr.x.max(), 10000))
for y in y_dims:
    # get the index of closest data
    x_closest = np.abs(arr_itp.isel(y=y) - target_value).argmin()
    print(arr_itp.isel(y=y, x=x_closest))

>>> <xarray.DataArray ()>
>>> array(0.99993199)
>>> Coordinates:
>>>     y        int64 1
>>>     x        float64 1.034
>>> <xarray.DataArray ()>
>>> array(1.00003)
>>> Coordinates:
>>>     y        int64 2
>>>     x        float64 1.321

---

Хотя это работает, это не очень эффективный способ решения проблемы , и здесь есть две причины почему бы и нет:

1. Использование xr.interp делает кусочную интерполяцию всего DataArray. Однако нам всегда требуется интерполяция только между двумя точками, ближайшими к целевому значению.
2. Здесь интерполяция - это прямая линия между двумя точками. Но если мы знаем одну координату точки на этой линии (y = 1,00), то мы можем просто вычислить другую координату, решив линейное уравнение прямой линии, и проблема решается с помощью нескольких арифметических операций c.

С учетом этих причин мы можем разработать более эффективное решение вашей проблемы:

# solution of linear function between two points (2. reason)
def lin_itp(p1,p2,tv):
    """Get x coord of point on line

    Determine the x coord. of a point (x, target_value) on the line
    through the points p1, p2.

    Approach:
      - parametrize x, y between p1 and p2: 
          x = p1[0] + t*(p2[0]-p1[0])
          y = p1[1] + t*(p2[1]-p1[1])
      - set y = tv and resolve 2nd eqt for t
          t = (tv - p1[1]) / (p2[1] - p1[1])
      - replace t in 1st eqt with solution for t
          x = p1[0] + (tv - p1[1])*(p2[0] - p1[0])/(p2[1] - p1[1])
    """
    return float(p1[0] + (tv - p1[1])*(p2[0] - p1[0])/(p2[1] - p1[1])) 

# target value:
t_v = 1.0
for y in [0, 1]:
    arr_sd = arr.isel(y=y)
    # get index for the value closest to the target value (but smaller)
    s_udim = int(xr.where(arr_sd - t_v <=0, arr_sd, arr_sd.min()).argmax())
    # I'm explicitly defining the two points here
    ps_itp = arr_sd[s_udim:s_udim+2]
    p1, p2 = (ps_itp.x[0], ps_itp[0]), (ps_itp.x[1], ps_itp[1])
    print(lin_itp(p1,p2,t_v))

>>> 1.0344827586206897
>>> 1.3214285714285714

Answer 2

Проблема, с которой я столкнулся при ответе Джоджо, состоит в том, что его трудно расширить во многих измерениях и сохранить структуру xarray. Поэтому я решил изучить это подробнее. Я использовал некоторые идеи из кода Jojo, чтобы сделать ответ ниже.

Я делаю два массива, один с условием, что значения меньше, чем я ищу, и один с условием, что они должны быть больше. Я сдвигаю второй в направлении x на минус 1. Теперь я объединяю их в формулу нормальной линейной интерполяции. Два массива имеют только значения, перекрывающиеся на «краю» условия. Если не сместить -1, никакие значения не будут перекрываться. В последней строке я суммирую по x-направлению и, поскольку все другие значения NaN, я извлекаю правильное значение и удаляю x-направление из DataArray в процессе.

def interpolate_dimension_x(arr, target_value, step):
    M0 = arr.where(arr - target_value <= 0)
    M1 = arr.where(arr - target_value > 0).shift(x=-1)

    work_mat = M0.x + step * (target_value - M0) / (M1 - M0)

    return work_mat.sum(dim='x')

interpolate_dimension_x(arr, 1, 0.25)

>>> <xarray.DataArray (y: 2)>
array([1.034483, 1.321429])
Coordinates:
  * y        (y) int32 1 2

У меня есть некоторые недостатки с моим кодом. Код работает, только если M0 и M1 находят значение, которое удовлетворяет условию. В противном случае все значения в этой строке будут установлены на NaN. Чтобы избежать проблем с M0, я решил, что значения х должны начинаться с 0, поскольку мое целевое значение всегда больше 0. Чтобы избежать проблем с M1, я выбираю значения х достаточно большими, чтобы знать, что мои значения там , Естественно, это не идеальные решения и могут нарушить код. Если я получу немного больше опыта с xarray и python, я мог бы переписать. В итоге у меня есть следующие пункты, которые я хотел бы решить:

• Как экстраполировать значения вне x-диапазона? В настоящее время я просто гарантирую, что мой x-диапазон достаточно большой, чтобы ответы попадали в него.
• Как сделать код устойчивым для переменного размера шага?
• Как сделать код так, чтобы мое измерение можно было выбирать динамически (теперь это работает только для 'x')
• Любые оптимизации приветствуются.