Revision d745e74710ab581d489e815095d0dd4ee91e9c35 authored by Bryna Hazelton on 15 September 2025, 18:00:43 UTC, committed by Jonathan Pober on 15 September 2025, 18:31:58 UTC
1 parent ea19da5
test_coordinates.py
# Copyright (c) 2024 Radio Astronomy Software Group
# Licensed under the 2-clause BSD License
"""Tests for coordinate utility functions."""
import os
import re
import numpy as np
import pytest
from astropy import units
from astropy.coordinates import EarthLocation
from pyuvdata import utils
from pyuvdata.data import DATA_PATH
from pyuvdata.testing import check_warnings
selenoids = ["SPHERE", "GSFC", "GRAIL23", "CE-1-LAM-GEO"]
try:
from lunarsky import MoonLocation
frame_selenoid = [["itrs", None]]
for snd in selenoids:
frame_selenoid.append(["mcmf", snd])
except ImportError:
frame_selenoid = [["itrs", None]]
# Earth
ref_latlonalt = (-26.7 * np.pi / 180.0, 116.7 * np.pi / 180.0, 377.8)
ref_xyz = (-2562123.42683, 5094215.40141, -2848728.58869)
# Moon
ref_latlonalt_moon = (0.6875 * np.pi / 180.0, 24.433 * np.pi / 180.0, 0.3)
ref_xyz_moon = {
"SPHERE": (1581421.43506347, 718463.12201783, 20843.2071012),
"GSFC": (1582332.08831085, 718876.84524219, 20805.18709001),
"GRAIL23": (1581855.3916402, 718660.27490195, 20836.2107652),
"CE-1-LAM-GEO": (1581905.99108228, 718683.26297605, 20806.77965693),
}
@pytest.fixture(scope="module")
def enu_ecef_info():
"""Some setup info for ENU/ECEF calculations."""
center_lat = -30.7215261207 * np.pi / 180.0
center_lon = 21.4283038269 * np.pi / 180.0
center_alt = 1051.7
# fmt: off
lats = (np.array([-30.72218216, -30.72138101, -30.7212785, -30.7210011,
-30.72159853, -30.72206199, -30.72174614, -30.72188775,
-30.72183915, -30.72100138])
* np.pi / 180.0)
lons = (np.array([21.42728211, 21.42811727, 21.42814544, 21.42795736,
21.42686739, 21.42918772, 21.42785662, 21.4286408,
21.42750933, 21.42896567])
* np.pi / 180.0)
alts = np.array([1052.25, 1051.35, 1051.2, 1051., 1051.45, 1052.04, 1051.68,
1051.87, 1051.77, 1051.06])
# used pymap3d, which implements matlab code, as a reference.
x = [5109327.46674067, 5109339.76407785, 5109344.06370947,
5109365.11297147, 5109372.115673, 5109266.94314734,
5109329.89620962, 5109295.13656657, 5109337.21810468,
5109329.85680612]
y = [2005130.57953031, 2005221.35184577, 2005225.93775268,
2005214.8436201, 2005105.42364036, 2005302.93158317,
2005190.65566222, 2005257.71335575, 2005157.78980089,
2005304.7729239]
z = [-3239991.24516348, -3239914.4185286, -3239904.57048431,
-3239878.02656316, -3239935.20415493, -3239979.68381865,
-3239949.39266985, -3239962.98805772, -3239958.30386264,
-3239878.08403833]
east = [-97.87631659, -17.87126443, -15.17316938, -33.19049252, -137.60520964,
84.67346748, -42.84049408, 32.28083937, -76.1094745, 63.40285935]
north = [-72.7437482, 16.09066646, 27.45724573, 58.21544651, -8.02964511,
-59.41961437, -24.39698388, -40.09891961, -34.70965816, 58.18410876]
up = [0.54883333, -0.35004539, -0.50007736, -0.70035299, -0.25148791, 0.33916067,
-0.02019057, 0.16979185, 0.06945155, -0.64058124]
# fmt: on
yield (
center_lat,
center_lon,
center_alt,
lats,
lons,
alts,
x,
y,
z,
east,
north,
up,
)
@pytest.fixture(scope="module")
def enu_mcmf_info():
center_lat, center_lon, center_alt = [
0.6875 * np.pi / 180.0,
24.433 * np.pi / 180.0,
0.3,
]
# Creating a test pattern of a circle of antennas, radius 500 m in ENU coordinates.
angs = np.linspace(0, 2 * np.pi, 10, endpoint=False)
enus = 500 * np.array([np.cos(angs), np.sin(angs), [0] * angs.size])
east = enus[0].tolist()
north = enus[1].tolist()
up = enus[2].tolist()
# fmt: off
lats = {
"SPHERE": np.deg2rad(
[
0.68749997, 0.69719361, 0.70318462, 0.70318462, 0.69719361,
0.68749997, 0.67780635, 0.67181538, 0.67181538, 0.67780635
]
),
"GSFC": np.deg2rad(
[
0.68749997, 0.69721132, 0.70321328, 0.70321328, 0.69721132,
0.68749997, 0.67778864, 0.67178672, 0.67178672, 0.67778864
]
),
"GRAIL23": np.deg2rad(
[
0.68749997, 0.69719686, 0.70318988, 0.70318988, 0.69719686,
0.68749997, 0.6778031 , 0.67181011, 0.67181011, 0.6778031
]
),
"CE-1-LAM-GEO": np.deg2rad(
[
0.68749997, 0.69721058, 0.70321207, 0.70321207, 0.69721058,
0.68749997, 0.67778938, 0.67178792, 0.67178792, 0.67778938
]
),
}
lons = {
"SPHERE": np.deg2rad(
[
24.44949297, 24.44634312, 24.43809663, 24.42790337, 24.41965688,
24.41650703, 24.41965693, 24.42790341, 24.43809659, 24.44634307
]
),
"GSFC": np.deg2rad(
[
24.44948348, 24.44633544, 24.43809369, 24.42790631, 24.41966456,
24.41651652, 24.41966461, 24.42790634, 24.43809366, 24.44633539
]
),
"GRAIL23": np.deg2rad(
[
24.44948845, 24.44633946, 24.43809523, 24.42790477, 24.41966054,
24.41651155, 24.41966059, 24.42790481, 24.43809519, 24.44633941
]
),
"CE-1-LAM-GEO": np.deg2rad(
[
24.44948792, 24.44633904, 24.43809507, 24.42790493, 24.41966096,
24.41651208, 24.41966102, 24.42790497, 24.43809503, 24.44633898
]
),
}
alts = {
"SPHERE": [
0.371959, 0.371959, 0.371959, 0.371959, 0.371959, 0.371959,
0.371959, 0.371959, 0.371959, 0.371959
],
"GSFC": [
0.37191758, 0.37197732, 0.37207396, 0.37207396, 0.37197732,
0.37191758, 0.37197732, 0.37207396, 0.37207396, 0.37197732
],
"GRAIL23": [
0.37193926, 0.37195442, 0.37197896, 0.37197896, 0.37195442,
0.37193926, 0.37195442, 0.37197896, 0.37197896, 0.37195442
],
"CE-1-LAM-GEO": [
0.37193696, 0.37198809, 0.37207083, 0.37207083, 0.37198809,
0.37193696, 0.37198809, 0.37207083, 0.37207083, 0.37198809
],
}
x = {
"SPHERE": [
1581214.62062477, 1581250.9080965 , 1581352.33107362,
1581480.14942611, 1581585.54088769, 1581628.24950218,
1581591.96203044, 1581490.53905332, 1581362.72070084,
1581257.32923925
],
"GSFC": [
1582125.27387214, 1582161.56134388, 1582262.984321,
1582390.80267348, 1582496.19413507, 1582538.90274956,
1582502.61527782, 1582401.1923007 , 1582273.37394822,
1582167.98248663
],
"GRAIL23": [
1581648.57720149, 1581684.86467323, 1581786.28765035,
1581914.10600283, 1582019.49746442, 1582062.2060789 ,
1582025.91860717, 1581924.49563005, 1581796.67727756,
1581691.28581598
],
"CE-1-LAM-GEO": [
1581699.17664357, 1581735.46411531, 1581836.88709243,
1581964.70544491, 1582070.0969065 , 1582112.80552098,
1582076.51804925, 1581975.09507213, 1581847.27671964,
1581741.88525806
]
}
y = {
"SPHERE": [
718918.34480718, 718829.94638063, 718601.4335154 , 718320.09035913,
718093.38043501, 718007.89922848, 718096.29765503, 718324.81052027,
718606.15367654, 718832.86360065
],
"GSFC": [
719332.06803154, 719243.66960499, 719015.15673976, 718733.81358349,
718507.10365937, 718421.62245284, 718510.02087939, 718738.53374463,
719019.8769009 , 719246.58682501
],
"GRAIL23": [
719115.4976913 , 719027.09926475, 718798.58639952, 718517.24324325,
718290.53331913, 718205.0521126 , 718293.45053915, 718521.96340439,
718803.30656066, 719030.01648477
],
"CE-1-LAM-GEO": [
719138.4857654 , 719050.08733885, 718821.57447362, 718540.23131734,
718313.52139323, 718228.0401867 , 718316.43861325, 718544.95147849,
718826.29463476, 719053.00455887
],
}
z = {
"SPHERE": [
20843.2071012 , 21137.07857037, 21318.70112664, 21318.70112664,
21137.07857037, 20843.2071012 , 20549.33563204, 20367.71307577,
20367.71307577, 20549.33563204
],
"GSFC": [
20805.18709001, 21099.05855918, 21280.68111545, 21280.68111545,
21099.05855918, 20805.18709001, 20511.31562084, 20329.69306457,
20329.69306457, 20511.31562084
],
"GRAIL23": [
20836.2107652 , 21130.08223437, 21311.70479064, 21311.70479064,
21130.08223437, 20836.2107652 , 20542.33929603, 20360.71673976,
20360.71673976, 20542.33929603
],
"CE-1-LAM-GEO": [
20806.77965693, 21100.6511261 , 21282.27368237, 21282.27368237,
21100.6511261 , 20806.77965693, 20512.90818776, 20331.28563149,
20331.28563149, 20512.90818776
],
}
# fmt: on
yield (
center_lat,
center_lon,
center_alt,
lats,
lons,
alts,
x,
y,
z,
east,
north,
up,
)
def test_XYZ_from_LatLonAlt():
"""Test conversion from lat/lon/alt to ECEF xyz with reference values."""
out_xyz = utils.XYZ_from_LatLonAlt(
ref_latlonalt[0], ref_latlonalt[1], ref_latlonalt[2]
)
# Got reference by forcing http://www.oc.nps.edu/oc2902w/coord/llhxyz.htm
# to give additional precision.
np.testing.assert_allclose(ref_xyz, out_xyz, rtol=0, atol=1e-3)
# test error checking
with pytest.raises(
ValueError,
match="latitude, longitude and altitude must all have the same length",
):
utils.XYZ_from_LatLonAlt(
ref_latlonalt[0],
ref_latlonalt[1],
np.array([ref_latlonalt[2], ref_latlonalt[2]]),
)
with pytest.raises(
ValueError,
match="latitude, longitude and altitude must all have the same length",
):
utils.XYZ_from_LatLonAlt(
ref_latlonalt[0],
np.array([ref_latlonalt[1], ref_latlonalt[1]]),
ref_latlonalt[2],
)
def test_LatLonAlt_from_XYZ():
"""Test conversion from ECEF xyz to lat/lon/alt with reference values."""
out_latlonalt = utils.LatLonAlt_from_XYZ(ref_xyz)
# Got reference by forcing http://www.oc.nps.edu/oc2902w/coord/llhxyz.htm
# to give additional precision.
np.testing.assert_allclose(ref_latlonalt, out_latlonalt, rtol=0, atol=1e-3)
pytest.raises(ValueError, utils.LatLonAlt_from_XYZ, ref_latlonalt)
# test passing multiple values
xyz_mult = np.stack((np.array(ref_xyz), np.array(ref_xyz)))
lat_vec, lon_vec, alt_vec = utils.LatLonAlt_from_XYZ(xyz_mult)
np.testing.assert_allclose(
ref_latlonalt, (lat_vec[1], lon_vec[1], alt_vec[1]), rtol=0, atol=1e-3
)
# check error if array transposed
with pytest.raises(
ValueError,
match=re.escape("The expected shape of ECEF xyz array is (Npts, 3)."),
):
utils.LatLonAlt_from_XYZ(xyz_mult.T)
# check error if only 2 coordinates
with pytest.raises(
ValueError,
match=re.escape("The expected shape of ECEF xyz array is (Npts, 3)."),
):
utils.LatLonAlt_from_XYZ(xyz_mult[:, 0:2])
# test error checking
pytest.raises(ValueError, utils.LatLonAlt_from_XYZ, ref_xyz[0:1])
@pytest.mark.parametrize("selenoid", selenoids)
def test_XYZ_from_LatLonAlt_mcmf(selenoid):
"""Test MCMF lat/lon/alt to xyz with reference values."""
pytest.importorskip("lunarsky")
lat, lon, alt = ref_latlonalt_moon
out_xyz = utils.XYZ_from_LatLonAlt(lat, lon, alt, frame="mcmf", ellipsoid=selenoid)
np.testing.assert_allclose(ref_xyz_moon[selenoid], out_xyz, rtol=0, atol=1e-3)
# test default ellipsoid
if selenoid == "SPHERE":
out_xyz = utils.XYZ_from_LatLonAlt(lat, lon, alt, frame="mcmf")
np.testing.assert_allclose(ref_xyz_moon[selenoid], out_xyz, rtol=0, atol=1e-3)
# Test errors with invalid frame
with pytest.raises(
ValueError, match="No cartesian to spherical transform defined for frame"
):
utils.XYZ_from_LatLonAlt(lat, lon, alt, frame="undef")
@pytest.mark.parametrize("selenoid", selenoids)
def test_LatLonAlt_from_XYZ_mcmf(selenoid):
"""Test MCMF xyz to lat/lon/alt with reference values."""
pytest.importorskip("lunarsky")
out_latlonalt = utils.LatLonAlt_from_XYZ(
ref_xyz_moon[selenoid], frame="mcmf", ellipsoid=selenoid
)
np.testing.assert_allclose(ref_latlonalt_moon, out_latlonalt, rtol=0, atol=1e-3)
# test default ellipsoid
if selenoid == "SPHERE":
out_latlonalt = utils.LatLonAlt_from_XYZ(ref_xyz_moon[selenoid], frame="mcmf")
np.testing.assert_allclose(ref_latlonalt_moon, out_latlonalt, rtol=0, atol=1e-3)
# Test errors with invalid frame
with pytest.raises(
ValueError, match="Cannot check acceptability for unknown frame"
):
out_latlonalt = utils.LatLonAlt_from_XYZ(ref_xyz_moon[selenoid], frame="undef")
with pytest.raises(
ValueError, match="No spherical to cartesian transform defined for frame"
):
utils.LatLonAlt_from_XYZ(
ref_xyz_moon[selenoid], frame="undef", check_acceptability=False
)
@pytest.mark.skipif(
len(frame_selenoid) > 1, reason="Test only when lunarsky not installed."
)
def test_no_moon():
"""Check errors when calling functions with MCMF without lunarsky."""
msg = "Need to install `lunarsky` package to work with MCMF frame."
with pytest.raises(ImportError, match=msg):
utils.get_lst_for_time(
[2451545.0], latitude=0, longitude=0, altitude=0, frame="mcmf"
)
msg = "Need to install `lunarsky` package to work with selenoids or MCMF frame."
with pytest.raises(ImportError, match=msg):
utils.coordinates.get_selenoids()
with pytest.raises(ImportError, match=msg):
utils.LatLonAlt_from_XYZ(ref_xyz_moon["SPHERE"], frame="mcmf")
lat, lon, alt = ref_latlonalt_moon
with pytest.raises(ImportError, match=msg):
utils.XYZ_from_LatLonAlt(lat, lon, alt, frame="mcmf")
with pytest.raises(ImportError, match=msg):
utils.ENU_from_ECEF(
None, latitude=0.0, longitude=1.0, altitude=10.0, frame="mcmf"
)
with pytest.raises(ImportError, match=msg):
utils.ECEF_from_ENU(
None, latitude=0.0, longitude=1.0, altitude=10.0, frame="mcmf"
)
msg = "Need to install `lunarsky` package to work with MoonLocations."
with pytest.raises(ImportError, match=msg):
utils.coordinates.check_surface_based_positions(
telescope_frame="mcmf", telescope_loc="foo"
)
def test_lla_xyz_lla_roundtrip():
"""Test roundtripping an array will yield the same values."""
np.random.seed(0)
lats = -30.721 + np.random.normal(0, 0.0005, size=30)
lons = 21.428 + np.random.normal(0, 0.0005, size=30)
alts = np.random.uniform(1051, 1054, size=30)
lats *= np.pi / 180.0
lons *= np.pi / 180.0
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
lats_new, lons_new, alts_new = utils.LatLonAlt_from_XYZ(xyz)
np.testing.assert_allclose(lats_new, lats, rtol=0, atol=utils.RADIAN_TOL)
np.testing.assert_allclose(lons_new, lons, rtol=0, atol=utils.RADIAN_TOL)
np.testing.assert_allclose(alts_new, alts, rtol=0, atol=1e-3)
def test_xyz_from_latlonalt(enu_ecef_info):
"""Test calculating xyz from lat lot alt."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
np.testing.assert_allclose(np.stack((x, y, z), axis=1), xyz, atol=1e-3)
def test_enu_from_ecef(enu_ecef_info):
"""Test calculating ENU from ECEF coordinates."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
enu = utils.ENU_from_ECEF(
xyz, latitude=center_lat, longitude=center_lon, altitude=center_alt
)
np.testing.assert_allclose(
np.stack((east, north, up), axis=1), enu, rtol=0, atol=1e-3
)
enu2 = utils.ENU_from_ECEF(
xyz,
center_loc=EarthLocation.from_geodetic(
lat=center_lat * units.rad,
lon=center_lon * units.rad,
height=center_alt * units.m,
),
)
np.testing.assert_allclose(enu, enu2, rtol=0, atol=1e-3)
@pytest.mark.parametrize("selenoid", selenoids)
def test_enu_from_mcmf(enu_mcmf_info, selenoid):
pytest.importorskip("lunarsky")
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_mcmf_info
)
xyz = utils.XYZ_from_LatLonAlt(
lats[selenoid], lons[selenoid], alts[selenoid], frame="mcmf", ellipsoid=selenoid
)
enu = utils.ENU_from_ECEF(
xyz,
latitude=center_lat,
longitude=center_lon,
altitude=center_alt,
frame="mcmf",
ellipsoid=selenoid,
)
np.testing.assert_allclose(
np.stack((east, north, up), axis=1), enu, rtol=0, atol=1e-3
)
enu2 = utils.ENU_from_ECEF(
xyz,
center_loc=MoonLocation.from_selenodetic(
lat=center_lat * units.rad,
lon=center_lon * units.rad,
height=center_alt * units.m,
ellipsoid=selenoid,
),
)
np.testing.assert_allclose(enu, enu2, rtol=0, atol=1e-3)
def test_invalid_frame():
"""Test error is raised when an invalid frame name is passed in."""
with pytest.raises(
ValueError, match='No ENU_from_ECEF transform defined for frame "UNDEF".'
):
utils.ENU_from_ECEF(
np.zeros((2, 3)), latitude=0.0, longitude=0.0, altitude=0.0, frame="undef"
)
with pytest.raises(
ValueError, match='No ECEF_from_ENU transform defined for frame "UNDEF".'
):
utils.ECEF_from_ENU(
np.zeros((2, 3)), latitude=0.0, longitude=0.0, altitude=0.0, frame="undef"
)
with pytest.raises(
ValueError, match="center_loc is not a supported type. It must be one of "
):
utils.ENU_from_ECEF(
np.zeros((2, 3)), center_loc=units.Quantity(np.array([0, 0, 0]) * units.m)
)
with pytest.raises(
ValueError, match="center_loc is not a supported type. It must be one of "
):
utils.ECEF_from_ENU(
np.zeros((2, 3)), center_loc=units.Quantity(np.array([0, 0, 0]) * units.m)
)
@pytest.mark.parametrize("shape_type", ["transpose", "Nblts,2", "Nblts,1"])
def test_enu_from_ecef_shape_errors(enu_ecef_info, shape_type):
"""Test ENU_from_ECEF input shape errors."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
if shape_type == "transpose":
xyz = xyz.T.copy()
elif shape_type == "Nblts,2":
xyz = xyz.copy()[:, 0:2]
elif shape_type == "Nblts,1":
xyz = xyz.copy()[:, 0:1]
# check error if array transposed
with pytest.raises(
ValueError,
match=re.escape("The expected shape of ECEF xyz array is (Npts, 3)."),
):
utils.ENU_from_ECEF(
xyz, longitude=center_lat, latitude=center_lon, altitude=center_alt
)
def test_enu_from_ecef_magnitude_error(enu_ecef_info):
"""Test ENU_from_ECEF input magnitude errors."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
# error checking
with pytest.raises(
ValueError,
match="itrs position vector magnitudes must be on the order of the "
"radius of Earth",
):
utils.ENU_from_ECEF(
xyz / 2.0, latitude=center_lat, longitude=center_lon, altitude=center_alt
)
def test_enu_from_ecef_error():
# check error no center location info passed
with pytest.raises(
ValueError,
match="Either center_loc or all of latitude, longitude and altitude "
"must be passed.",
):
utils.ENU_from_ECEF(np.array([0, 0, 0]))
with pytest.raises(
ValueError,
match="Either center_loc or all of latitude, longitude and altitude "
"must be passed.",
):
utils.ECEF_from_ENU(np.array([0, 0, 0]))
@pytest.mark.parametrize(["frame", "selenoid"], frame_selenoid)
def test_ecef_from_enu_roundtrip(enu_ecef_info, enu_mcmf_info, frame, selenoid):
"""Test ECEF_from_ENU values."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info if frame == "itrs" else enu_mcmf_info
)
if frame == "mcmf":
lats = lats[selenoid]
lons = lons[selenoid]
alts = alts[selenoid]
loc_obj = MoonLocation.from_selenodetic(
lat=center_lat * units.rad,
lon=center_lon * units.rad,
height=center_alt * units.m,
ellipsoid=selenoid,
)
else:
loc_obj = EarthLocation.from_geodetic(
lat=center_lat * units.rad,
lon=center_lon * units.rad,
height=center_alt * units.m,
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts, frame=frame, ellipsoid=selenoid)
enu = utils.ENU_from_ECEF(
xyz,
latitude=center_lat,
longitude=center_lon,
altitude=center_alt,
frame=frame,
ellipsoid=selenoid,
)
# check that a round trip gives the original value.
xyz_from_enu = utils.ECEF_from_ENU(
enu,
latitude=center_lat,
longitude=center_lon,
altitude=center_alt,
frame=frame,
ellipsoid=selenoid,
)
np.testing.assert_allclose(xyz, xyz_from_enu, rtol=0, atol=1e-3)
xyz_from_enu2 = utils.ECEF_from_ENU(enu, center_loc=loc_obj)
np.testing.assert_allclose(xyz_from_enu, xyz_from_enu2, rtol=0, atol=1e-3)
if selenoid == "SPHERE":
enu = utils.ENU_from_ECEF(
xyz,
latitude=center_lat,
longitude=center_lon,
altitude=center_alt,
frame=frame,
)
# check that a round trip gives the original value.
xyz_from_enu = utils.ECEF_from_ENU(
enu,
latitude=center_lat,
longitude=center_lon,
altitude=center_alt,
frame=frame,
)
np.testing.assert_allclose(xyz, xyz_from_enu, rtol=0, atol=1e-3)
@pytest.mark.parametrize("shape_type", ["transpose", "Nblts,2", "Nblts,1"])
def test_ecef_from_enu_shape_errors(enu_ecef_info, shape_type):
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
enu = utils.ENU_from_ECEF(
xyz, latitude=center_lat, longitude=center_lon, altitude=center_alt
)
if shape_type == "transpose":
enu = enu.copy().T
elif shape_type == "Nblts,2":
enu = enu.copy()[:, 0:2]
elif shape_type == "Nblts,1":
enu = enu.copy()[:, 0:1]
# check error if array transposed
with pytest.raises(
ValueError, match=re.escape("The expected shape of the ENU array is (Npts, 3).")
):
utils.ECEF_from_ENU(
enu, latitude=center_lat, longitude=center_lon, altitude=center_alt
)
def test_ecef_from_enu_single(enu_ecef_info):
"""Test single coordinate transform."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
# check passing a single value
enu_single = utils.ENU_from_ECEF(
xyz[0, :], latitude=center_lat, longitude=center_lon, altitude=center_alt
)
np.testing.assert_allclose(
np.array((east[0], north[0], up[0])), enu_single, rtol=0, atol=1e-3
)
def test_ecef_from_enu_single_roundtrip(enu_ecef_info):
"""Test single coordinate roundtrip."""
(center_lat, center_lon, center_alt, lats, lons, alts, x, y, z, east, north, up) = (
enu_ecef_info
)
xyz = utils.XYZ_from_LatLonAlt(lats, lons, alts)
# check passing a single value
enu = utils.ENU_from_ECEF(
xyz, latitude=center_lat, longitude=center_lon, altitude=center_alt
)
enu_single = utils.ENU_from_ECEF(
xyz[0, :], latitude=center_lat, longitude=center_lon, altitude=center_alt
)
np.testing.assert_allclose(
np.array((east[0], north[0], up[0])), enu[0, :], rtol=0, atol=1e-3
)
xyz_from_enu = utils.ECEF_from_ENU(
enu_single, latitude=center_lat, longitude=center_lon, altitude=center_alt
)
np.testing.assert_allclose(xyz[0, :], xyz_from_enu, rtol=0, atol=1e-3)
def test_mwa_ecef_conversion():
"""
Test based on comparing the antenna locations in a Cotter uvfits file to
the antenna locations in MWA_tools.
"""
test_data_file = os.path.join(DATA_PATH, "mwa128_ant_layouts.npz")
f = np.load(test_data_file)
# From the STABXYZ table in a cotter-generated uvfits file, obsid = 1066666832
xyz = f["stabxyz"]
# From the East/North/Height columns in a cotter-generated metafits file,
# obsid = 1066666832
enh = f["ENH"]
# From a text file antenna_locations.txt in MWA_Tools/scripts
txt_topo = f["txt_topo"]
# From the unphased uvw coordinates of obsid 1066666832, positions relative
# to antenna 0
# these aren't used in the current test, but are interesting and might help
# with phasing diagnosis in the future
uvw_topo = f["uvw_topo"]
# Sky coordinates are flipped for uvw derived values
uvw_topo = -uvw_topo
uvw_topo += txt_topo[0]
# transpose these arrays to get them into the right shape
txt_topo = txt_topo.T
uvw_topo = uvw_topo.T
# ARRAYX, ARRAYY, ARRAYZ in ECEF frame from Cotter file
arrcent = f["arrcent"]
lat, lon, alt = utils.LatLonAlt_from_XYZ(arrcent)
# The STABXYZ coordinates are defined with X through the local meridian,
# so rotate back to the prime meridian
new_xyz = utils.ECEF_from_rotECEF(xyz.T, lon)
# add in array center to get real ECEF
ecef_xyz = new_xyz + arrcent
enu = utils.ENU_from_ECEF(ecef_xyz, latitude=lat, longitude=lon, altitude=alt)
np.testing.assert_allclose(enu, enh, rtol=0, atol=1e-3)
# test other direction of ECEF rotation
rot_xyz = utils.rotECEF_from_ECEF(new_xyz, lon)
np.testing.assert_allclose(rot_xyz.T, xyz, rtol=0, atol=1e-3)
def test_hpx_latlon_az_za():
zenith_angle = np.deg2rad(np.linspace(0, 90, 10))
azimuth = np.deg2rad(np.linspace(0, 360, 36, endpoint=False))
az_mesh, za_mesh = np.meshgrid(azimuth, zenith_angle)
hpx_lat = np.deg2rad(np.linspace(90, 0, 10))
hpx_lon = np.deg2rad(np.linspace(0, 360, 36, endpoint=False))
lon_mesh, lat_mesh = np.meshgrid(hpx_lon, hpx_lat)
with pytest.raises(
ValueError, match="shapes of zenith_angle and azimuth values must match."
):
utils.coordinates.zenithangle_azimuth_to_hpx_latlon(zenith_angle, azimuth)
calc_lat, calc_lon = utils.coordinates.zenithangle_azimuth_to_hpx_latlon(
za_mesh, az_mesh
)
np.testing.assert_allclose(calc_lat, lat_mesh, rtol=0, atol=utils.RADIAN_TOL)
np.testing.assert_allclose(calc_lon, lon_mesh, rtol=0, atol=utils.RADIAN_TOL)
with pytest.raises(
ValueError, match="shapes of hpx_lat and hpx_lon values must match."
):
utils.coordinates.hpx_latlon_to_zenithangle_azimuth(hpx_lat, hpx_lon)
calc_za, calc_az = utils.coordinates.hpx_latlon_to_zenithangle_azimuth(
lat_mesh, lon_mesh
)
np.testing.assert_allclose(calc_za, za_mesh, rtol=0, atol=utils.RADIAN_TOL)
np.testing.assert_allclose(calc_az, az_mesh, rtol=0, atol=utils.RADIAN_TOL)
@pytest.mark.parametrize("err_state", ["err", "warn", "none"])
@pytest.mark.parametrize("tel_loc", ["Center", "Moon", "Earth", "Space"])
@pytest.mark.parametrize("check_frame", ["Moon", "Earth"])
@pytest.mark.parametrize("del_tel_loc", [False, None, True])
def test_check_surface_based_positions(err_state, tel_loc, check_frame, del_tel_loc):
tel_loc_dict = {
"Center": np.array([0, 0, 0]),
"Moon": np.array([0, 0, 1.737e6]),
"Earth": np.array([0, 6.37e6, 0]),
"Space": np.array([4.22e7, 0, 0]),
}
tel_frame_dict = {"Moon": "mcmf", "Earth": "itrs"}
ant_pos = np.array(
[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
)
if del_tel_loc:
ant_pos += tel_loc_dict[tel_loc]
fail_type = err_msg = err_type = None
err_check = check_warnings
if (tel_loc != check_frame) and (err_state != "none"):
if tel_loc == "Center":
fail_type = "below"
elif tel_loc == "Space":
fail_type = "above"
else:
fail_type = "above" if tel_loc == "Earth" else "below"
if fail_type is not None:
err_msg = (
f"{tel_frame_dict[check_frame]} position vector magnitudes must be "
f"on the order of the radius of {check_frame} -- they appear to lie well "
f"{fail_type} this."
)
if err_state == "err":
err_type = ValueError
err_check = pytest.raises
else:
err_type = UserWarning
with err_check(err_type, match=err_msg):
status = utils.coordinates.check_surface_based_positions(
telescope_loc=None if (del_tel_loc) else tel_loc_dict[tel_loc],
antenna_positions=None if (del_tel_loc is None) else ant_pos,
telescope_frame=tel_frame_dict[check_frame],
raise_error=err_state == "err",
raise_warning=err_state == "warn",
)
assert (err_state == "err") or (status == (tel_loc == check_frame))
@pytest.mark.parametrize("tel_loc", ["Earth", "Moon"])
@pytest.mark.parametrize("check_frame", ["Earth", "Moon"])
def test_check_surface_based_positions_earthmoonloc(tel_loc, check_frame):
frame = "mcmf" if (check_frame == "Moon") else "itrs"
if tel_loc == "Earth":
loc = EarthLocation.from_geodetic(0, 0, 0)
else:
pytest.importorskip("lunarsky")
loc = MoonLocation.from_selenodetic(0, 0, 0)
if tel_loc == check_frame:
assert utils.coordinates.check_surface_based_positions(
telescope_loc=loc, telescope_frame=frame
)
else:
with pytest.raises(ValueError, match=(f"{frame} position vector")):
utils.coordinates.check_surface_based_positions(
telescope_loc=[loc.x.value, loc.y.value, loc.z.value],
telescope_frame=frame,
)
Computing file changes ...
