I have a line and some points that are on that are on that line in 3D space. I know there is a certain amount of error in the the point, but the error only extends perpendicular to the line. To view this I'd like to have disks with a radius of the error that are centered on the line and orthogonal to the direction of the line. I found this solution but I can't get it to quite work.
If I run the code and state I want the output normal to the 'z' axis I get waht I woudl expect. Disks with a given radius on the line and oriented on z axis.
pathpatch_2d_to_3d(p, z=z,normal='z')
I would like the disks rotated. In order to find the well vector at that point I'm using a point close by using that vector. This is the vector that I am putting as normal=(vx,vy,vz) but when I do that, the disks are not even on the chart. I'm not certain where I'm going wrong. Does anyone have any advice?
Here is my code.
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, PathPatch
from mpl_toolkits.mplot3d import Axes3D
import mpl_toolkits.mplot3d.art3d as art3d
import numpy as np
from scipy.interpolate import interp1d
md,wellz,wellx,welly=np.genfromtxt("./well.csv",delimiter=",",unpack=True)
# Building interpolation function that map a measured depth to its repsective x,y,z coordinates
fz = interp1d(md,wellz)
fx = interp1d(md,wellx)
fy = interp1d(md,welly)
pointDepth = np.array([15790,15554,15215,14911,14274,13927,13625,13284,12983,12640,12345,12004,11704,11361,11061,10717,10418,10080,9771])
def rotation_matrix(d):
"""
Calculates a rotation matrix given a vector d. The direction of d
corresponds to the rotation axis. The length of d corresponds to
the sin of the angle of rotation.
Variant of: http://mail.scipy.org/pipermail/numpy-discussion/2009-March/040806.html
"""
sin_angle = np.linalg.norm(d)
if sin_angle == 0:
return np.identity(3)
d = d/sin_angle
eye = np.eye(3)
ddt = np.outer(d, d)
skew = np.array([[ 0, d[2], -d[1]],
[-d[2], 0, d[0]],
[d[1], -d[0], 0]], dtype=np.float64)
M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew
return M
def pathpatch_2d_to_3d(pathpatch, z = 0, normal = 'z'):
"""
Transforms a 2D Patch to a 3D patch using the given normal vector.
The patch is projected into they XY plane, rotated about the origin
and finally translated by z.
"""
if type(normal) is str: #Translate strings to normal vectors
index = "xyz".index(normal)
normal = np.roll((1,0,0), index)
normal = normal/np.linalg.norm(normal) #Make sure the vector is normalised
path = pathpatch.get_path() #Get the path and the associated transform
trans = pathpatch.get_patch_transform()
path = trans.transform_path(path) #Apply the transform
pathpatch.__class__ = art3d.PathPatch3D #Change the class
pathpatch._code3d = path.codes #Copy the codes
pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color
verts = path.vertices #Get the vertices in 2D
d = np.cross(normal, (0, 0, 1)) #Obtain the rotation vector
M = rotation_matrix(d) #Get the rotation matrix
pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts])
def pathpatch_translate(pathpatch, delta):
"""
Translates the 3D pathpatch by the amount delta.
"""
pathpatch._segment3d += delta
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(wellx,welly,wellz,c='k')
for n,pd in enumerate(pointDepth):
x,y,z = fx(pd),fy(pd),fz(pd)
# figure out a vector from the point
vx,vy,vz = (fx(pd-10)-x),(fy(pd-10)-y),(fz(pd-10)-z)
#Draw Circle
p = Circle((x,y), 100)
ax.add_patch(p)
pathpatch_2d_to_3d(p, z=z,normal=(vx,vy,vz))
pathpatch_translate(p,(0,0,0))
ax.set_xlim3d(np.min(wellx),np.max(wellx))
ax.set_ylim3d(np.min(welly), np.max(welly))
ax.set_zlim3d(np.min(wellz), np.max(wellz))
plt.show()
