Better handing of discrepancies for LSTs and telescope position in check

CHANGELOG.md

@@ -3,6 +3,20 @@ All notable changes to this project will be documented in this file.
 
 ## [Unreleased]
 
+### Added
+- Added a `check_surface_based_positions` positions method for verifying antenna
+positions are near surface of whatever celestial body their positions are referenced to
+(either the Earth or Moon, currently).
+
+### Changed
+- Increased the tolerance to 75 mas (equivalent to 5 ms time error) for a warning about
+values in `lst_array` not conforming to expectations for `UVData`, `UVCal`, and `UVFlag`
+(was 1 mas) inside of `check`. Additionally, added a keyword to `check` enable the
+tolerance value to be user-specified.
+- Changed the behavior of checking of telescope location to look at the combination of
+`antenna_positions` and `telescope_location` together for `UVData`, `UVCal`, and `UVFlag`.
+Additionally, failing this check results in a warning (was an error).
+
 ## [2.4.1] - 2023-10-13
 
 ### Added

pyuvdata/parameter.py

@@ -895,8 +895,6 @@ class LocationParameter(UVParameter):
 
     """
 
-    _range_dict = {"itrs": (6.35e6, 6.39e6), "mcmf": (1717100.0, 1757100.0)}
-
     def __init__(
         self,
         name,
@@ -905,7 +903,7 @@ def __init__(
         spoof_val=None,
         description="",
         frame="itrs",
-        acceptable_range=-1,
+        acceptable_range=None,
         tols=1e-3,
     ):
         super(LocationParameter, self).__init__(
@@ -921,21 +919,6 @@ def __init__(
         )
         self.frame = frame
 
-        # If acceptable_range is set, set it now. Has to be done after frame is set
-        # because setting the frame changes the acceptable_range.
-        if acceptable_range != -1:
-            self.acceptable_range = acceptable_range
-
-    def __setattr__(self, name, value):
-        """Ensure that acceptable_range is set properly when frame is changed."""
-        if name == "frame":
-            if value in self._range_dict.keys():
-                self.acceptable_range = self._range_dict[value]
-            else:
-                self.acceptable_range = None
-
-        return super().__setattr__(name, value)
-
     def lat_lon_alt(self):
         """Get value in (latitude, longitude, altitude) tuple in radians."""
         if self.value is None:

pyuvdata/tests/test_parameter.py

@@ -397,6 +397,17 @@ def test_location_set_lat_lon_alt_degrees_none():
     assert param1.value is None
 
 
+def test_location_acceptability():
+    """Test check_acceptability with LocationParameters"""
+    val = np.array([0.5, 0.5, 0.5])
+    param1 = uvp.LocationParameter("p1", value=val, acceptable_range=[0, 1])
+    assert param1.check_acceptability()[0]
+
+    val += 0.5
+    param1 = uvp.LocationParameter("p1", value=val, acceptable_range=[0, 1])
+    assert not param1.check_acceptability()[0]
+
+
 def test_location_acceptable_none():
     param1 = uvp.LocationParameter(name="p2", value=1, acceptable_range=None)
 
@@ -408,33 +419,6 @@ def test_location_acceptable_none():
     assert param1.check_acceptability()[0]
 
 
-def test_location_mcmf_acceptability():
-    loc = np.array([1717200.0, 0.0, 0.0])
-    param1 = uvp.LocationParameter(name="p2", value=loc, frame="mcmf")
-    assert param1.acceptable_range == (1717100.0, 1757100.0)
-    assert param1.check_acceptability()[0]
-
-    # test that the acceptable range is changed when the frame is changed
-    param1.frame = "itrs"
-    assert param1.acceptable_range == (6.35e6, 6.39e6)
-    assert not param1.check_acceptability()[0]
-
-    param1.frame = "foo"
-    assert param1.acceptable_range is None
-    assert param1.check_acceptability()[0]
-
-    param1.frame = "mcmf"
-    assert param1.acceptable_range == (1717100.0, 1757100.0)
-    assert param1.check_acceptability()[0]
-
-    # check that if you set the acceptable range on init that it is used
-    param1 = uvp.LocationParameter(
-        name="p2", value=loc, acceptable_range=(1717100.0, 1717150.0), frame="mcmf"
-    )
-    assert param1.acceptable_range == (1717100.0, 1717150.0)
-    assert not param1.check_acceptability()[0]
-
-
 @pytest.mark.parametrize(
     "sky2",
     [

pyuvdata/tests/test_utils.py

@@ -4261,3 +4261,78 @@ def test_uvw_track_generator_moon():
 
     # Check that the total lengths all match 1
     assert np.allclose((gen_results["uvw"] ** 2.0).sum(1), 2.0)
+
+
+@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 = uvtest.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 = uvutils.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.skipif(not hasmoon, reason="lunarsky not installed")
+@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:
+        loc = MoonLocation.from_selenodetic(0, 0, 0)
+
+    if tel_loc == check_frame:
+        assert uvutils.check_surface_based_positions(
+            telescope_loc=loc, telescope_frame=frame
+        )
+    else:
+        with pytest.raises(ValueError, match=(f"{frame} position vector")):
+            uvutils.check_surface_based_positions(
+                telescope_loc=loc, telescope_frame=frame
+            )

pyuvdata/utils.py

@@ -75,6 +75,11 @@
 
 # standard angle tolerance: 1 mas in radians.
 RADIAN_TOL = 1 * 2 * np.pi * 1e-3 / (60.0 * 60.0 * 360.0)
+# standard lst time tolerance: 5 ms (75 mas in radians), based on an expected RMS
+# accuracy of 1 ms at 7 days out from issuance of Bulletin A (which are issued once a
+# week with rapidly determined parameters and forecasted values of DUT1), the exact
+# formula for which is t_err = 0.00025 (MJD-<Bulletin A Release Data>)**0.75 (in secs).
+LST_RAD_TOL = 2 * np.pi * 5e-3 / (86400.0)
 
 # fmt: off
 # polarization constants
@@ -115,6 +120,11 @@
                     "rr": ["r", "r"], "ll": ["l", "l"],
                     "rl": ["r", "l"], "lr": ["l", "r"]}
 
+_range_dict = {
+    "itrs": (6.35e6, 6.39e6, "Earth"), "mcmf": (1717100.0, 1757100.0, "Moon")
+}
+
+
 # fmt: on
 
 
@@ -4078,6 +4088,89 @@ def check_lsts_against_times(
         )
 
 
+def check_surface_based_positions(
+    *,
+    telescope_loc=None,
+    telescope_frame="itrs",
+    antenna_positions=None,
+    raise_error=True,
+    raise_warning=True,
+):
+    """
+    Check that antenna positions are consistent with ground-based values.
+
+    Check that the antenna position, telescope location, or combination of both produces
+    locations that are consistent with surface-based positions. If supplying both
+    antenna position and telescope location, the check will be run against the sum total
+    of both. For the Earth, the permitted range of values is betwen 6350 and 6390 km,
+    whereas for theMoon the range is 1717.1 to 1757.1 km.
+
+    telescope_loc : tuple or EarthLocation or MoonLocation
+        Telescope location, specified as a 3-element tuple (specifying geo/selenocentric
+        position in meters) or as an astropy EarthLocation (or lunarsky MoonLocation).
+    telescope_frame : str, optional
+        Reference frame for latitude/longitude/altitude. Options are itrs (default) or
+        mcmf. Only used if telescope_loc is not an EarthLocation or MoonLocation.
+    antenna_positions : ndarray of float
+        List of antenna positions relative to array center in ECEF coordinates,
+        required if not providing `uvw_array`. Shape is (Nants, 3). If no telescope_loc
+        is specified, these values will be assumed to be relative to geocenter.
+    raise_error : bool
+        If True, an error is raised if telescope_loc and/or telescope_loc do not conform
+        to expectations for a surface-based telescope. Default is True.
+    raise_warning : bool
+        If True, a warning is raised if telescope_loc and/or telescope_loc do not
+        conform to expectations for a surface-based telescope. Default is True, only
+        used if `raise_error` is set to False.
+
+    Returns
+    -------
+    valid : bool
+        If True, the antenna_positions and/or telescope_loc conform to expectations for
+        a surface-based telescope. Otherwise returns false.
+
+    """
+    if antenna_positions is None:
+        antenna_positions = np.zeros((1, 3))
+
+    if isinstance(telescope_loc, EarthLocation) or (
+        hasmoon and isinstance(telescope_loc, MoonLocation)
+    ):
+        antenna_positions = antenna_positions + (
+            telescope_loc.x.to("m").value,
+            telescope_loc.y.to("m").value,
+            telescope_loc.z.to("m").value,
+        )
+    elif telescope_loc is not None:
+        antenna_positions = antenna_positions + telescope_loc
+
+    low_lim, hi_lim, world = _range_dict[telescope_frame]
+
+    err_type = None
+    if np.any(np.sum(antenna_positions**2.0, axis=1) < low_lim**2.0):
+        err_type = "below"
+    elif np.any(np.sum(antenna_positions**2.0, axis=1) > hi_lim**2.0):
+        err_type = "above"
+
+    if err_type is None:
+        return True
+
+    err_msg = (
+        f"{telescope_frame} position vector magnitudes must be on the order of "
+        f"the radius of {world} -- they appear to lie well {err_type} this."
+    )
+
+    # If desired, raise an error
+    if raise_error:
+        raise ValueError(err_msg)
+
+    # Otherwise, if desired, raise a warning instead
+    if raise_warning:
+        warnings.warn(err_msg)
+
+    return False
+
+
 def uvw_track_generator(
     *,
     lon_coord=None,

pyuvdata/uvcal/tests/test_uvcal.py

@@ -3990,13 +3990,7 @@ def test_init_from_uvdata(
     # of precision in the processing pipeline.
     assert uvc_new._time_array == uvc2._time_array
     uvc_new.time_array = uvc2.time_array
-    with uvtest.check_warnings(
-        UserWarning,
-        match="The lst_array is not self-consistent with the time_array and "
-        "telescope location. Consider recomputing with the "
-        "`set_lsts_from_time_array` method.",
-    ):
-        uvc_new.check()
+    uvc_new.check()
 
     uvc_new.set_lsts_from_time_array()
 

pyuvdata/uvcal/uvcal.py

@@ -102,11 +102,7 @@ def __init__(self):
             "telescope_location_lat_lon_alt_degrees properties"
         )
         self._telescope_location = uvp.LocationParameter(
-            "telescope_location",
-            description=desc,
-            acceptable_range=(6.35e6, 6.39e6),
-            tols=1e-3,
-            required=True,
+            "telescope_location", description=desc, tols=1e-3, required=True
         )
 
         desc = (
@@ -1194,7 +1190,11 @@ def _check_freq_spacing(self, raise_errors=True):
         )
 
     def check(
-        self, check_extra=True, run_check_acceptability=True, check_freq_spacing=False
+        self,
+        check_extra=True,
+        run_check_acceptability=True,
+        check_freq_spacing=False,
+        lst_tol=uvutils.LST_RAD_TOL,
     ):
         """
         Add some extra checks on top of checks on UVBase class.
