I have three 2D-arrays of same shape, lets call them theta, phi and A. Let theta and phi be angles to the normal vector seen from different distances on a surface:
size = 100 # this value is fixed
x = np.arange(-size, size)
y = np.arange(-size, size)
xx, yy = np.meshgrid(xx, yy)
theta = np.arctan((xx**2+yy**2)**0.5 / 100) # angle from distance 100
phi = np.arctan((xx**2+yy**2)**0.5 / 1000) # angle from distance 1000
And let A be a 2D map of measured values where the x axis is theta and the y axis is phi in known and linear steps (which is actually not the same shape as theta and phi). What I need is the values of A(theta,phi) expressed as A(x,y). It seems I can't figure out how to convert A(theta,phi) to A(x,y), even though I know both theta(x,y) and phi(x,y).
What I tried: Via scipy.interpolate.interp2d I can map A to the same number of rows and columns as theta and phi. Now I can iterate over the indices and guess/round the best matching indices in my array
B = np.zeros(A.shape)
for i in range(0,A.shape[0]):
for j in range(0,A.shape[1]):
B[i,j] = A[int(f_theta*theta[i,j]),int(f_phi*phi[i,j])]
Where f_theta and f_phi are prefactors determined by the measured step length of a index step. This looks like very bad and inefficient coding to me, and like a crude approximation of what I actually want to do (which is an inverse interpolation mapping?). It reminds me of lookup-tables, coordinate transforms and interpolation, but with none of those keywords I found something appropriate to solve the problem. My python experience shouts out that there will be a module/function for this which I do not know.
Edit about restrictions: the range of the axis in A(theta, phi) is larger than the range of theta(x,y) and phi(x,y), such that a mapped value always exists. I do not need to map B back to A, so there are no issues with missing values. Many values in the map A(theta, phi) would never be used.
Edit about clarity: I will give an example with small matrices, hoping to clarify things:
# phi given in degrees
phi = np.array([
[2,1,2],
[1,0,1],
[2,1,2],
])
# theta given in degrees
theta = np.array([
[6,4,6],
[4,0,5],
[6,5,6],
])
# A[0,0] is the value at phi=0deg, theta=0deg
# A[0,1] is the value at phi=1deg, theta=0deg
# A[1,1] is the value at phi=1deg, theta=1deg etc
# this is a toy example, the actual A cannot be constructed by a simple rule
A = np.array([
[0.0,0.1,0.2],
[1.0,1.1,1.2],
[2.0,2.1,2.2],
[3.0,3.1,3.2],
[4.0,4.1,4.2],
[5.0,5.1,5.2],
[6.0,6.1,6.2],
])
# what I want to reach:
B = [[6.2,4.1,6.2],
[4.1,0.0,5.1],
[6.2,5.1,6.2]]
I need to clarify that I did some simplifications here:
1) For a given theta I can check the corresponding phi by looking in the table: theta[i,j] corresponds to phi[i,j]. But the construction of the example is too simplified, they do not e.g. share the same origin, it is noisy data and I therefore cannot give an analytical expression theta(phi) or phi(theta)
2) The values in my actual theta and phi are floats, and my actual A measures in non-integer steps too (e.g. 0.45 deg per step in theta direction, 0.2 deg per step in phi direction)
3) in principle, as there is a strict relation between theta and phi, I only need a specific 1D "trace" of values A to find B, but I do not understand how to find this trace, nor how to create B out of the trace. This trace in the example would be [A[0,0],A[4,1],A[5,1],A[6,2]] = [0.0,4.1,5.1,6.2]
You can do, e.g., a bilinear interpolation:
from scipy.interpolate import interpn
delta = [1.0, 1.0] # theta, phi
points = [np.arange(s)*d for s, d in zip(A.shape, delta)]
xi = np.stack((theta, phi), axis = -1)
B = interpn(points, A, xi)
This gives:
print(B)
[[6.2 4.1 6.2]
[4.1 0. 5.1]
[6.2 5.1 6.2]]