Qfyilyi

I have xy coordinates that represents a subject over a given space. It is referenced from another point and is therefore off centre. As in the longitudinal axes is not aligned along the x-axis.

The randomly generated ellipse below provides an indication of this:

import numpy as np

from matplotlib.pyplot import scatter



xx = np.array([-0.51, 51.2])

yy = np.array([0.33, 51.6])

means = [xx.mean(), yy.mean()]  

stds = [xx.std() / 3, yy.std() / 3]

corr = 0.8         # correlation

covs = [[stds[0]**2          , stds[0]*stds[1]*corr], 

    [stds[0]*stds[1]*corr,           stds[1]**2]] 



m = np.random.multivariate_normal(means, covs, 1000).T

scatter(m[0], m[1])

To straighten the coordinates I was thinking of applying the vector to a rotation matrix.

Would something like this work?

angle = 65.

theta = (angle/180.) * np.pi



rotMatrix = np.array([[np.cos(theta), -np.sin(theta)], 

                     [np.sin(theta),  np.cos(theta)]])

This may also seem like a silly question but is there a way to determine if the resulting vector of xy coordinates is perpendicular? Or will you just have to play around with the rotation angle?

You can use sklearn.decomposition.PCA (principal component analysis) with n_components=2 to extract the smallest angle required to rotate the point cloud such that its major axis is horizontal.

Runnable example

import numpy as np

import matplotlib.pyplot as plt

from sklearn.decomposition import PCA



np.random.seed(1)



xx = np.array([-0.51, 51.2])

yy = np.array([0.33, 51.6])

means = [xx.mean(), yy.mean()]  

stds = [xx.std() / 3, yy.std() / 3]

corr = 0.8         # correlation

covs = [[stds[0]**2,       stds[0]*stds[1]*corr], 

        [stds[0]*stds[1]*corr, stds[1]**2]]



m = np.random.multivariate_normal(means, covs, 1000)



pca = PCA(2)



# This was in my first answer attempt: fit_transform works fine, but it randomly 

# flips (mirrors) points across one of the principal axes.

# m2 = pca.fit_transform(m)



# Workaround: get the rotation angle from the PCA components and manually 

# build the rotation matrix.



# Fit the PCA object, but do not transform the data

pca.fit(m)



# pca.components_ : array, shape (n_components, n_features)

# cos theta

ct = pca.components_[0, 0]

# sin theta

st = pca.components_[0, 1]



# One possible value of theta that lies in [0, pi]

t = np.arccos(ct)



# If t is in quadrant 1, rotate CLOCKwise by t

if ct > 0 and st > 0:

    t *= -1

# If t is in Q2, rotate COUNTERclockwise by the complement of theta

elif ct < 0 and st > 0:

    t = np.pi - t

# If t is in Q3, rotate CLOCKwise by the complement of theta

elif ct < 0 and st < 0:

    t = -(np.pi - t)

# If t is in Q4, rotate COUNTERclockwise by theta, i.e., do nothing

elif ct > 0 and st < 0:

    pass



# Manually build the ccw rotation matrix

rotmat = np.array([[np.cos(t), -np.sin(t)], 

                   [np.sin(t),  np.cos(t)]])



# Apply rotation to each row of m

m2 = (rotmat @ m.T).T



# Center the rotated point cloud at (0, 0)

m2 -= m2.mean(axis=0)



fig, ax = plt.subplots()

plot_kws = {'alpha': '0.75',

            'edgecolor': 'white',

            'linewidths': 0.75}

ax.scatter(m[:, 0], m[:, 1], **plot_kws)

ax.scatter(m2[:, 0], m2[:, 1], **plot_kws)

Output

Warning: pca.fit_transform() sometimes flips (mirrors) the point cloud

The principal components can randomly come out as either positive or negative. In some cases, your point cloud may appear flipped upside down or even mirrored across one of its principal axes. (To test this, change the random seed and re-run the code until you observe flipping.) There's an in-depth discussion here (based in R, but the math is relevant). To correct this, you'd have to replace the fit_transform line with manual flipping of one or both components' signs, then multiply the sign-flipped component matrix by the point cloud array.

Indeed a very useful concept here is a linear transformation of a vector v performed by a matrix A. If you treat your scatter points as the tip of vectors originating from (0,0), then is very easy to rotate them any angle theta. A matrix that performs such rotation of theta would be

A = [[cos(theta) -sin(theta]

     [sin(theta)  cos(theta)]]

Evidently, when theta is 90 degrees this results into

A = [[ 0 1]

     [-1 0]]

And to apply the rotation you would only need to perform the matrix multiplication w = A v

With this, the current goal is to perform a matrix multiplication of the vectors stored in m with x,y tips as m[0],m[1]. The rotated vector are gonna be stored in m2. Below is the relevant code to do so. Note that I have transposed m for an easier computation of the matrix multiplication (performed with @) and that the rotation angle is 90 degress counterclockwise.

import numpy as np

import matplotlib.pyplot as plt



xx = np.array([-0.51, 51.2])

yy = np.array([0.33, 51.6])

means = [xx.mean(), yy.mean()]  

stds = [xx.std() / 3, yy.std() / 3]

corr = 0.8         # correlation

covs = [[stds[0]**2          , stds[0]*stds[1]*corr], 

    [stds[0]*stds[1]*corr,           stds[1]**2]] 



m = np.random.multivariate_normal(means, covs, 1000).T

plt.scatter(m[0], m[1])



theta_deg = 90

theta_rad = np.deg2rad(theta_deg)

A = np.matrix([[np.cos(theta_rad), -np.sin(theta_rad)],

               [np.sin(theta_rad), np.cos(theta_rad)]])

m2 = np.zeros(m.T.shape)



for i,v in enumerate(m.T):

  w = A @ v.T

  m2[i] = w

m2 = m2.T



plt.scatter(m2[0], m2[1])

This leads to the rotated scatter plot:

You can be sure that the rotated version is exactly 90 degrees counterclockwise with the linear transformation.

Edit

To find the rotation angle you need to apply in order for the scatter plot to be aligned with the x axis a good approach is to find the linear approximation of the scattered data with numpy.polyfit. This yields to a linear function by providing the slope and the intercept of the y axis b. Then get the rotation angle with the arctan function of the slope and compute the transformation matrix as before. You can do this by adding the following part to the code

slope, b = np.polyfit(m[1], m[0], 1)

x = np.arange(min(m[0]), max(m[0]), 1)

y_line = slope*x + b

plt.plot(x, y_line, color='r')

theta_rad = -np.arctan(slope)

And result to the plot you were seeking

Edit 2

Because @Peter Leimbigler pointed out that numpy.polyfit does not find the correct global direction of the scattered data, I have thought that you can get the average slope by averaging the x and y parts of the data. This is to find another slope, called slope2 (depicted in green now) to apply the rotation. So simply,

slope, b = np.polyfit(m[1], m[0], 1)

x = np.arange(min(m[0]), max(m[0]), 1)

y_line = slope*x + b

slope2 = np.mean(m[1])/np.mean(m[0])

y_line2 = slope2*x + b

plt.plot(x, y_line, color='r')

plt.plot(x, y_line2, color='g')

theta_rad = -np.arctan(slope2)

And by applying the linear transformation with the rotation matrix you get

If the slope of the two lines multiplied together is equal to -1 than they are perpendicular.
The other case this is true, is when one slope is 0 and the other is undefined (a perfectly horizontal line and a perfectly vertical line).