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.
