xarray reverse interpolation (on coordinate, not on data)

xarray has a very handy function for this: xr.interp which will do a piecewise linear interpolation of an xarray.

In your case, you can use it to obtain a piecewise interpolation of the (x, y1) and (x, y1) points. Once this is done, the only thing that remains to do is getting the value of your interpolated x array that is associated to the closes value of your interpolated y1/y2/.. array to the target number (1.00 in your example).

Here is how this could look like:

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

---

While this works, it is not a really efficient way to approach the problem and here are 2 reasons why not:

1. Using xr.interp does a piecewise interpolation of the entire DataArray. Howerver, we only ever need the interpolation between the two points closest to your target value.
2. Here, an interpolation is a straight line between 2 points. But if we know one coordinate of a point on that line (y=1.00) then we can simply calculate the other coordinate by resolving the linear equation of the straight line and the problem is solved in a few arithmetic operations.

Taking these reasons into account we can develop a more efficient solution to your problem:

# 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