[GIS] 3d interpolation between two xyz coordinates

geomaticainterpolationnumpypython

I need to find the x and y coordinate on a known z coordinate based on two known xyz coordinates. In more detail: I have a csv-file with x y and z values and I need to find the place where the 0, z value is crossed.

The data looks like this:

z,x,y
0.822,174269.8332,608470.8717
0.692,174269.5191,608480.8668
0.502,174269.205,608490.8618
0.382,174268.8909,608500.8569
0.272,174268.5768,608510.852
0.172,174268.2627,608520.847
0.022,174267.9486,608530.8421
-0.262,174267.6345,608540.8372
-0.772,174267.3203,608550.8322
-1.052,174267.0062,608560.8273
-0.792,174266.6921,608570.8224
-0.382,174266.378,608580.8174

The (very basic…) illustration below illustrates what I mean. The green and the red lines are the known points and the blue line between those points is a slope that I could calculate I suspect. The yellow point with the circle around t is the x and y coordinate that I have to find:

Illustration of finding the x and y of a known z location. The green and the red lines are the known points and the blue line between those points is a slope that I could calculate I suspect. The yellow point with the circle around t is the x and y coordinate that I have to find

I have been looking at scipy.interpolate.interp2d but that returns the error:

Traceback (most recent call last):
  File "D:\ARNO\JarkusPython\ARNO\Create_Data\createStraightData.py", line 50, in <module>
    some = scipy.interpolate.interp2d(xx, yy, zz)
  File "C:\Python27\ArcGIS10.5\lib\site-packages\scipy\interpolate\interpolate.py", line 229, in __init__
    self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
  File "C:\Python27\ArcGIS10.5\lib\site-packages\scipy\interpolate\fitpack.py", line 929, in bisplrep
    raise TypeError('m >= (kx+1)(ky+1) must hold')
TypeError: m >= (kx+1)(ky+1) must hold

Which means (I believe) something as that there's not enough coordinates. Also, I think I have to feed the function the x and y and not the z value. Is this the mathematical problem that I run into?

And does anybody know a solution to this?

Best Answer

You can interpolate the two axes independently like:

Code

def find_xy(p1, p2, z):

    x1, y1, z1 = p1
    x2, y2, z2 = p2
    if z2 < z1:
        return find_xy(p2, p1, z)

    x = np.interp(z, (z1, z2), (x1, x2))
    y = np.interp(z, (z1, z2), (y1, y2))

    return x, y

Test Code:

import numpy as np

# arbitrary z
z1, x1, y1 = 0.822, 174269.8332, 608470.8717
z2, x2, y2 = 0.692, 174269.5191, 608480.8668

print(find_xy((x1, y1, z1), (x2, y2, z2), .7))
print(find_xy((x1, y1, z1), (x2, y2, z2), .8))

# z = 0
z1, x1, y1 = 0.022, 174267.9486, 608530.8421
z2, x2, y2 = -0.262, 174267.6345, 608540.8372

print(find_xy((x1, y1, z1), (x2, y2, z2), 0))

Results:

(174269.53842923077, 608480.2517169231)
(174269.78004461539, 608472.5631784615)

(174267.92426830987, 608531.6163683098)
Related Question