Python >> Tutoriel Python >  >> Python Tag >> Array

interpolation inverse xarray (sur les coordonnées, pas sur les données)

xarray a une fonction très pratique pour cela :xr.interp qui fera une interpolation linéaire par morceaux d'un xarray.

Dans votre cas, vous pouvez l'utiliser pour obtenir une interpolation par morceaux des points (x, y1) et (x, y1). Une fois cela fait, il ne reste plus qu'à obtenir la valeur de votre x tableau associé à la valeur de fermeture de votre y1/y2/.. interpolé tableau au nombre cible (1,00 dans votre exemple).

Voici à quoi cela pourrait ressembler :

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


Bien que cela fonctionne, ce n'est pas un moyen vraiment efficace d'aborder le problème et voici 2 raisons pourquoi pas :

  1. L'utilisation de xr.interp effectue une interpolation par morceaux de l'ensemble du DataArray. Cependant, nous n'avons besoin que de l'interpolation entre les deux points les plus proches de votre valeur cible.
  2. Ici, une interpolation est une ligne droite entre 2 points. Mais si nous connaissons une coordonnée d'un point sur cette ligne (y =1,00), nous pouvons simplement calculer l'autre coordonnée en résolvant l'équation linéaire de la ligne droite et le problème est résolu en quelques opérations arithmétiques.

En tenant compte de ces raisons, nous pouvons développer une solution plus efficace à votre problème :

# 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