Proper tilting of input antenna positions when z-component is non-zero #24
Conversation
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## main #24 +/- ##
===========================================
+ Coverage 80.59% 93.12% +12.52%
===========================================
Files 5 4 -1
Lines 201 189 -12
Branches 47 44 -3
===========================================
+ Hits 162 176 +14
+ Misses 31 7 -24
+ Partials 8 6 -2
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. |
| def get_plane_to_xy_rotation_matrix(antvecs): | ||
| """ | ||
| Compute the rotation matrix that projects the antenna positions onto the xy-plane. | ||
| This function is used to rotate the antenna positions so that they lie in the xy-plane. | ||
|
|
||
| Parameters: | ||
| ---------- | ||
| antvecs: np.array | ||
| Array of antenna positions in the form (Nants, 3). | ||
|
|
||
| Returns: | ||
| ------- | ||
| rotation_matrix: np.array | ||
| Rotation matrix that projects the antenna positions onto the xy-plane of shape (3, 3). | ||
| """ | ||
| # Fit a plane to the antenna positions | ||
| antx, anty, antz = antvecs.T | ||
| basis = np.array([antx, anty, np.ones_like(antz)]).T | ||
| plane, res, rank, s = linalg.lstsq(basis, antz) | ||
|
|
||
| # Project the antenna positions onto the plane | ||
| slope_x, slope_y, z_offset = plane | ||
|
|
||
| # Plane is already approximately aligned with the xy-axes, | ||
| # return identity rotation matrix | ||
| if np.isclose(slope_x, 0) and np.isclose(slope_y, 0.0): | ||
| return np.eye(3) | ||
|
|
||
| # Normalize the normal vector |
There was a problem hiding this comment.
This looks like it probably works, but I also had gotten a working solution that might be a little more compact (maybe faster, not sure):
from scipy.spatial.transform import Rotation
mean = np.mean(antvecs, axis=0)
svd = np.linalg.svd(antvecs.T - mean[None].T)
# Get the normal vector of the plane
normal = svd[0][:, -1]
cross = np.cross(normal, np.array([0,0,1]))
rot = Rotation.from_rotvec(cross)
return rot.as_matrix()There was a problem hiding this comment.
Or you could just return rot itself, and use rot.apply(antvecs) to get new antenna vectors that are rotated to the plane (they will still have whatever offset mean they started with).
There was a problem hiding this comment.
I like the compactness of the solution, but I'm finding that it's not reducing the z-component as the existing solution. Maybe I'm doing something wrong here?
There was a problem hiding this comment.
I would think the mean abs of the z-component would actually be expected to be about 12, after rotation, because at least in the compact method, it only applies a rotation around the mean point, which means the output points should all have z-components equivalent to that mean. Here it looks like that mean should be in the ballpark of 12. You can also apply an offset to bring the points back down to the z=0 plane, but I don't think that is necessary.
There was a problem hiding this comment.
I think this is good! I made one comment about potentially making the rotation matrix function simpler. Other than that, it would be good to see some more testing. Definitely I'd like to see some unit-tests of the rotation matrix function (like, put some co-planar but rotated positions in, and make sure the output has zero z-component).
| def get_plane_to_xy_rotation_matrix(antvecs): | ||
| """ | ||
| Compute the rotation matrix that projects the antenna positions onto the xy-plane. | ||
| This function is used to rotate the antenna positions so that they lie in the xy-plane. | ||
|
|
||
| Parameters: | ||
| ---------- | ||
| antvecs: np.array | ||
| Array of antenna positions in the form (Nants, 3). | ||
|
|
||
| Returns: | ||
| ------- | ||
| rotation_matrix: np.array | ||
| Rotation matrix that projects the antenna positions onto the xy-plane of shape (3, 3). | ||
| """ | ||
| # Fit a plane to the antenna positions | ||
| antx, anty, antz = antvecs.T | ||
| basis = np.array([antx, anty, np.ones_like(antz)]).T | ||
| plane, res, rank, s = linalg.lstsq(basis, antz) | ||
|
|
||
| # Project the antenna positions onto the plane | ||
| slope_x, slope_y, z_offset = plane | ||
|
|
||
| # Plane is already approximately aligned with the xy-axes, | ||
| # return identity rotation matrix | ||
| if np.isclose(slope_x, 0) and np.isclose(slope_y, 0.0): | ||
| return np.eye(3) | ||
|
|
||
| # Normalize the normal vector |
There was a problem hiding this comment.
Or you could just return rot itself, and use rot.apply(antvecs) to get new antenna vectors that are rotated to the plane (they will still have whatever offset mean they started with).
| # Check that method is robust to errors | ||
| rng = np.random.default_rng(42) | ||
| random_antvecs = antvecs.copy() | ||
| random_antvecs[:, -1] += rng.standard_normal(100) | ||
|
|
||
| # Rotate the array to the xy-plane | ||
| rotation_matrix = utils.get_plane_to_xy_rotation_matrix( | ||
| random_antvecs | ||
| ) | ||
| rm_antvecs = np.dot(rotation_matrix.T, antvecs.T) | ||
|
|
||
| # Check that all elements of the z-axis within 5-sigma of zero | ||
| np.testing.assert_array_less(np.abs(rm_antvecs[-1]), 5) |
There was a problem hiding this comment.
I would make this a new test with a different name.
This PR implements a better algorithm for flattening an array that is sloped. I've added a new function with determines the slope of the array and rotates the array and sky accordingly.