@@ -1212,6 +1212,15 @@ def check(
             Option to check if frequencies are evenly spaced and the spacing is
             equal to their channel_width. This is not required for UVCal
             objects in general but is required to write to calfits files.
+        lst_tol : float or None
+            Tolerance level at which to test LSTs against their expected values. If
+            provided as a float, must be in units of radians. If set to None, the
+            default precision tolerance from the `lst_array` parameter is used (1 mas).
+            Default value is 75 mas,  which is set by the predictive uncertainty in IERS
+            calculations of DUT1 (RMS is of order 1 ms, with with a 5-sigma threshold
+            for detection is used to prevent false issues from being reported), which
+            for some observatories sets the precision with which these values are
+            written.
 
         Returns
         -------
@@ -1309,14 +1318,22 @@ def check(
             self._check_freq_spacing()
 
         if run_check_acceptability:
+            # Check antenna positions
+            uvutils.check_surface_based_positions(
+                antenna_positions=self.antenna_positions,
+                telescope_loc=self.telescope_location,
+                telescope_frame=self._telescope_location.frame,
+                raise_error=False,
+            )
+
             lat, lon, alt = self.telescope_location_lat_lon_alt_degrees
             uvutils.check_lsts_against_times(
                 jd_array=self.time_array,
                 lst_array=self.lst_array,
                 latitude=lat,
                 longitude=lon,
                 altitude=alt,
-                lst_tols=self._lst_array.tols,
+                lst_tols=self._lst_array.tols if lst_tol is None else [0, lst_tol],
                 frame=self._telescope_location.frame,
             )