Is there a python xarray function, like interp, to lookup the closest coordinate value that matches a variable value?

Question

I am looking for a way to return the interpolated coordinate value of a variable at a query variable value. I think it is the opposite use to a normal interpolation/table-look up.

For example, here is a simple, normal interpolation, From this very good tutorial:

x_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
f = xr.DataArray(x_data**2, dims=['x'], coords={'x': x_data})
f.plot(marker='o')

Which produces this plot:

So a standard lookup (e.g. f.interp(x=4.5)) would look for the value of the variable at some coordinate value of x. But I want to look for the value of the x where the variable is equal to some value. For example in the above plot, what is the value of x where the variable equals 20.

I am ultimately trying to do this for a 3D field of temperature (longitude x latitude x depth). So, how do I find the value of depth, where the temperature field is closest to a value that I specify.

Many thanks in advance.

David_G · Accepted Answer · 2022-04-01 02:28:31Z

I have found a simple method for solving this using scipy.interpolate which works, though is relatively slow to compute.

from scipy import interpolate

ii = interpolate.interp1d(np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])**2,np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]))
ii(20)

Which returns array(4.44444444).

I think this is an answer, which might help someone, but not the ideal answer. Ideally this could be done in a pythonic manner, perhaps using the xarray module.

To extend this to 3D, I had to loop it within 2 for-loops, like this:

def calc_depthToTempROMS(tempfield,depthfield,queryfield,xxlims,eelims):
    ## the two input fields are the tempfield (equivalent to the coordinates) and the depthfield (equivalent to the function to be interpolated). 
    ## queryfield is the values of temperature to interpolate for depth.
    ## xxlims and eelims are the for-loop bounds, which are model specific
    for xx in xxlims:
        for ee in eelims:
            ii = interpolate.interp1d(tempfield.isel(xi_rho=xx,eta_rho=ee),depthfield.isel(xi_rho=xx,eta_rho=ee)) # make the interpolator
            depthTC[dict(eta_rho=ee,xi_rho=xx)] = ii(queryfield.isel(xi_rho=xx,eta_rho=ee)) # interpolate for each location
        if xx.values in np.round(np.linspace(0,len(xxlims.values),10)):
            print(xx.values) # print every 10th index to monitor progress
    return depthTC


depthTC = (truth.tempTC.mean(dim='ocean_time')*np.nan).load() # initialise an empty matrix
tempfield = truth.temp.mean(dim='ocean_time').load() # coordinates
depthfield = truth.z_rho0.load() # function to interpolate
queryfield = truth.tempTC.mean(dim='ocean_time').load() # target coordinate
xxlims = truth.xi_rho.load()
eelims = truth.eta_rho.load()

depthTC = calc_depthToTempROMS(tempfield,depthfield,queryfield,xxlims,eelims)

Someone else using xarray for processing earth science model data may find the above code useful.

I will leave this unaccepted for now - hopefully someone can come up with a better, xarray-based method for achieving this.

JoSeis · Accepted Answer · 2025-11-10 23:36:13Z

I'm in earth sciences, and I just implemented a similar solution with xarray


import xarray as xr
from scipy.ndimage import gaussian_filter

def extract_isosurface_depth(model:xr.DataArray, elev:xr.DataArray, threshold:float, clip:bool|None=None,smooth:int|None=1.2):
    """
    Given a DataArray for an attribute e.g. vp,vs,rho ... Extract the depth horizon 
    for an iso-value. This is useful for extracting depth to basement 
    e.g vp <= 4.5 km/s

    Parameters:
        model: Input data array attribute e.g. temp in deg celsius
        elev: Elevation attribute. Used to infill null values where the isosurface 
            does not exist
        threshold: iso-surface target value. Unit should be the same as model
        clip (optional): Can be used to avoid spurious depth values. e.g low velocity 
            zones that can cause spikes on the depth map. Unit is same as z- attribute 
        smooth (optional): Convolve result with a gaussian kernel to give a more 
            continuous result

    Example:
        bdm = extract_isosurface_depth(data["temp"],data["elevation"],25,3)
    """
    filt = model.where(model <= threshold)
    # Find the shallowest (maximum z) coordinate where the condition is met
    depth_map = filt["z"].where(filt.notnull()).max(dim="z", skipna=True)

    if elev is not None:
        # Fill the NaN values in the basement map with surface elevation
        depth_map = depth_map.fillna(elev)

    if clip:
        depth_map = depth_map.clip(max=clip)
    
    if smooth:
        depth_map.values = gaussian_filter(depth_map.values,smooth)
    
    return depth_map

In your example, if you have a temperature DataArray, you can set a single temperature value, and extract all the depth coordinates for the target value. For x-y points where the value is not detected, it returns nan by default, so I filled all the nan values with surface elevation. Operation is also vectorized so it's fast

Collectives™ on Stack Overflow

Is there a python xarray function, like interp, to lookup the closest coordinate value that matches a variable value?

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related