interpolação reversa de xarray (na coordenada, não nos dados)

Eu tenho o seguinte DataArray

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

Isso fornece a seguinte saída

<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

ou classificados abaixo com xe saída (z) um ao lado do outro por conveniência.

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

Os dados que tenho são o resultado de vários valores de entrada. Um deles é o valor x. Existem várias outras dimensões (como y) para outros valores de entrada. Quero saber quando meu valor de saída (z) está crescendo acima de 1,00, mantendo as outras dimensões fixas e variando o valor x. No exemplo bidimensional acima, eu gostaria de obter a resposta [1.03 1.32]. Como um valor de 1,03 para x me dará 1,00 para z quando y = 1 e um valor de 1,32 para x me dará 1,00 para z quando y = 2.

edit: Como a saída z aumentará com o aumento de x, existe apenas um ponto em que z terá 1,0 como saída.

Existe alguma maneira eficiente de conseguir isso com o xarray? Minha tabela atual é muito maior e possui 4 entradas (dimensões).

Obrigado por qualquer ajuda!

O xarray tem uma função muito útil para isso: o xr.interpque fará uma interpolação linear por partes de um xarray.

No seu caso, você pode usá-lo para obter uma interpolação por partes dos pontos (x, y1) e (x, y1). Feito isso, a única coisa que resta a fazer é obter o valor do seu xarray interpolado associado ao valor de fechamento do seu y1/y2/..array interpolado para o número de destino (1,00 no seu exemplo).

Aqui está como isso pode parecer:

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

Enquanto isso funciona, não é uma maneira realmente eficiente de abordar o problema e aqui estão duas razões pelas quais não:

1. O uso do xr.interp faz uma interpolação por partes de todo o DataArray. No entanto, sempre precisamos da interpolação entre os dois pontos mais próximos do seu valor-alvo.
2. Aqui, uma interpolação é uma linha reta entre 2 pontos. Mas se conhecemos uma coordenada de um ponto nessa linha (y = 1,00), podemos simplesmente calcular a outra coordenada resolvendo a equação linear da linha reta e o problema é resolvido em algumas operações aritméticas.

Levando em consideração esses motivos, podemos desenvolver uma solução mais eficiente para o seu problema:

# 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

O problema que tive com a resposta do jojo é que é difícil expandi-lo em várias dimensões e manter a estrutura de raio-x. Por isso, decidi aprofundar isso. Eu usei algumas idéias do código de jojo para fazer a resposta abaixo.

Eu faço duas matrizes, uma com a condição de que os valores sejam menores do que o que procuro e outra com a condição de que eles precisam ser maiores. Mudo o segundo na direção x em menos 1. Agora os combino em uma fórmula de interpolação linear normal. As duas matrizes têm apenas valores sobrepostos na 'borda' da condição. Se não for deslocado por -1, nenhum valor será sobreposto. Na linha final, NaNsomarei a direção x e, como todos os outros valores são , extraio o valor correto e removo a direção x do DataArray no processo.

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

Eu tenho algumas desvantagens no meu código. O código funciona apenas se M0 e M1 encontrarem um valor que atenda à condição. Caso contrário, todos os valores nessa linha serão definidos como NaN. Para evitar problemas com M0, decidi que os valores-x iniciassem em 0, pois meu valor-alvo é sempre maior que 0. Para evitar problemas com M1, escolho meus valores de x grandes o suficiente para saber que meus valores estão lá. . Naturalmente, essas não são soluções ideais e podem quebrar o código. Se eu tiver um pouco mais de experiência com xarray e python, posso reescrever. Em resumo, tenho os seguintes itens que gostaria de resolver:

• Como extrapolar valores fora do intervalo x? Atualmente, estou apenas garantindo que meu intervalo x seja grande o suficiente para que as respostas caiam dentro dele.
• Como tornar o código robusto para uma variável de tamanho escalonado?
• Como criar o código para que minha dimensão possa ser escolhida dinamicamente (agora funciona apenas para 'x')
• Quaisquer otimizações são apreciadas.