From a4d89efe6f1b2994cd83e1ca932ab2732148629c Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Sat, 5 Aug 2023 18:06:02 +0200 Subject: [PATCH 01/12] Implement support for targets defined by TLE (Two-Line Element set) * Add TLETarget class that uses skyfield internally to return RA/Dec values for specific times * The get_skycoord() function now takes a 'time' argument because the returned SkyCoord RA/Dec values are time-dependent for TLETargets * get_skycoord() behavior for non-TLETargets and multiple times can be controlled with the 'backwards_compatible' argument * Only perform 'grid_times_targets' when targets and times are not broadcastable * Update the schedule table so it also displays sensible info for TLETargets * Update Tests --- astroplan/constraints.py | 13 +- astroplan/exceptions.py | 4 + astroplan/observer.py | 4 +- astroplan/scheduling.py | 42 +++-- astroplan/target.py | 250 +++++++++++++++++++++++++++- astroplan/tests/test_constraints.py | 11 +- astroplan/tests/test_scheduling.py | 48 +++++- astroplan/tests/test_target.py | 128 +++++++++++++- pyproject.toml | 6 +- 9 files changed, 468 insertions(+), 38 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 17b74f84..3be2fdeb 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -254,16 +254,17 @@ def __call__(self, observer, targets, times=None, time_resolution=time_grid_resolution) if grid_times_targets: - targets = get_skycoord(targets) + targets = get_skycoord(targets, times) # TODO: these broadcasting operations are relatively slow # but there is potential for huge speedup if the end user # disables gridding and re-shapes the coords themselves # prior to evaluating multiple constraints. - if targets.isscalar: - # ensure we have a (1, 1) shape coord - targets = SkyCoord(np.tile(targets, 1))[:, np.newaxis] - else: - targets = targets[..., np.newaxis] + if not observer._is_broadcastable(targets.shape, times.shape): + if targets.isscalar: + # ensure we have a (1, 1) shape coord + targets = SkyCoord(np.tile(targets, 1))[:, np.newaxis] + else: + targets = targets[..., np.newaxis] times, targets = observer._preprocess_inputs(times, targets, grid_times_targets=False) result = self.compute_constraint(times, observer, targets) diff --git a/astroplan/exceptions.py b/astroplan/exceptions.py index 5221b367..61bad2c0 100644 --- a/astroplan/exceptions.py +++ b/astroplan/exceptions.py @@ -40,3 +40,7 @@ class PlotBelowHorizonWarning(PlotWarning): class MissingConstraintWarning(AstroplanWarning): """Triggered when a constraint is expected but not supplied""" pass + +class InvalidTLEDataWarning(AstroplanWarning): + """TLE data invalid for the requested time""" + pass diff --git a/astroplan/observer.py b/astroplan/observer.py index a86b4793..def4d2cc 100644 --- a/astroplan/observer.py +++ b/astroplan/observer.py @@ -515,8 +515,8 @@ def _preprocess_inputs(self, time, target=None, grid_times_targets=False): return time, None # convert any kind of target argument to non-scalar SkyCoord - target = get_skycoord(target) - if grid_times_targets: + target = get_skycoord(target, time) + if grid_times_targets and not self._is_broadcastable(target.shape, time.shape): if target.isscalar: # ensure we have a (1, 1) shape coord target = SkyCoord(np.tile(target, 1))[:, np.newaxis] diff --git a/astroplan/scheduling.py b/astroplan/scheduling.py index 8fda4a2b..debce628 100644 --- a/astroplan/scheduling.py +++ b/astroplan/scheduling.py @@ -14,7 +14,7 @@ from .utils import time_grid_from_range, stride_array from .constraints import AltitudeConstraint -from .target import get_skycoord +from .target import get_skycoord, FixedTarget, TLETarget __all__ = ['ObservingBlock', 'TransitionBlock', 'Schedule', 'Slot', 'Scheduler', 'SequentialScheduler', 'PriorityScheduler', @@ -118,7 +118,6 @@ def __init__(self, blocks, observer, schedule, global_constraints=[]): self.observer = observer self.schedule = schedule self.global_constraints = global_constraints - self.targets = get_skycoord([block.target for block in self.blocks]) def create_score_array(self, time_resolution=1*u.minute): """ @@ -139,6 +138,7 @@ def create_score_array(self, time_resolution=1*u.minute): start = self.schedule.start_time end = self.schedule.end_time times = time_grid_from_range((start, end), time_resolution) + targets = get_skycoord([block.target for block in self.blocks], times) score_array = np.ones((len(self.blocks), len(times))) for i, block in enumerate(self.blocks): # TODO: change the default constraints from None to [] @@ -148,7 +148,7 @@ def create_score_array(self, time_resolution=1*u.minute): times=times) score_array[i] *= applied_score for constraint in self.global_constraints: - score_array *= constraint(self.observer, self.targets, times, + score_array *= constraint(self.observer, targets, times, grid_times_targets=True) return score_array @@ -270,25 +270,32 @@ def to_table(self, show_transitions=True, show_unused=False): start_times = [] end_times = [] durations = [] - ra = [] - dec = [] + coordiante_type = [] + coordinate_info = [] config = [] for slot in self.slots: if hasattr(slot.block, 'target'): start_times.append(slot.start.iso) end_times.append(slot.end.iso) - durations.append(slot.duration.to(u.minute).value) + durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) target_names.append(slot.block.target.name) - ra.append(u.Quantity(slot.block.target.ra)) - dec.append(u.Quantity(slot.block.target.dec)) + if isinstance(slot.block.target, FixedTarget): + coordiante_type.append("RA/Dec") + coordinate_info.append(slot.block.target.coord.to_string('hmsdms')) + elif isinstance(slot.block.target, TLETarget): + coordiante_type.append("TLE") + coordinate_info.append( + f"#{slot.block.target.satellite.model.satnum} " + f"epoch {slot.block.target.satellite.epoch.utc_strftime(format='%Y-%m-%d %H:%M:%S')}" + ) config.append(slot.block.configuration) elif show_transitions and slot.block: start_times.append(slot.start.iso) end_times.append(slot.end.iso) - durations.append(slot.duration.to(u.minute).value) + durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) target_names.append('TransitionBlock') - ra.append('') - dec.append('') + coordiante_type.append('') + coordinate_info.append('') changes = list(slot.block.components.keys()) if 'slew_time' in changes: changes.remove('slew_time') @@ -296,14 +303,14 @@ def to_table(self, show_transitions=True, show_unused=False): elif slot.block is None and show_unused: start_times.append(slot.start.iso) end_times.append(slot.end.iso) - durations.append(slot.duration.to(u.minute).value) + durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) target_names.append('Unused Time') - ra.append('') - dec.append('') + coordiante_type.append('') + coordinate_info.append('') config.append('') - return Table([target_names, start_times, end_times, durations, ra, dec, config], + return Table([target_names, start_times, end_times, durations, coordiante_type, coordinate_info, config], names=('target', 'start time (UTC)', 'end time (UTC)', - 'duration (minutes)', 'ra', 'dec', 'configuration')) + 'duration (min)', 'type', 'coordinates/info', 'configuration')) def new_slots(self, slot_index, start_time, end_time): """ @@ -996,9 +1003,8 @@ def __call__(self, oldblock, newblock, start_time, observer): # use the constraints cache for now, but should move that machinery # to observer from .constraints import _get_altaz - from .target import get_skycoord if oldblock.target != newblock.target: - targets = get_skycoord([oldblock.target, newblock.target]) + targets = get_skycoord([oldblock.target, newblock.target], start_time) aaz = _get_altaz(start_time, observer, targets)['altaz'] sep = aaz[0].separation(aaz[1]) if sep/self.slew_rate > 1 * u.second: diff --git a/astroplan/target.py b/astroplan/target.py index e82c7e5e..cf2a9165 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -2,12 +2,27 @@ # Standard library from abc import ABCMeta +import warnings # Third-party import astropy.units as u -from astropy.coordinates import SkyCoord, ICRS, UnitSphericalRepresentation +from astropy.time import Time +from astropy.coordinates import SkyCoord, ICRS, UnitSphericalRepresentation, AltAz, EarthLocation +import numpy as np +try: + from sgp4.io import twoline2rv + from sgp4.earth_gravity import wgs84 as sgp4_wgs84 + from skyfield.api import load, wgs84 + from skyfield.sgp4lib import EarthSatellite + skyfield_available = True +except ImportError: + skyfield_available = False -__all__ = ["Target", "FixedTarget", "NonFixedTarget"] +# Package +from .exceptions import InvalidTLEDataWarning + + +__all__ = ["Target", "FixedTarget", "NonFixedTarget", "TLETarget"] # Docstring code examples include printed SkyCoords, but the format changed # in astropy 1.3. Thus the doctest needs astropy >=1.3 and this is the @@ -188,31 +203,250 @@ class NonFixedTarget(Target): """ -def get_skycoord(targets): +class TLETarget(Target): + """ + A target defined by TLE (Two-Line Element set) for satellites. + """ + def __init__(self, line1, line2, name=None, observer=None, skip_tle_check=False): + """ + Parameters + ---------- + line1 : str + The first line of the TLE set + + line2 : str + The second line of the TLE set + + name : str, optional + Name of the target, used for plotting and representing the target + as a string + + observer : `~astropy.coordinates.EarthLocation`, `~astroplan.Observer`, optional + The location of observer. If `None`, the observer is assumed to be at sea level at the equator. + + skip_tle_check : bool, optional + Whether to skip TLE validation + """ + if not skyfield_available: + raise ImportError("Please install the skyfield package to use the TLETarget class.") + + if not skip_tle_check: + twoline2rv(line1, line2, sgp4_wgs84) # Raises ValueError if TLE is invalid + + self.name = name + self.satellite = EarthSatellite(line1, line2, name, load.timescale()) + + if observer is None: + self.observer = wgs84.latlon(0, 0, 0) + else: + observer = getattr(observer, 'location', observer) + longitude, latitude, height = observer.to_geodetic() + self.observer = wgs84.latlon(latitude.to(u.deg).value, + longitude.to(u.deg).value, + height.to(u.m).value) + + @classmethod + def from_string(cls, tle_string, name=None, *args, **kwargs): + """ + Creates a TLETarget instance from a complete TLE string. + + Parameters + ---------- + tle_string : str + String to be parsed, expected to contain 2 or 3 newline-separated lines. + + name : str, optional + Name of the target. If not provided and the tle_string contains 3 lines, + the first line will be used as the name. + + args, kwargs : tuple, dict, optional + Additional arguments and keyword arguments to be passed to the TLETarget class constructor. + """ + lines = tle_string.strip().splitlines() + + if len(lines) not in (2, 3): + raise ValueError(f"Expected TLE string to contain 2 or 3 lines, got {len(lines)}") + + if len(lines) == 3: + line1, line2, name = lines[1], lines[2], name or lines[0] + else: # len(lines) == 2 + line1, line2 = lines + return cls(line1, line2, name, *args, **kwargs) + + @property + def ra(self): + raise NotImplementedError("Satellite RA changes rapidly, compute it at a specific time with self.coord(time)") + + @property + def dec(self): + raise NotImplementedError("Satellite Dec changes rapidly, compute it at a specific time with self.coord(time)") + + def _compute_topocentric(self, time=None): + """ + Compute the topocentric coordinates (relative to observer) at a particular time. + + Parameters + ---------- + time : `~astropy.time.Time`, optional + The time(s) to use in the calculation. + + Returns + ------- + topocentric : `skyfield.positionlib.ICRF` + The topocentric object representing the relative coordinates of the target. + """ + if time is None: + time = Time.now() + ts = load.timescale() + t = ts.from_astropy(time) + + topocentric = (self.satellite - self.observer).at(t) + + # Check for invalid TLE data. A non-None usually message means the computation went beyond the physically + # sensible point. Details: https://rhodesmill.org/skyfield/earth-satellites.html#detecting-propagation-errors + if ((topocentric.message is not None and not isinstance(topocentric.message, list)) or + (isinstance(topocentric.message, list) and not all(x == None for x in topocentric.message))): + warnings.warn(f"Invalid TLE Data: {topocentric.message}", InvalidTLEDataWarning) + return topocentric + + def coord(self, time=None): + """ + Get the coordinates of the target at a particular time. + + Parameters + ---------- + time : `~astropy.time.Time`, optional + The time(s) to use in the calculation. + + Returns + ------- + coord : `~astropy.coordinates.SkyCoord` + A single SkyCoord object, which may be non-scalar, representing the target's + RA/Dec coordinates at the specified time(s). Might return np.nan and output a + warning for times where the elements stop making physical sense. + """ + topocentric = self._compute_topocentric(time) + ra, dec, distance = topocentric.radec() + # Don't add distance to SkyCoord: in SkyCoord, distance is from frame origin, but here, it's from observer. + return SkyCoord(ra.hours*u.hourangle, dec.degrees*u.deg, obstime=time, frame='icrs') + + def altaz(self, time=None): + """ + Get the altitude and azimuth of the target at a particular time. + + Parameters + ---------- + time : `~astropy.time.Time`, optional + The time(s) to use in the calculation. + + Returns + ------- + altaz_coord : `~astropy.coordinates.SkyCoord` + A SkyCoord object representing the target's altitude and azimuth at the specified time(s). + """ + topocentric = self._compute_topocentric(time) + alt, az, distance = topocentric.altaz() + + earth_location = EarthLocation(lat=self.observer.latitude.degrees*u.deg, + lon=self.observer.longitude.degrees*u.deg, + height=self.observer.elevation.m*u.m) + + altaz = AltAz(alt=alt.degrees*u.deg, az=az.degrees*u.deg, obstime=time, location=earth_location) + return SkyCoord(altaz) + + def __repr__(self): + return f'<{self.__class__.__name__} "{self.name}">' + + def __str__(self): + return self.name + + + +def repeat_skycoord(coord, times): + """ + Repeats the coordinates of a SkyCoord object 'times.size' number of times. + + Parameters + ---------- + coord : `~astropy.coordinates.SkyCoord` + The original SkyCoord object whose coordinates need to be repeated. + + times : `~astropy.time.Time` + The size of times determines the number of times the coordinates should be repeated. + + Returns + -------- + SkyCoord : `~astropy.coordinates.SkyCoord` + A new SkyCoord object with the coordinates of the original object + repeated 'times.size' number of times. If the SkyCoord object is scalar + or 'times' is None or a scalar, this function returns the + original SkyCoord object. + """ + if coord.size != 1 or times is None or times.size == 1: + return coord + return SkyCoord( + ra=coord.ra.repeat(times.size), + dec=coord.dec.repeat(times.size), + distance=None if coord.distance.unit is u.one else coord.distance.repeat(times.size), + frame=coord.frame, + obstime=times + ) + +def get_skycoord(targets, time=None, backwards_compatible=True): """ Return an `~astropy.coordinates.SkyCoord` object. When performing calculations it is usually most efficient to have a single `~astropy.coordinates.SkyCoord` object, rather than a - list of `FixedTarget` or `~astropy.coordinates.SkyCoord` objects. + list of `Target` or `~astropy.coordinates.SkyCoord` objects. This is a convenience routine to do that. Parameters ----------- - targets : list, `~astropy.coordinates.SkyCoord`, `Fixedtarget` + targets : list, `~astropy.coordinates.SkyCoord`, `Target` either a single target or a list of targets + time : `~astropy.time.Time`, optional + The time(s) to use in the calculation. + + backwards_compatible : bool, optional + Controls the output format when only FixedTarget or SkyCoord targets are combined with a time argument. + If False, it will return (targets x times), where all coordinates per target are the same. + If True, it will return one coordinate per target (default is True). + Returns -------- coord : `~astropy.coordinates.SkyCoord` a single SkyCoord object, which may be non-scalar """ + + # Note on backwards_compatible: + # This method will always return (targets x times) when there is a TLETarget in targets, because RA/Dec changes with time. + # Do we want to be 100% backwards compatible, or do we prefer consistent output, + # for FixedTarget or SkyCoord targets combined with multiple times? + # backwards_compatible = True will continue to return one coordinate per target + # backwards_compatible = False will now return (targets x times), where all coordinates per target are the same + + # Early exit for single target if not isinstance(targets, list): - return getattr(targets, 'coord', targets) + if isinstance(targets, TLETarget): + return targets.coord(time) + else: + if backwards_compatible: + return getattr(targets, 'coord', targets) + else: + return repeat_skycoord(getattr(targets, 'coord', targets), time) + + # Identify if any of the targets is not FixedTarget or SkyCoord + has_non_fixed_target = any(not isinstance(target, (FixedTarget, SkyCoord)) for target in targets) + + # Get the SkyCoord object itself + coords = [target.coord(time) if isinstance(target, TLETarget) else getattr(target, 'coord', target) for target in targets] - # get the SkyCoord object itself - coords = [getattr(target, 'coord', target) for target in targets] + # Fill time dimension for SkyCoords that only have a single coordinate + if (backwards_compatible and has_non_fixed_target or not backwards_compatible) and time is not None: + coords = [repeat_skycoord(coord, time) for coord in coords] # are all SkyCoordinate's in equivalent frames? If not, convert to ICRS convert_to_icrs = not all( diff --git a/astroplan/tests/test_constraints.py b/astroplan/tests/test_constraints.py index e4db748a..c0915c30 100644 --- a/astroplan/tests/test_constraints.py +++ b/astroplan/tests/test_constraints.py @@ -461,7 +461,16 @@ def test_caches_shapes(): targets = get_skycoord([m31, ippeg, htcas]) observer = Observer.at_site('lapalma') ac = AltitudeConstraint(min=30*u.deg) - assert ac(observer, targets, times, grid_times_targets=True).shape == (3, 3) + + # When time and targets have the same size, they are broadcastable so grid_times_targets is ignored + assert ac(observer, targets, times, grid_times_targets=True).shape == (3,) + targets = get_skycoord([m31, ippeg]) + assert ac(observer, targets, times, grid_times_targets=True).shape == (2, 3) + + # When time and targets don't have the same size this fails with grid_times_targets=False + with pytest.raises(ValueError): + ac(observer, targets, times, grid_times_targets=False) + targets = get_skycoord([m31, ippeg, htcas]) assert ac(observer, targets, times, grid_times_targets=False).shape == (3,) diff --git a/astroplan/tests/test_scheduling.py b/astroplan/tests/test_scheduling.py index 5bfcd874..25bd3b22 100644 --- a/astroplan/tests/test_scheduling.py +++ b/astroplan/tests/test_scheduling.py @@ -5,15 +5,21 @@ import pytest from astropy.coordinates import SkyCoord, EarthLocation from astropy.time import Time +try: + import skyfield # noqa + HAS_SKYFIELD = True +except ImportError: + HAS_SKYFIELD = False from astroplan.utils import time_grid_from_range from astroplan.observer import Observer -from astroplan.target import FixedTarget, get_skycoord +from astroplan.target import FixedTarget, TLETarget, get_skycoord from astroplan.constraints import (AirmassConstraint, AtNightConstraint, _get_altaz, MoonIlluminationConstraint, PhaseConstraint) from astroplan.periodic import EclipsingSystem from astroplan.scheduling import (ObservingBlock, PriorityScheduler, SequentialScheduler, Transitioner, TransitionBlock, Schedule, Slot, Scorer) +from astroplan.exceptions import InvalidTLEDataWarning vega = FixedTarget(coord=SkyCoord(ra=279.23473479 * u.deg, dec=38.78368896 * u.deg), name="Vega") @@ -329,3 +335,43 @@ def test_scorer(): scores = scorer.create_score_array(time_resolution=20 * u.minute) # the ``global_constraint``: constraint2 should have applied to the blocks assert np.array_equal(c2, scores) + + +@pytest.mark.skipif('not HAS_SKYFIELD') +def test_priority_scheduler_TLETarget(): + line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" + iss = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=apo) + constraints = [AirmassConstraint(3, boolean_constraint=False)] + blocks = [ObservingBlock(t, 5*u.minute, i) for i, t in enumerate(targets)] + blocks.append(ObservingBlock(iss, 0.5*u.minute, 4)) + start_time = Time('2016-02-06 03:00:00') + end_time = start_time + 1*u.hour + scheduler = PriorityScheduler(transitioner=default_transitioner, + constraints=constraints, observer=apo, + time_resolution=0.5*u.minute) + schedule = Schedule(start_time, end_time) + scheduler(blocks, schedule) + assert len(schedule.observing_blocks) == 3 + assert all(np.abs(block.end_time - block.start_time - block.duration) < + 1*u.second for block in schedule.scheduled_blocks) + assert all([schedule.observing_blocks[0].target == polaris, + schedule.observing_blocks[1].target == rigel, + schedule.observing_blocks[2].target == iss]) + + # test that the scheduler does not error when called with a partially + # filled schedule + scheduler(blocks, schedule) + scheduler(blocks, schedule) + + # Time too far in the future where elements stop making physical sense + with pytest.warns(): + start_time = Time('2035-08-02 10:00:00') + end_time = start_time + 1*u.hour + schedule = Schedule(start_time, end_time) + with pytest.warns(InvalidTLEDataWarning): + scheduler(blocks, schedule) + assert len(schedule.observing_blocks) == 3 + assert all([schedule.observing_blocks[0].target == vega, + schedule.observing_blocks[1].target == rigel, + schedule.observing_blocks[2].target == polaris]) diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index ba2dde47..ca77bdc9 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -5,10 +5,18 @@ import pytest from astropy.coordinates import SkyCoord, GCRS, ICRS from astropy.time import Time +import numpy as np +try: + import skyfield # noqa + HAS_SKYFIELD = True +except ImportError: + HAS_SKYFIELD = False # Package -from astroplan.target import FixedTarget, get_skycoord +from astroplan.target import FixedTarget, TLETarget, get_skycoord from astroplan.observer import Observer +from astroplan.utils import time_grid_from_range +from astroplan.exceptions import InvalidTLEDataWarning @pytest.mark.remote_data @@ -44,6 +52,79 @@ def test_FixedTarget_ra_dec(): assert vega.coord.dec == vega_coords.dec == vega.dec, ('Retrieve Dec from ' 'SkyCoord') +@pytest.mark.skipif('not HAS_SKYFIELD') +def test_TLETarget(): + tle_string=("ISS (ZARYA)\n" + "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" + "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") + line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" + subaru = Observer.at_site('subaru') + time = Time("2023-08-02 10:00", scale='utc') + times = time_grid_from_range([time, time + 3.1*u.hour], + time_resolution=1 * u.hour) + + tle_target1 = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=subaru) + tle_target2 = TLETarget.from_string(tle_string=tle_string, observer=subaru) + + assert tle_target1.name == "ISS (ZARYA)" + assert tle_target2.name == "ISS (ZARYA)" + assert repr(tle_target1) == repr(tle_target2) + assert str(tle_target1) == str(tle_target2) + + # Single time + ra_dec1 = tle_target1.coord(time) + ra_dec2 = tle_target2.coord(time) + ra_dec_from_horizon = SkyCoord("08h29m26.36s +07d31m24.0s") + + assert ra_dec1.separation(ra_dec_from_horizon) < 30*u.arcsec + assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') + + # Multiple times + ra_dec1 = tle_target1.coord(times) + ra_dec2 = tle_target2.coord(times) + ra_dec_from_horizon = SkyCoord(["08h29m26.36s +07d31m24.0s", + "06h25m46.81s -54d32m20.2s", + "13h52m09.08s +04d26m40.1s", + "09h20m04.74s -00d51m19.2s"]) + + assert all(list(ra_dec1.separation(ra_dec_from_horizon) < 30*u.arcsec)) + assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') + + # TLE Check + line1_invalid ="1 25544U 98067A .00041610 00000-0 73103-3 0 9990" + with pytest.raises(ValueError): + TLETarget(name="ISS (ZARYA)", line1=line1_invalid, line2=line2, observer=subaru) + TLETarget(name="ISS (ZARYA)", line1=line1_invalid, line2=line2, observer=subaru, skip_tle_check=True) + + # from_string + tle_string=("1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" + "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") + tle_target3 = TLETarget.from_string(tle_string=tle_string, observer=subaru, name="ISS") + assert tle_target3.name == "ISS" + + tle_string = "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + with pytest.raises(ValueError): + TLETarget.from_string(tle_string=tle_string, observer=subaru) + + + # AltAz (This is slow) + altaz_observer = subaru.altaz(time, tle_target1) + altaz_skyfield = tle_target1.altaz(time) + + assert altaz_observer.separation(altaz_skyfield) < 20*u.arcsec + + # Time too far in the future where elements stop making physical sense + with pytest.warns(): # ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6) + time_invalid = Time("2035-08-02 10:00", scale='utc') + times_list = list(times) + times_list[2] = Time("2035-08-02 10:00", scale='utc') + times_invalid = Time(times_list) + with pytest.warns(InvalidTLEDataWarning): + assert np.isnan(tle_target1.coord(time_invalid).ra) + with pytest.warns(InvalidTLEDataWarning): + assert np.isnan(tle_target1.coord(times_invalid)[2].ra) + @pytest.mark.remote_data def test_get_skycoord(): @@ -80,3 +161,48 @@ def test_get_skycoord(): coo = get_skycoord([m31_gcrs, m31_gcrs_with_distance]) assert coo.is_equivalent_frame(m31_gcrs.frame) assert len(coo) == 2 + + +@pytest.mark.skipif('not HAS_SKYFIELD') +def test_get_skycoord_with_TLETarget(): + skycoord_targed = SkyCoord(10.6847083*u.deg, 41.26875*u.deg) + fixed_target1 = FixedTarget(name="fixed1", coord=SkyCoord(279.23458, 38.78369, unit='deg')) + line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" + subaru = Observer.at_site('subaru') + tle_target1 = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=subaru) + + time = Time("2023-08-02 10:00", scale='utc') + times = time_grid_from_range([time, time + 3.1*u.hour], + time_resolution=1 * u.hour) + + tle_output = get_skycoord(tle_target1, time) + assert tle_output.size == 1 + tle_output = get_skycoord(tle_target1, times) + assert tle_output.shape == (4,) + tle_output = get_skycoord([tle_target1, tle_target1], time) + assert tle_output.shape == (2,) + tle_output = get_skycoord([tle_target1, tle_target1], times) + assert tle_output.shape == (2,4) + + mixed_output = get_skycoord([skycoord_targed, fixed_target1, tle_target1], time) + assert mixed_output.shape == (3,) + mixed_output = get_skycoord([skycoord_targed, fixed_target1, tle_target1], times) + assert mixed_output.shape == (3,4) + + # backwards_compatible + tle_output = get_skycoord(skycoord_targed, time, backwards_compatible=False) + assert tle_output.size == 1 + + tle_output = get_skycoord(skycoord_targed, times) + assert tle_output.size == 1 + tle_output = get_skycoord(skycoord_targed, times, backwards_compatible=False) + assert tle_output.shape == (4,) + + tle_output = get_skycoord([skycoord_targed, skycoord_targed], time, backwards_compatible=False) + assert tle_output.shape == (2,) + + tle_output = get_skycoord([skycoord_targed, skycoord_targed], times) + assert tle_output.shape == (2,) + tle_output = get_skycoord([skycoord_targed, skycoord_targed], times, backwards_compatible=False) + assert tle_output.shape == (2,4) diff --git a/pyproject.toml b/pyproject.toml index c8ad7c76..e4110c49 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,7 +34,8 @@ dependencies = [ [project.optional-dependencies] all = [ "matplotlib>=3.6.0", - "astroquery" + "astroquery", + "skyfield" ] test = [ "pytest-astropy", @@ -48,6 +49,9 @@ docs = [ plotting = [ "astroplan[all]" ] +tle = [ + "skyfield" +] [project.readme] file = "README.rst" From a9339543fe3eaddd038d65c412c0fe2514e93d67 Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Sat, 5 Aug 2023 19:00:51 +0200 Subject: [PATCH 02/12] Fix code style checks --- astroplan/exceptions.py | 1 + astroplan/scheduling.py | 16 ++++---- astroplan/target.py | 58 ++++++++++++++++++----------- astroplan/tests/test_constraints.py | 3 +- astroplan/tests/test_scheduling.py | 4 +- astroplan/tests/test_target.py | 35 ++++++++--------- 6 files changed, 69 insertions(+), 48 deletions(-) diff --git a/astroplan/exceptions.py b/astroplan/exceptions.py index 61bad2c0..18642cce 100644 --- a/astroplan/exceptions.py +++ b/astroplan/exceptions.py @@ -41,6 +41,7 @@ class MissingConstraintWarning(AstroplanWarning): """Triggered when a constraint is expected but not supplied""" pass + class InvalidTLEDataWarning(AstroplanWarning): """TLE data invalid for the requested time""" pass diff --git a/astroplan/scheduling.py b/astroplan/scheduling.py index debce628..867fdb88 100644 --- a/astroplan/scheduling.py +++ b/astroplan/scheduling.py @@ -278,15 +278,16 @@ def to_table(self, show_transitions=True, show_unused=False): start_times.append(slot.start.iso) end_times.append(slot.end.iso) durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) - target_names.append(slot.block.target.name) - if isinstance(slot.block.target, FixedTarget): + target = slot.block.target + target_names.append(target.name) + if isinstance(target, FixedTarget): coordiante_type.append("RA/Dec") - coordinate_info.append(slot.block.target.coord.to_string('hmsdms')) - elif isinstance(slot.block.target, TLETarget): + coordinate_info.append(target.coord.to_string('hmsdms')) + elif isinstance(target, TLETarget): coordiante_type.append("TLE") coordinate_info.append( - f"#{slot.block.target.satellite.model.satnum} " - f"epoch {slot.block.target.satellite.epoch.utc_strftime(format='%Y-%m-%d %H:%M:%S')}" + f"#{target.satellite.model.satnum} " + f"epoch {target.satellite.epoch.utc_strftime(format='%Y-%m-%d %H:%M:%S')}" ) config.append(slot.block.configuration) elif show_transitions and slot.block: @@ -308,7 +309,8 @@ def to_table(self, show_transitions=True, show_unused=False): coordiante_type.append('') coordinate_info.append('') config.append('') - return Table([target_names, start_times, end_times, durations, coordiante_type, coordinate_info, config], + return Table([target_names, start_times, end_times, + durations, coordiante_type, coordinate_info, config], names=('target', 'start time (UTC)', 'end time (UTC)', 'duration (min)', 'type', 'coordinates/info', 'configuration')) diff --git a/astroplan/target.py b/astroplan/target.py index cf2a9165..8a020ab1 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -8,7 +8,6 @@ import astropy.units as u from astropy.time import Time from astropy.coordinates import SkyCoord, ICRS, UnitSphericalRepresentation, AltAz, EarthLocation -import numpy as np try: from sgp4.io import twoline2rv from sgp4.earth_gravity import wgs84 as sgp4_wgs84 @@ -222,7 +221,8 @@ def __init__(self, line1, line2, name=None, observer=None, skip_tle_check=False) as a string observer : `~astropy.coordinates.EarthLocation`, `~astroplan.Observer`, optional - The location of observer. If `None`, the observer is assumed to be at sea level at the equator. + The location of observer. + If `None`, the observer is assumed to be at sea level at the equator. skip_tle_check : bool, optional Whether to skip TLE validation @@ -260,7 +260,7 @@ def from_string(cls, tle_string, name=None, *args, **kwargs): the first line will be used as the name. args, kwargs : tuple, dict, optional - Additional arguments and keyword arguments to be passed to the TLETarget class constructor. + Additional arguments and keyword arguments to be passed to the TLETarget constructor. """ lines = tle_string.strip().splitlines() @@ -269,17 +269,17 @@ def from_string(cls, tle_string, name=None, *args, **kwargs): if len(lines) == 3: line1, line2, name = lines[1], lines[2], name or lines[0] - else: # len(lines) == 2 + else: # len(lines) == 2 line1, line2 = lines return cls(line1, line2, name, *args, **kwargs) @property def ra(self): - raise NotImplementedError("Satellite RA changes rapidly, compute it at a specific time with self.coord(time)") + raise NotImplementedError("Compute satellite RA at a specific time with self.coord(time)") @property def dec(self): - raise NotImplementedError("Satellite Dec changes rapidly, compute it at a specific time with self.coord(time)") + raise NotImplementedError("Compute satellite RA at a specific time with self.coord(time)") def _compute_topocentric(self, time=None): """ @@ -302,11 +302,16 @@ def _compute_topocentric(self, time=None): topocentric = (self.satellite - self.observer).at(t) - # Check for invalid TLE data. A non-None usually message means the computation went beyond the physically - # sensible point. Details: https://rhodesmill.org/skyfield/earth-satellites.html#detecting-propagation-errors - if ((topocentric.message is not None and not isinstance(topocentric.message, list)) or - (isinstance(topocentric.message, list) and not all(x == None for x in topocentric.message))): - warnings.warn(f"Invalid TLE Data: {topocentric.message}", InvalidTLEDataWarning) + # Check for invalid TLE data. A non-None usually message means the computation went beyond + # the physically sensible point. Details: + # https://rhodesmill.org/skyfield/earth-satellites.html#detecting-propagation-errors + message = topocentric.message + if ( + (message is not None and not isinstance(message, list)) or + (isinstance(message, list) and not all(x is None for x in message)) + ): + warnings.warn(f"Invalid TLE Data: {message}", InvalidTLEDataWarning) + return topocentric def coord(self, time=None): @@ -327,7 +332,7 @@ def coord(self, time=None): """ topocentric = self._compute_topocentric(time) ra, dec, distance = topocentric.radec() - # Don't add distance to SkyCoord: in SkyCoord, distance is from frame origin, but here, it's from observer. + # No distance, in SkyCoord, distance is from frame origin, but here, it's from observer. return SkyCoord(ra.hours*u.hourangle, dec.degrees*u.deg, obstime=time, frame='icrs') def altaz(self, time=None): @@ -342,7 +347,7 @@ def altaz(self, time=None): Returns ------- altaz_coord : `~astropy.coordinates.SkyCoord` - A SkyCoord object representing the target's altitude and azimuth at the specified time(s). + SkyCoord object representing the target's altitude and azimuth at the specified time(s) """ topocentric = self._compute_topocentric(time) alt, az, distance = topocentric.altaz() @@ -351,7 +356,8 @@ def altaz(self, time=None): lon=self.observer.longitude.degrees*u.deg, height=self.observer.elevation.m*u.m) - altaz = AltAz(alt=alt.degrees*u.deg, az=az.degrees*u.deg, obstime=time, location=earth_location) + altaz = AltAz(alt=alt.degrees*u.deg, az=az.degrees*u.deg, + obstime=time, location=earth_location) return SkyCoord(altaz) def __repr__(self): @@ -361,7 +367,6 @@ def __str__(self): return self.name - def repeat_skycoord(coord, times): """ Repeats the coordinates of a SkyCoord object 'times.size' number of times. @@ -392,6 +397,7 @@ def repeat_skycoord(coord, times): obstime=times ) + def get_skycoord(targets, time=None, backwards_compatible=True): """ Return an `~astropy.coordinates.SkyCoord` object. @@ -411,7 +417,7 @@ def get_skycoord(targets, time=None, backwards_compatible=True): The time(s) to use in the calculation. backwards_compatible : bool, optional - Controls the output format when only FixedTarget or SkyCoord targets are combined with a time argument. + Controls output format when FixedTarget or SkyCoord targets are used with a time argument. If False, it will return (targets x times), where all coordinates per target are the same. If True, it will return one coordinate per target (default is True). @@ -422,11 +428,11 @@ def get_skycoord(targets, time=None, backwards_compatible=True): """ # Note on backwards_compatible: - # This method will always return (targets x times) when there is a TLETarget in targets, because RA/Dec changes with time. + # Method always returns (targets x times) with TLETarget in targets, as RA/Dec changes with time # Do we want to be 100% backwards compatible, or do we prefer consistent output, # for FixedTarget or SkyCoord targets combined with multiple times? # backwards_compatible = True will continue to return one coordinate per target - # backwards_compatible = False will now return (targets x times), where all coordinates per target are the same + # backwards_compatible = False, returns (targets x times) with identical coordinates per target # Early exit for single target if not isinstance(targets, list): @@ -439,13 +445,23 @@ def get_skycoord(targets, time=None, backwards_compatible=True): return repeat_skycoord(getattr(targets, 'coord', targets), time) # Identify if any of the targets is not FixedTarget or SkyCoord - has_non_fixed_target = any(not isinstance(target, (FixedTarget, SkyCoord)) for target in targets) + has_non_fixed_target = any( + not isinstance(target, (FixedTarget, SkyCoord)) + for target in targets + ) # Get the SkyCoord object itself - coords = [target.coord(time) if isinstance(target, TLETarget) else getattr(target, 'coord', target) for target in targets] + coords = [ + target.coord(time) if isinstance(target, TLETarget) + else getattr(target, 'coord', target) + for target in targets + ] # Fill time dimension for SkyCoords that only have a single coordinate - if (backwards_compatible and has_non_fixed_target or not backwards_compatible) and time is not None: + if ( + (backwards_compatible and has_non_fixed_target or not backwards_compatible) and + time is not None + ): coords = [repeat_skycoord(coord, time) for coord in coords] # are all SkyCoordinate's in equivalent frames? If not, convert to ICRS diff --git a/astroplan/tests/test_constraints.py b/astroplan/tests/test_constraints.py index c0915c30..cb599d1f 100644 --- a/astroplan/tests/test_constraints.py +++ b/astroplan/tests/test_constraints.py @@ -462,7 +462,8 @@ def test_caches_shapes(): observer = Observer.at_site('lapalma') ac = AltitudeConstraint(min=30*u.deg) - # When time and targets have the same size, they are broadcastable so grid_times_targets is ignored + # When time and targets are the same size, + # they're broadcastable and grid_times_targets is ignored assert ac(observer, targets, times, grid_times_targets=True).shape == (3,) targets = get_skycoord([m31, ippeg]) assert ac(observer, targets, times, grid_times_targets=True).shape == (2, 3) diff --git a/astroplan/tests/test_scheduling.py b/astroplan/tests/test_scheduling.py index 25bd3b22..e40ca1ab 100644 --- a/astroplan/tests/test_scheduling.py +++ b/astroplan/tests/test_scheduling.py @@ -339,8 +339,8 @@ def test_scorer(): @pytest.mark.skipif('not HAS_SKYFIELD') def test_priority_scheduler_TLETarget(): - line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" - line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" + line1 = "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + line2 = "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" iss = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=apo) constraints = [AirmassConstraint(3, boolean_constraint=False)] blocks = [ObservingBlock(t, 5*u.minute, i) for i, t in enumerate(targets)] diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index ca77bdc9..046f5670 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -52,11 +52,12 @@ def test_FixedTarget_ra_dec(): assert vega.coord.dec == vega_coords.dec == vega.dec, ('Retrieve Dec from ' 'SkyCoord') + @pytest.mark.skipif('not HAS_SKYFIELD') def test_TLETarget(): - tle_string=("ISS (ZARYA)\n" - "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" - "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") + tle_string = ("ISS (ZARYA)\n" + "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" + "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" subaru = Observer.at_site('subaru') @@ -84,22 +85,23 @@ def test_TLETarget(): ra_dec1 = tle_target1.coord(times) ra_dec2 = tle_target2.coord(times) ra_dec_from_horizon = SkyCoord(["08h29m26.36s +07d31m24.0s", - "06h25m46.81s -54d32m20.2s", - "13h52m09.08s +04d26m40.1s", - "09h20m04.74s -00d51m19.2s"]) + "06h25m46.81s -54d32m20.2s", + "13h52m09.08s +04d26m40.1s", + "09h20m04.74s -00d51m19.2s"]) assert all(list(ra_dec1.separation(ra_dec_from_horizon) < 30*u.arcsec)) assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') # TLE Check - line1_invalid ="1 25544U 98067A .00041610 00000-0 73103-3 0 9990" + line1_invalid = "1 25544U 98067A .00041610 00000-0 73103-3 0 9990" with pytest.raises(ValueError): TLETarget(name="ISS (ZARYA)", line1=line1_invalid, line2=line2, observer=subaru) - TLETarget(name="ISS (ZARYA)", line1=line1_invalid, line2=line2, observer=subaru, skip_tle_check=True) + TLETarget(name="ISS (ZARYA)", line1=line1_invalid, line2=line2, observer=subaru, + skip_tle_check=True) # from_string - tle_string=("1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" - "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") + tle_string = ("1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" + "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") tle_target3 = TLETarget.from_string(tle_string=tle_string, observer=subaru, name="ISS") assert tle_target3.name == "ISS" @@ -107,7 +109,6 @@ def test_TLETarget(): with pytest.raises(ValueError): TLETarget.from_string(tle_string=tle_string, observer=subaru) - # AltAz (This is slow) altaz_observer = subaru.altaz(time, tle_target1) altaz_skyfield = tle_target1.altaz(time) @@ -115,7 +116,7 @@ def test_TLETarget(): assert altaz_observer.separation(altaz_skyfield) < 20*u.arcsec # Time too far in the future where elements stop making physical sense - with pytest.warns(): # ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6) + with pytest.warns(): # ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6) time_invalid = Time("2035-08-02 10:00", scale='utc') times_list = list(times) times_list[2] = Time("2035-08-02 10:00", scale='utc') @@ -167,8 +168,8 @@ def test_get_skycoord(): def test_get_skycoord_with_TLETarget(): skycoord_targed = SkyCoord(10.6847083*u.deg, 41.26875*u.deg) fixed_target1 = FixedTarget(name="fixed1", coord=SkyCoord(279.23458, 38.78369, unit='deg')) - line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" - line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" + line1 = "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + line2 = "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" subaru = Observer.at_site('subaru') tle_target1 = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=subaru) @@ -183,12 +184,12 @@ def test_get_skycoord_with_TLETarget(): tle_output = get_skycoord([tle_target1, tle_target1], time) assert tle_output.shape == (2,) tle_output = get_skycoord([tle_target1, tle_target1], times) - assert tle_output.shape == (2,4) + assert tle_output.shape == (2, 4) mixed_output = get_skycoord([skycoord_targed, fixed_target1, tle_target1], time) assert mixed_output.shape == (3,) mixed_output = get_skycoord([skycoord_targed, fixed_target1, tle_target1], times) - assert mixed_output.shape == (3,4) + assert mixed_output.shape == (3, 4) # backwards_compatible tle_output = get_skycoord(skycoord_targed, time, backwards_compatible=False) @@ -205,4 +206,4 @@ def test_get_skycoord_with_TLETarget(): tle_output = get_skycoord([skycoord_targed, skycoord_targed], times) assert tle_output.shape == (2,) tle_output = get_skycoord([skycoord_targed, skycoord_targed], times, backwards_compatible=False) - assert tle_output.shape == (2,4) + assert tle_output.shape == (2, 4) From 614c56fb3e4f0329d7ccbdcfccd1c5adb852953f Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Sat, 5 Aug 2023 20:12:54 +0200 Subject: [PATCH 03/12] Fix code style checks II --- astroplan/tests/test_target.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index 046f5670..a710b243 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -58,8 +58,8 @@ def test_TLETarget(): tle_string = ("ISS (ZARYA)\n" "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990\n" "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") - line1="1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" - line2="2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" + line1 = "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" + line2 = "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" subaru = Observer.at_site('subaru') time = Time("2023-08-02 10:00", scale='utc') times = time_grid_from_range([time, time + 3.1*u.hour], From 7ce8e6f0b43046e93a678fd494e5f944f59bac76 Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Thu, 10 Aug 2023 08:36:57 +0200 Subject: [PATCH 04/12] Tests: add more detailed tests to verify accuracy of TLETarget coordinates --- astroplan/tests/test_target.py | 54 +++++++++++++++++++++++++++------- 1 file changed, 44 insertions(+), 10 deletions(-) diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index a710b243..46666f94 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -60,7 +60,7 @@ def test_TLETarget(): "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092") line1 = "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" line2 = "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" - subaru = Observer.at_site('subaru') + subaru = Observer.at_site('subaru') # (lon, lat, el)=(-155.476111 deg, 19.825555 deg, 4139.0 m) time = Time("2023-08-02 10:00", scale='utc') times = time_grid_from_range([time, time + 3.1*u.hour], time_resolution=1 * u.hour) @@ -73,23 +73,57 @@ def test_TLETarget(): assert repr(tle_target1) == repr(tle_target2) assert str(tle_target1) == str(tle_target2) - # Single time - ra_dec1 = tle_target1.coord(time) + # Single time (Below Horizon) + ra_dec1 = tle_target1.coord(time) # '08h29m26.00003243s +07d31m36.65950907s' ra_dec2 = tle_target2.coord(time) - ra_dec_from_horizon = SkyCoord("08h29m26.36s +07d31m24.0s") - assert ra_dec1.separation(ra_dec_from_horizon) < 30*u.arcsec + # Comparison with the JPL Horizons System + ra_dec_horizon_icrf = SkyCoord("08h29m27.029117s +07d31m28.35610s") + # ICRF: Compensated for the down-leg light-time delay aberration + # ra_dec1.separation(ra_dec_horizon_icrf).to(u.arcsec) # Difference: 17.41″ + # ~ 2 * tan(17,41/2/3600) * 11801,56 = 57 km + + ra_dec_horizon_ref_apparent = SkyCoord("08h30m54.567398s +08d05m32.72764s") + # Refracted Apparent: In an equatorial coordinate system with all compensations + # ra_dec1.separation(ra_dec_horizon_ref_apparent).to(u.arcsec) # Difference: 2424.44″ + + ra_dec_horizon_icrf_ref_apparent = SkyCoord("08h29m37.373866s +08d10m14.78811s") + # ICRF Refracted Apparent: In the ICRF reference frame with all compensations + # ra_dec1.separation(ra_dec_horizon_icrf_ref_apparent).to(u.arcsec) # Difference: 2324.28″ + + # Skyfield appears to use no compensations. According to this, it's not even recommended to + # compensate for light travel time. Compensating changes the difference to ra_dec_horizon_icrf + # to 20.05 arcsec. + # https://rhodesmill.org/skyfield/earth-satellites.html#avoid-calling-the-observe-method + + assert ra_dec1.separation(ra_dec_horizon_icrf) < 20*u.arcsec assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') + # Single time (Above Horizon) + time_ah = Time("2023-08-02 07:20", scale='utc') + ra_dec_ah = tle_target1.coord(time_ah) # '11h19m48.53631001s +44d49m45.22194611s' + + ra_dec_ah_horizon_icrf = SkyCoord("11h19m49.660349s +44d49m34.65875s") # 15.95″ + ra_dec_ah_horizon_ref_apparent = SkyCoord("11h21m34.102381s +44d43m40.06899s") # 1181.84″ + ra_dec_ah_horizon_icrf_ref_apparent = SkyCoord("11h20m16.627261s +44d51m24.25337s") # 314.75″ + + # Default is WGS72 for Skyfield. Coordinates with WGS84 gravity model that Horizon uses: + # '11h19m48.28084569s +44d49m46.33649241s' - 18.75″ + # See 'Build a satellite with a specific gravity model' in Skyfield's Earth Satellites docu + + # There are many potential sources of inaccuracies, and it's not all super precise. + # Should the accuracy be better than < 25*u.arcsec when compared to the JPL Horizons System? + assert ra_dec_ah.separation(ra_dec_ah_horizon_icrf) < 20*u.arcsec + # Multiple times ra_dec1 = tle_target1.coord(times) ra_dec2 = tle_target2.coord(times) - ra_dec_from_horizon = SkyCoord(["08h29m26.36s +07d31m24.0s", - "06h25m46.81s -54d32m20.2s", - "13h52m09.08s +04d26m40.1s", - "09h20m04.74s -00d51m19.2s"]) + ra_dec_from_horizon = SkyCoord(["08h29m27.029117s +07d31m28.35610s", # 17.41″ + "06h25m46.672661s -54d32m16.77533s", # 22.05″ + "13h52m08.854291s +04d26m49.56432s", # 3.20″ + "09h20m04.872215s -00d51m21.17432s"]) # 17.55″ - assert all(list(ra_dec1.separation(ra_dec_from_horizon) < 30*u.arcsec)) + assert all(list(ra_dec1.separation(ra_dec_from_horizon) < 25*u.arcsec)) assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') # TLE Check From 8a0b473f0804415e565681d84b921bd97efe45f8 Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Thu, 10 Aug 2023 08:51:47 +0200 Subject: [PATCH 05/12] Fix code style checks III --- astroplan/tests/test_target.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index 46666f94..5bacbb33 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -76,36 +76,37 @@ def test_TLETarget(): # Single time (Below Horizon) ra_dec1 = tle_target1.coord(time) # '08h29m26.00003243s +07d31m36.65950907s' ra_dec2 = tle_target2.coord(time) + assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') # Comparison with the JPL Horizons System ra_dec_horizon_icrf = SkyCoord("08h29m27.029117s +07d31m28.35610s") # ICRF: Compensated for the down-leg light-time delay aberration - # ra_dec1.separation(ra_dec_horizon_icrf).to(u.arcsec) # Difference: 17.41″ - # ~ 2 * tan(17,41/2/3600) * 11801,56 = 57 km + assert ra_dec1.separation(ra_dec_horizon_icrf) < 20*u.arcsec # 17.41″ + # Distance estimation: ~ 2 * tan(17,41/2/3600) * 11801,56 = 57 km ra_dec_horizon_ref_apparent = SkyCoord("08h30m54.567398s +08d05m32.72764s") # Refracted Apparent: In an equatorial coordinate system with all compensations - # ra_dec1.separation(ra_dec_horizon_ref_apparent).to(u.arcsec) # Difference: 2424.44″ + assert ra_dec1.separation(ra_dec_horizon_ref_apparent) > 2000*u.arcsec # 2424.44″ ra_dec_horizon_icrf_ref_apparent = SkyCoord("08h29m37.373866s +08d10m14.78811s") # ICRF Refracted Apparent: In the ICRF reference frame with all compensations - # ra_dec1.separation(ra_dec_horizon_icrf_ref_apparent).to(u.arcsec) # Difference: 2324.28″ + assert ra_dec1.separation(ra_dec_horizon_icrf_ref_apparent) > 2000*u.arcsec # 2324.28″ # Skyfield appears to use no compensations. According to this, it's not even recommended to # compensate for light travel time. Compensating changes the difference to ra_dec_horizon_icrf # to 20.05 arcsec. # https://rhodesmill.org/skyfield/earth-satellites.html#avoid-calling-the-observe-method - assert ra_dec1.separation(ra_dec_horizon_icrf) < 20*u.arcsec - assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') - # Single time (Above Horizon) time_ah = Time("2023-08-02 07:20", scale='utc') ra_dec_ah = tle_target1.coord(time_ah) # '11h19m48.53631001s +44d49m45.22194611s' - ra_dec_ah_horizon_icrf = SkyCoord("11h19m49.660349s +44d49m34.65875s") # 15.95″ - ra_dec_ah_horizon_ref_apparent = SkyCoord("11h21m34.102381s +44d43m40.06899s") # 1181.84″ - ra_dec_ah_horizon_icrf_ref_apparent = SkyCoord("11h20m16.627261s +44d51m24.25337s") # 314.75″ + ra_dec_ah_horizon_icrf = SkyCoord("11h19m49.660349s +44d49m34.65875s") + assert ra_dec_ah.separation(ra_dec_ah_horizon_icrf) < 20*u.arcsec # 15.95″ + ra_dec_ah_horizon_ref_apparent = SkyCoord("11h21m34.102381s +44d43m40.06899s") + assert ra_dec_ah.separation(ra_dec_ah_horizon_ref_apparent) > 1000*u.arcsec # 1181.84″ + ra_dec_ah_horizon_icrf_ref_apparent = SkyCoord("11h20m16.627261s +44d51m24.25337s") + assert ra_dec_ah.separation(ra_dec_ah_horizon_icrf_ref_apparent) > 300*u.arcsec # 314.75″ # Default is WGS72 for Skyfield. Coordinates with WGS84 gravity model that Horizon uses: # '11h19m48.28084569s +44d49m46.33649241s' - 18.75″ @@ -113,16 +114,15 @@ def test_TLETarget(): # There are many potential sources of inaccuracies, and it's not all super precise. # Should the accuracy be better than < 25*u.arcsec when compared to the JPL Horizons System? - assert ra_dec_ah.separation(ra_dec_ah_horizon_icrf) < 20*u.arcsec # Multiple times ra_dec1 = tle_target1.coord(times) ra_dec2 = tle_target2.coord(times) - ra_dec_from_horizon = SkyCoord(["08h29m27.029117s +07d31m28.35610s", # 17.41″ - "06h25m46.672661s -54d32m16.77533s", # 22.05″ - "13h52m08.854291s +04d26m49.56432s", # 3.20″ - "09h20m04.872215s -00d51m21.17432s"]) # 17.55″ - + ra_dec_from_horizon = SkyCoord(["08h29m27.029117s +07d31m28.35610s", + "06h25m46.672661s -54d32m16.77533s", + "13h52m08.854291s +04d26m49.56432s", + "09h20m04.872215s -00d51m21.17432s"]) + # 17.41″, 22.05″, 3.20″, 17.55″ assert all(list(ra_dec1.separation(ra_dec_from_horizon) < 25*u.arcsec)) assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') From 8e1f83a8b6057216981ce9390cbd5ef09bc70a81 Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Wed, 26 Feb 2025 13:52:38 +0100 Subject: [PATCH 06/12] Tests: warnings for invalid TLE changed; implement fix --- astroplan/tests/test_scheduling.py | 5 ++++- astroplan/tests/test_target.py | 10 +++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/astroplan/tests/test_scheduling.py b/astroplan/tests/test_scheduling.py index e40ca1ab..c2c9cafa 100644 --- a/astroplan/tests/test_scheduling.py +++ b/astroplan/tests/test_scheduling.py @@ -369,7 +369,10 @@ def test_priority_scheduler_TLETarget(): start_time = Time('2035-08-02 10:00:00') end_time = start_time + 1*u.hour schedule = Schedule(start_time, end_time) - with pytest.warns(InvalidTLEDataWarning): + + # InvalidTLEDataWarning/AstropyWarning and + # ErfaWarning: ERFA function "utctai" yielded 121 of "dubious year (Note 3)" + with pytest.warns(): scheduler(blocks, schedule) assert len(schedule.observing_blocks) == 3 assert all([schedule.observing_blocks[0].target == vega, diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index 5bacbb33..170d5993 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -150,14 +150,18 @@ def test_TLETarget(): assert altaz_observer.separation(altaz_skyfield) < 20*u.arcsec # Time too far in the future where elements stop making physical sense - with pytest.warns(): # ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6) + with pytest.warns(): # ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6) time_invalid = Time("2035-08-02 10:00", scale='utc') times_list = list(times) times_list[2] = Time("2035-08-02 10:00", scale='utc') times_invalid = Time(times_list) - with pytest.warns(InvalidTLEDataWarning): + # InvalidTLEDataWarning and + # ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" + with pytest.warns(): assert np.isnan(tle_target1.coord(time_invalid).ra) - with pytest.warns(InvalidTLEDataWarning): + # InvalidTLEDataWarning and + # ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" + with pytest.warns(): assert np.isnan(tle_target1.coord(times_invalid)[2].ra) From 4c43d3da74b7c4253758aeb3bc1947fa73591667 Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Wed, 26 Feb 2025 14:45:12 +0100 Subject: [PATCH 07/12] TLETarget: support only astroplan Observer and use its properties for atmospheric refraction in `ICRF.altaz()` - Update tests with an example without an observer and one with temperature and pressure. Also, use time above the horizon to test `TLETarget.altaz()`, and test `TLETarget.altaz()` with multiple times. - Use `@pytest.mark.remote_data` for TLE tests since they use `Observer.at_site()`. --- astroplan/target.py | 56 +++++++++++++++++++++++----------- astroplan/tests/test_target.py | 28 +++++++++++++++-- 2 files changed, 64 insertions(+), 20 deletions(-) diff --git a/astroplan/target.py b/astroplan/target.py index 8a020ab1..61407ac3 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -220,7 +220,7 @@ def __init__(self, line1, line2, name=None, observer=None, skip_tle_check=False) Name of the target, used for plotting and representing the target as a string - observer : `~astropy.coordinates.EarthLocation`, `~astroplan.Observer`, optional + observer : `~astroplan.Observer`, optional The location of observer. If `None`, the observer is assumed to be at sea level at the equator. @@ -237,13 +237,16 @@ def __init__(self, line1, line2, name=None, observer=None, skip_tle_check=False) self.satellite = EarthSatellite(line1, line2, name, load.timescale()) if observer is None: - self.observer = wgs84.latlon(0, 0, 0) + # Prevent circular import and usually not used + from .observer import Observer + self.observer = Observer(latitude=0*u.deg, longitude=0*u.deg, elevation=0*u.m) else: - observer = getattr(observer, 'location', observer) - longitude, latitude, height = observer.to_geodetic() - self.observer = wgs84.latlon(latitude.to(u.deg).value, - longitude.to(u.deg).value, - height.to(u.m).value) + self.observer = observer + + longitude, latitude, height = self.observer.location.to_geodetic() + self.geographic_position = wgs84.latlon(latitude.to(u.deg).value, + longitude.to(u.deg).value, + height.to(u.m).value) @classmethod def from_string(cls, tle_string, name=None, *args, **kwargs): @@ -300,7 +303,7 @@ def _compute_topocentric(self, time=None): ts = load.timescale() t = ts.from_astropy(time) - topocentric = (self.satellite - self.observer).at(t) + topocentric = (self.satellite - self.geographic_position).at(t) # Check for invalid TLE data. A non-None usually message means the computation went beyond # the physically sensible point. Details: @@ -333,7 +336,13 @@ def coord(self, time=None): topocentric = self._compute_topocentric(time) ra, dec, distance = topocentric.radec() # No distance, in SkyCoord, distance is from frame origin, but here, it's from observer. - return SkyCoord(ra.hours*u.hourangle, dec.degrees*u.deg, obstime=time, frame='icrs') + return SkyCoord( + ra.hours*u.hourangle, + dec.degrees*u.deg, + obstime=time, + frame='icrs', + location=self.observer.location, + ) def altaz(self, time=None): """ @@ -350,15 +359,28 @@ def altaz(self, time=None): SkyCoord object representing the target's altitude and azimuth at the specified time(s) """ topocentric = self._compute_topocentric(time) - alt, az, distance = topocentric.altaz() - - earth_location = EarthLocation(lat=self.observer.latitude.degrees*u.deg, - lon=self.observer.longitude.degrees*u.deg, - height=self.observer.elevation.m*u.m) - altaz = AltAz(alt=alt.degrees*u.deg, az=az.degrees*u.deg, - obstime=time, location=earth_location) - return SkyCoord(altaz) + temperature_C = None + pressure_mbar = None + if self.observer.temperature is not None: + temperature_C = self.observer.temperature.to_value(u.deg_C) + if self.observer.pressure is not None: + pressure_mbar = self.observer.pressure.to_value(u.mbar) + + alt, az, distance = topocentric.altaz( + temperature_C=temperature_C, + pressure_mbar=pressure_mbar, + ) + + # 'relative_humidity' and 'obswl' were not used in coordinate calculation + altaz_frame = AltAz( + location=self.observer.location, + obstime=time, + pressure=self.observer.pressure, + temperature=self.observer.temperature, + ) + + return SkyCoord(alt=alt.degrees*u.deg, az=az.degrees*u.deg, frame=altaz_frame) def __repr__(self): return f'<{self.__class__.__name__} "{self.name}">' diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index 170d5993..50626ac7 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -53,6 +53,7 @@ def test_FixedTarget_ra_dec(): 'SkyCoord') +@pytest.mark.remote_data @pytest.mark.skipif('not HAS_SKYFIELD') def test_TLETarget(): tle_string = ("ISS (ZARYA)\n" @@ -61,18 +62,27 @@ def test_TLETarget(): line1 = "1 25544U 98067A 23215.27256123 .00041610 00000-0 73103-3 0 9990" line2 = "2 25544 51.6403 95.2411 0000623 157.9606 345.0624 15.50085581409092" subaru = Observer.at_site('subaru') # (lon, lat, el)=(-155.476111 deg, 19.825555 deg, 4139.0 m) + subaru_temp_pressure = Observer.at_site('subaru', temperature=-10*u.Celsius, pressure=0.9*u.bar) time = Time("2023-08-02 10:00", scale='utc') times = time_grid_from_range([time, time + 3.1*u.hour], time_resolution=1 * u.hour) tle_target1 = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=subaru) tle_target2 = TLETarget.from_string(tle_string=tle_string, observer=subaru) + tle_target_no_observer = TLETarget(name="ISS (ZARYA)", line1=line1, line2=line2, observer=None) + tle_target_temp_pressure = TLETarget( + name="ISS (ZARYA)", line1=line1, line2=line2, observer=subaru_temp_pressure + ) assert tle_target1.name == "ISS (ZARYA)" assert tle_target2.name == "ISS (ZARYA)" assert repr(tle_target1) == repr(tle_target2) assert str(tle_target1) == str(tle_target2) + assert isinstance(tle_target_no_observer.observer, Observer) + assert abs(tle_target_no_observer.observer.location.lat) < 0.001*u.deg + tle_target_no_observer.coord(time) # Just needs to work + # Single time (Below Horizon) ra_dec1 = tle_target1.coord(time) # '08h29m26.00003243s +07d31m36.65950907s' ra_dec2 = tle_target2.coord(time) @@ -144,10 +154,21 @@ def test_TLETarget(): TLETarget.from_string(tle_string=tle_string, observer=subaru) # AltAz (This is slow) - altaz_observer = subaru.altaz(time, tle_target1) - altaz_skyfield = tle_target1.altaz(time) + altaz_observer = subaru.altaz(time_ah, tle_target1) + altaz_skyfield = tle_target1.altaz(time_ah) + + assert altaz_observer.separation(altaz_skyfield) < 21*u.arcsec - assert altaz_observer.separation(altaz_skyfield) < 20*u.arcsec + # AltAz with atmospheric refraction (temperature and pressure) + altaz_observer = subaru_temp_pressure.altaz(time_ah, tle_target_temp_pressure) + altaz_skyfield = tle_target_temp_pressure.altaz(time_ah) + + assert altaz_observer.separation(altaz_skyfield) < 26*u.arcsec + + # AltAz with multiple times + #subaru.altaz(times, tle_target1) + altaz_multiple = tle_target1.altaz(times) + assert len(altaz_multiple.obstime) == len(altaz_multiple) == len(times) # Time too far in the future where elements stop making physical sense with pytest.warns(): # ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6) @@ -202,6 +223,7 @@ def test_get_skycoord(): assert len(coo) == 2 +@pytest.mark.remote_data @pytest.mark.skipif('not HAS_SKYFIELD') def test_get_skycoord_with_TLETarget(): skycoord_targed = SkyCoord(10.6847083*u.deg, 41.26875*u.deg) From 4a0e93f8d56806a880c91b6fa620c31c1717d5ba Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Mon, 2 Feb 2026 09:57:59 +0100 Subject: [PATCH 08/12] Fix code style checks IIII --- astroplan/target.py | 6 +++--- astroplan/tests/test_scheduling.py | 3 +-- astroplan/tests/test_target.py | 3 +-- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/astroplan/target.py b/astroplan/target.py index 61407ac3..755d666d 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -7,7 +7,7 @@ # Third-party import astropy.units as u from astropy.time import Time -from astropy.coordinates import SkyCoord, ICRS, UnitSphericalRepresentation, AltAz, EarthLocation +from astropy.coordinates import SkyCoord, ICRS, UnitSphericalRepresentation, AltAz try: from sgp4.io import twoline2rv from sgp4.earth_gravity import wgs84 as sgp4_wgs84 @@ -366,12 +366,12 @@ def altaz(self, time=None): temperature_C = self.observer.temperature.to_value(u.deg_C) if self.observer.pressure is not None: pressure_mbar = self.observer.pressure.to_value(u.mbar) - + alt, az, distance = topocentric.altaz( temperature_C=temperature_C, pressure_mbar=pressure_mbar, ) - + # 'relative_humidity' and 'obswl' were not used in coordinate calculation altaz_frame = AltAz( location=self.observer.location, diff --git a/astroplan/tests/test_scheduling.py b/astroplan/tests/test_scheduling.py index c2c9cafa..8a65a0e1 100644 --- a/astroplan/tests/test_scheduling.py +++ b/astroplan/tests/test_scheduling.py @@ -19,7 +19,6 @@ from astroplan.periodic import EclipsingSystem from astroplan.scheduling import (ObservingBlock, PriorityScheduler, SequentialScheduler, Transitioner, TransitionBlock, Schedule, Slot, Scorer) -from astroplan.exceptions import InvalidTLEDataWarning vega = FixedTarget(coord=SkyCoord(ra=279.23473479 * u.deg, dec=38.78368896 * u.deg), name="Vega") @@ -369,7 +368,7 @@ def test_priority_scheduler_TLETarget(): start_time = Time('2035-08-02 10:00:00') end_time = start_time + 1*u.hour schedule = Schedule(start_time, end_time) - + # InvalidTLEDataWarning/AstropyWarning and # ErfaWarning: ERFA function "utctai" yielded 121 of "dubious year (Note 3)" with pytest.warns(): diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index 50626ac7..13d07145 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -16,7 +16,6 @@ from astroplan.target import FixedTarget, TLETarget, get_skycoord from astroplan.observer import Observer from astroplan.utils import time_grid_from_range -from astroplan.exceptions import InvalidTLEDataWarning @pytest.mark.remote_data @@ -166,7 +165,7 @@ def test_TLETarget(): assert altaz_observer.separation(altaz_skyfield) < 26*u.arcsec # AltAz with multiple times - #subaru.altaz(times, tle_target1) + # subaru.altaz(times, tle_target1) altaz_multiple = tle_target1.altaz(times) assert len(altaz_multiple.obstime) == len(altaz_multiple) == len(times) From ca30cb3b9fe667e94d93abcc8ed34e71d9b2cdfd Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Thu, 12 Feb 2026 13:26:01 +0100 Subject: [PATCH 09/12] TLETarget / time-dependent targets: make compatible with AltAzTarget --- astroplan/constraints.py | 42 +++-- astroplan/observer.py | 37 ++++- astroplan/scheduling.py | 96 +++++++---- astroplan/target.py | 246 +++++++++++++--------------- astroplan/tests/test_constraints.py | 12 +- astroplan/tests/test_target.py | 33 +--- 6 files changed, 237 insertions(+), 229 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 3be2fdeb..06d841ca 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -223,49 +223,45 @@ def __call__(self, observer, targets, times=None, time_range=None, time_grid_resolution=0.5*u.hour, grid_times_targets=False): """ - Compute the constraint for this class + Compute the constraint for this class. Parameters ---------- observer : `~astroplan.Observer` - the observation location from which to apply the constraints - targets : sequence of `~astroplan.Target` + The observation location from which to apply the constraints. + + targets : sequence of `~astroplan.Target` or `~astropy.coordinates.SkyCoord` The targets on which to apply the constraints. - times : `~astropy.time.Time` + times : `~astropy.time.Time` (optional) The times to compute the constraint. - WHAT HAPPENS WHEN BOTH TIMES AND TIME_RANGE ARE SET? - time_range : `~astropy.time.Time` (length = 2) + time_range : `~astropy.time.Time` (length = 2) (optional) Lower and upper bounds on time sequence. + Only used when `times` is not provided. time_grid_resolution : `~astropy.units.Quantity` - Time-grid spacing + Time-grid spacing. grid_times_targets : bool - if True, grids the constraint result with targets along the first + If True, grids the constraint result with targets along the first index and times along the second. Otherwise, we rely on broadcasting the shapes together using standard numpy rules. + Returns ------- constraint_result : 1D or 2D array of float or bool - The constraints. If 2D with targets along the first index and times along + The constraint values. If 2D with targets along the first index and times along the second. """ - if times is None and time_range is not None: times = time_grid_from_range(time_range, time_resolution=time_grid_resolution) - if grid_times_targets: - targets = get_skycoord(targets, times) - # TODO: these broadcasting operations are relatively slow - # but there is potential for huge speedup if the end user - # disables gridding and re-shapes the coords themselves - # prior to evaluating multiple constraints. - if not observer._is_broadcastable(targets.shape, times.shape): - if targets.isscalar: - # ensure we have a (1, 1) shape coord - targets = SkyCoord(np.tile(targets, 1))[:, np.newaxis] - else: - targets = targets[..., np.newaxis] - times, targets = observer._preprocess_inputs(times, targets, grid_times_targets=False) + # TODO: broadcasting operations in _preprocess_inputs are relatively slow + # but there is potential for huge speedup if the end user + # disables gridding and re-shapes the coords themselves + # prior to evaluating multiple constraints. + times, targets = observer._preprocess_inputs( + times, targets, grid_times_targets=grid_times_targets + ) + result = self.compute_constraint(times, observer, targets) # make sure the output has the same shape as would result from diff --git a/astroplan/observer.py b/astroplan/observer.py index def4d2cc..6cc4e9cd 100644 --- a/astroplan/observer.py +++ b/astroplan/observer.py @@ -507,22 +507,43 @@ def _preprocess_inputs(self, time, target=None, grid_times_targets=False): the shapes together using standard numpy rules. Useful for grid searches for rise/set times etc. """ - # make sure we have a non-scalar time if not isinstance(time, Time): time = Time(time) + # In grid mode, scalar time should still behave like a length-1 time axis + if grid_times_targets and time.isscalar: + time = time[None] # shape (1,) + if target is None: return time, None + # Remember whether target is a single time-dependent target + is_multiple_targets = ( + isinstance(target, (list, tuple)) or + (isinstance(target, SkyCoord) and not target.isscalar) + ) + is_target_time_dependent = ( + callable(getattr(target, "get_skycoord", None)) and not + hasattr(target, "coord") + ) + is_single_time_dependent_target = (not is_multiple_targets) and is_target_time_dependent + # convert any kind of target argument to non-scalar SkyCoord - target = get_skycoord(target, time) - if grid_times_targets and not self._is_broadcastable(target.shape, time.shape): + target = get_skycoord(target, times=time) + + if grid_times_targets: + # Only ambiguous case: a single time-dependent target produces shape == time.shape + # but grid mode requires a leading target axis (1, ...). + if is_single_time_dependent_target and (not target.isscalar) and (target.shape == time.shape): + target = target[np.newaxis, ...] + + # Ensure at least one targets axis for scalar targets if target.isscalar: - # ensure we have a (1, 1) shape coord - target = SkyCoord(np.tile(target, 1))[:, np.newaxis] - else: - while target.ndim <= time.ndim: - target = target[:, np.newaxis] + target = SkyCoord(np.tile(target, 1)) # shape (1,) + + # Make target have one more dim than time (first targets axis, then time axes) + while target.ndim < 1 + time.ndim: + target = target[..., np.newaxis] elif not self._is_broadcastable(target.shape, time.shape): raise ValueError('Time and Target arguments cannot be broadcast ' diff --git a/astroplan/scheduling.py b/astroplan/scheduling.py index 867fdb88..24220d68 100644 --- a/astroplan/scheduling.py +++ b/astroplan/scheduling.py @@ -14,7 +14,7 @@ from .utils import time_grid_from_range, stride_array from .constraints import AltitudeConstraint -from .target import get_skycoord, FixedTarget, TLETarget +from .target import get_skycoord __all__ = ['ObservingBlock', 'TransitionBlock', 'Schedule', 'Slot', 'Scheduler', 'SequentialScheduler', 'PriorityScheduler', @@ -118,6 +118,7 @@ def __init__(self, blocks, observer, schedule, global_constraints=[]): self.observer = observer self.schedule = schedule self.global_constraints = global_constraints + self.targets = [block.target for block in self.blocks] def create_score_array(self, time_resolution=1*u.minute): """ @@ -138,7 +139,6 @@ def create_score_array(self, time_resolution=1*u.minute): start = self.schedule.start_time end = self.schedule.end_time times = time_grid_from_range((start, end), time_resolution) - targets = get_skycoord([block.target for block in self.blocks], times) score_array = np.ones((len(self.blocks), len(times))) for i, block in enumerate(self.blocks): # TODO: change the default constraints from None to [] @@ -147,6 +147,7 @@ def create_score_array(self, time_resolution=1*u.minute): applied_score = constraint(self.observer, block.target, times=times) score_array[i] *= applied_score + targets = get_skycoord(self.targets, times=times) for constraint in self.global_constraints: score_array *= constraint(self.observer, targets, times, grid_times_targets=True) @@ -265,54 +266,81 @@ def open_slots(self): return [slot for slot in self.slots if not slot.occupied] def to_table(self, show_transitions=True, show_unused=False): - # TODO: allow different coordinate types + def _format_target_info(target): + if hasattr(target, "coord"): + try: + return target.coord.icrs.to_string("hmsdms") + except Exception: + return repr(target.coord) + if hasattr(target, "alt") and hasattr(target, "az"): + parts = [ + "alt={:.6f} deg".format(u.Quantity(target.alt).to_value(u.deg)), + "az={:.6f} deg".format(u.Quantity(target.az).to_value(u.deg)), + ] + if hasattr(target, "pressure"): + try: + p = u.Quantity(target.pressure) + if p.to_value(u.hPa) != 0.0: + parts.append("pressure={:.3f} hPa".format(p.to_value(u.hPa))) + except Exception: + pass + return ", ".join(parts) + if hasattr(target, "satellite"): + return ( + f"#{target.satellite.model.satnum} " + f"epoch {target.satellite.epoch.utc_strftime(format='%Y-%m-%d %H:%M:%S')}" + ) + return "" + target_names = [] start_times = [] end_times = [] durations = [] - coordiante_type = [] - coordinate_info = [] + target_types = [] + target_info = [] config = [] + for slot in self.slots: - if hasattr(slot.block, 'target'): + if hasattr(slot.block, "target"): start_times.append(slot.start.iso) end_times.append(slot.end.iso) - durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) - target = slot.block.target - target_names.append(target.name) - if isinstance(target, FixedTarget): - coordiante_type.append("RA/Dec") - coordinate_info.append(target.coord.to_string('hmsdms')) - elif isinstance(target, TLETarget): - coordiante_type.append("TLE") - coordinate_info.append( - f"#{target.satellite.model.satnum} " - f"epoch {target.satellite.epoch.utc_strftime(format='%Y-%m-%d %H:%M:%S')}" - ) + durations.append(slot.duration.to(u.minute).value) + target_names.append(slot.block.target.name) + target_types.append(slot.block.target.__class__.__name__) + target_info.append(_format_target_info(slot.block.target)) config.append(slot.block.configuration) elif show_transitions and slot.block: start_times.append(slot.start.iso) end_times.append(slot.end.iso) - durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) - target_names.append('TransitionBlock') - coordiante_type.append('') - coordinate_info.append('') + durations.append(slot.duration.to(u.minute).value) + target_names.append("TransitionBlock") + target_types.append("TransitionBlock") + target_info.append("") changes = list(slot.block.components.keys()) - if 'slew_time' in changes: - changes.remove('slew_time') + if "slew_time" in changes: + changes.remove("slew_time") config.append(changes) elif slot.block is None and show_unused: start_times.append(slot.start.iso) end_times.append(slot.end.iso) - durations.append('{:.4f}'.format(slot.duration.to(u.minute).value)) - target_names.append('Unused Time') - coordiante_type.append('') - coordinate_info.append('') - config.append('') - return Table([target_names, start_times, end_times, - durations, coordiante_type, coordinate_info, config], - names=('target', 'start time (UTC)', 'end time (UTC)', - 'duration (min)', 'type', 'coordinates/info', 'configuration')) + durations.append(slot.duration.to(u.minute).value) + target_names.append("Unused Time") + target_types.append("") + target_info.append("") + config.append("") + + return Table( + [target_names, start_times, end_times, durations, target_types, target_info, config], + names=( + "target", + "start time (UTC)", + "end time (UTC)", + "duration (minutes)", + "target type", + "target info", + "configuration", + ), + ) def new_slots(self, slot_index, start_time, end_time): """ @@ -1006,7 +1034,7 @@ def __call__(self, oldblock, newblock, start_time, observer): # to observer from .constraints import _get_altaz if oldblock.target != newblock.target: - targets = get_skycoord([oldblock.target, newblock.target], start_time) + targets = get_skycoord([oldblock.target, newblock.target], times=start_time) aaz = _get_altaz(start_time, observer, targets)['altaz'] sep = aaz[0].separation(aaz[1]) if sep/self.slew_rate > 1 * u.second: diff --git a/astroplan/target.py b/astroplan/target.py index 755d666d..0afc590c 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -5,6 +5,7 @@ import warnings # Third-party +import numpy as np import astropy.units as u from astropy.time import Time from astropy.coordinates import SkyCoord, ICRS, UnitSphericalRepresentation, AltAz @@ -21,7 +22,7 @@ from .exceptions import InvalidTLEDataWarning -__all__ = ["Target", "FixedTarget", "NonFixedTarget", "TLETarget"] +__all__ = ["Target", "FixedTarget", "TLETarget", "NonFixedTarget"] # Docstring code examples include printed SkyCoords, but the format changed # in astropy 1.3. Thus the doctest needs astropy >=1.3 and this is the @@ -216,15 +217,15 @@ def __init__(self, line1, line2, name=None, observer=None, skip_tle_check=False) line2 : str The second line of the TLE set - name : str, optional + name : str (optional) Name of the target, used for plotting and representing the target as a string - observer : `~astroplan.Observer`, optional + observer : `~astroplan.Observer` (optional) The location of observer. If `None`, the observer is assumed to be at sea level at the equator. - skip_tle_check : bool, optional + skip_tle_check : bool (optional) Whether to skip TLE validation """ if not skyfield_available: @@ -258,11 +259,11 @@ def from_string(cls, tle_string, name=None, *args, **kwargs): tle_string : str String to be parsed, expected to contain 2 or 3 newline-separated lines. - name : str, optional + name : str (optional) Name of the target. If not provided and the tle_string contains 3 lines, the first line will be used as the name. - args, kwargs : tuple, dict, optional + args, kwargs : tuple, dict (optional) Additional arguments and keyword arguments to be passed to the TLETarget constructor. """ lines = tle_string.strip().splitlines() @@ -276,21 +277,13 @@ def from_string(cls, tle_string, name=None, *args, **kwargs): line1, line2 = lines return cls(line1, line2, name, *args, **kwargs) - @property - def ra(self): - raise NotImplementedError("Compute satellite RA at a specific time with self.coord(time)") - - @property - def dec(self): - raise NotImplementedError("Compute satellite RA at a specific time with self.coord(time)") - - def _compute_topocentric(self, time=None): + def _compute_topocentric(self, times=None): """ Compute the topocentric coordinates (relative to observer) at a particular time. Parameters ---------- - time : `~astropy.time.Time`, optional + times : `~astropy.time.Time` (optional) The time(s) to use in the calculation. Returns @@ -298,10 +291,10 @@ def _compute_topocentric(self, time=None): topocentric : `skyfield.positionlib.ICRF` The topocentric object representing the relative coordinates of the target. """ - if time is None: - time = Time.now() + if times is None: + times = Time.now() ts = load.timescale() - t = ts.from_astropy(time) + t = ts.from_astropy(times) topocentric = (self.satellite - self.geographic_position).at(t) @@ -317,13 +310,13 @@ def _compute_topocentric(self, time=None): return topocentric - def coord(self, time=None): + def get_skycoord(self, times=None): """ Get the coordinates of the target at a particular time. Parameters ---------- - time : `~astropy.time.Time`, optional + times : `~astropy.time.Time` (optional) The time(s) to use in the calculation. Returns @@ -333,24 +326,24 @@ def coord(self, time=None): RA/Dec coordinates at the specified time(s). Might return np.nan and output a warning for times where the elements stop making physical sense. """ - topocentric = self._compute_topocentric(time) + topocentric = self._compute_topocentric(times) ra, dec, distance = topocentric.radec() # No distance, in SkyCoord, distance is from frame origin, but here, it's from observer. return SkyCoord( ra.hours*u.hourangle, dec.degrees*u.deg, - obstime=time, + obstime=times, frame='icrs', location=self.observer.location, ) - def altaz(self, time=None): + def altaz(self, times=None): """ Get the altitude and azimuth of the target at a particular time. Parameters ---------- - time : `~astropy.time.Time`, optional + times : `~astropy.time.Time` (optional) The time(s) to use in the calculation. Returns @@ -358,7 +351,7 @@ def altaz(self, time=None): altaz_coord : `~astropy.coordinates.SkyCoord` SkyCoord object representing the target's altitude and azimuth at the specified time(s) """ - topocentric = self._compute_topocentric(time) + topocentric = self._compute_topocentric(times) temperature_C = None pressure_mbar = None @@ -375,7 +368,7 @@ def altaz(self, time=None): # 'relative_humidity' and 'obswl' were not used in coordinate calculation altaz_frame = AltAz( location=self.observer.location, - obstime=time, + obstime=times, pressure=self.observer.pressure, temperature=self.observer.temperature, ) @@ -389,38 +382,7 @@ def __str__(self): return self.name -def repeat_skycoord(coord, times): - """ - Repeats the coordinates of a SkyCoord object 'times.size' number of times. - - Parameters - ---------- - coord : `~astropy.coordinates.SkyCoord` - The original SkyCoord object whose coordinates need to be repeated. - - times : `~astropy.time.Time` - The size of times determines the number of times the coordinates should be repeated. - - Returns - -------- - SkyCoord : `~astropy.coordinates.SkyCoord` - A new SkyCoord object with the coordinates of the original object - repeated 'times.size' number of times. If the SkyCoord object is scalar - or 'times' is None or a scalar, this function returns the - original SkyCoord object. - """ - if coord.size != 1 or times is None or times.size == 1: - return coord - return SkyCoord( - ra=coord.ra.repeat(times.size), - dec=coord.dec.repeat(times.size), - distance=None if coord.distance.unit is u.one else coord.distance.repeat(times.size), - frame=coord.frame, - obstime=times - ) - - -def get_skycoord(targets, time=None, backwards_compatible=True): +def get_skycoord(targets, times=None): """ Return an `~astropy.coordinates.SkyCoord` object. @@ -428,67 +390,61 @@ def get_skycoord(targets, time=None, backwards_compatible=True): a single `~astropy.coordinates.SkyCoord` object, rather than a list of `Target` or `~astropy.coordinates.SkyCoord` objects. - This is a convenience routine to do that. + This is a convenience routine to do that, and it also supports targets + that require evaluation at specific times (e.g., AltAz-defined targets). Parameters - ----------- - targets : list, `~astropy.coordinates.SkyCoord`, `Target` - either a single target or a list of targets - - time : `~astropy.time.Time`, optional - The time(s) to use in the calculation. + ---------- + targets : list, `~astropy.coordinates.SkyCoord`, `~astroplan.Target` + Either a single target or a list of targets. - backwards_compatible : bool, optional - Controls output format when FixedTarget or SkyCoord targets are used with a time argument. - If False, it will return (targets x times), where all coordinates per target are the same. - If True, it will return one coordinate per target (default is True). + times : `~astropy.time.Time` or time-like (optional) + Times at which to evaluate time-dependent targets. Required if any + target in ``targets`` needs evaluation at a time. Returns - -------- + ------- coord : `~astropy.coordinates.SkyCoord` - a single SkyCoord object, which may be non-scalar + A single SkyCoord object, which may be non-scalar. If ``times`` + is provided and any target is time-dependent, coordinates are broadcast + or evaluated across time along subsequent axes. """ - - # Note on backwards_compatible: - # Method always returns (targets x times) with TLETarget in targets, as RA/Dec changes with time - # Do we want to be 100% backwards compatible, or do we prefer consistent output, - # for FixedTarget or SkyCoord targets combined with multiple times? - # backwards_compatible = True will continue to return one coordinate per target - # backwards_compatible = False, returns (targets x times) with identical coordinates per target - - # Early exit for single target - if not isinstance(targets, list): - if isinstance(targets, TLETarget): - return targets.coord(time) - else: - if backwards_compatible: - return getattr(targets, 'coord', targets) - else: - return repeat_skycoord(getattr(targets, 'coord', targets), time) - - # Identify if any of the targets is not FixedTarget or SkyCoord - has_non_fixed_target = any( - not isinstance(target, (FixedTarget, SkyCoord)) - for target in targets + if times is not None and not isinstance(times, Time): + times = Time(times) + + def _is_time_dependent(obj): + return callable(getattr(obj, "get_skycoord", None)) and not hasattr(obj, "coord") + + def _as_coord(obj): + if hasattr(obj, "coord"): + return obj.coord + if callable(getattr(obj, "get_skycoord", None)): + return obj.get_skycoord(times) + return obj + + is_multiple_targets = ( + isinstance(targets, (list, tuple)) or + (isinstance(targets, SkyCoord) and not targets.isscalar) ) + if not is_multiple_targets: + return _as_coord(targets) + + coords = [_as_coord(t) for t in targets] + + # If any target is time dependent, broadcast fixed coords to match times.shape + time_dependent = (times is not None) and any(_is_time_dependent(t) for t in targets) + times_shape = times.shape if time_dependent else None - # Get the SkyCoord object itself - coords = [ - target.coord(time) if isinstance(target, TLETarget) - else getattr(target, 'coord', target) - for target in targets - ] - - # Fill time dimension for SkyCoords that only have a single coordinate - if ( - (backwards_compatible and has_non_fixed_target or not backwards_compatible) and - time is not None - ): - coords = [repeat_skycoord(coord, time) for coord in coords] - - # are all SkyCoordinate's in equivalent frames? If not, convert to ICRS + def _broadcast_quantity(q, shape): + """Broadcast quantity to target shape if needed.""" + if shape is None or q.shape == shape: + return q + return u.Quantity(np.broadcast_to(q.to_value(q.unit), shape), q.unit) + + # Are all SkyCoord's in equivalent frames? If not, convert to ICRS convert_to_icrs = not all( - [coord.frame.is_equivalent_frame(coords[0].frame) for coord in coords[1:]]) + [coord.frame.is_equivalent_frame(coords[0].frame) for coord in coords[1:]] + ) # we also need to be careful about handling mixtures of # UnitSphericalRepresentations and others @@ -503,10 +459,18 @@ def get_skycoord(targets, time=None, backwards_compatible=True): # mixture of frames for coordinate in coords: icrs_coordinate = coordinate.icrs - longitudes.append(icrs_coordinate.ra) - latitudes.append(icrs_coordinate.dec) + lon = icrs_coordinate.ra + lat = icrs_coordinate.dec + if times_shape is not None: + lon = _broadcast_quantity(lon, times_shape) + lat = _broadcast_quantity(lat, times_shape) + longitudes.append(lon) + latitudes.append(lat) if get_distances: - distances.append(icrs_coordinate.distance) + dist = icrs_coordinate.distance + if times_shape is not None: + dist = _broadcast_quantity(dist, times_shape) + distances.append(dist) frame = ICRS() else: # all the same frame, get the longitude and latitude names @@ -521,27 +485,53 @@ def get_skycoord(targets, time=None, backwards_compatible=True): frame = coords[0].frame for coordinate in coords: - longitudes.append(getattr(coordinate, lon_name)) - latitudes.append(getattr(coordinate, lat_name)) + lon = getattr(coordinate, lon_name) + lat = getattr(coordinate, lat_name) + if times_shape is not None: + lon = _broadcast_quantity(lon, times_shape) + lat = _broadcast_quantity(lat, times_shape) + longitudes.append(lon) + latitudes.append(lat) if get_distances: - distances.append(coordinate.distance) + dist = coordinate.distance + if times_shape is not None: + dist = _broadcast_quantity(dist, times_shape) + distances.append(dist) + + # Convert all longitude/latitude quantities to a common unit + # and plain ndarrays before stacking (robust across units/Quantity subclasses). + lon_unit = longitudes[0].unit + lat_unit = latitudes[0].unit + lon_vals = np.stack([lon.to_value(lon_unit) for lon in longitudes], axis=0) + lat_vals = np.stack([lat.to_value(lat_unit) for lat in latitudes], axis=0) + lon_q = u.Quantity(lon_vals, unit=lon_unit) + lat_q = u.Quantity(lat_vals, unit=lat_unit) # now let's deal with the fact that we may have a mixture of coords with distances and # coords with UnitSphericalRepresentations if all(targets_is_unitsphericalrep): - return SkyCoord(longitudes, latitudes, frame=frame) - elif not any(targets_is_unitsphericalrep): - return SkyCoord(longitudes, latitudes, distances, frame=frame) - else: - """ - We have a mixture of coords with distances and without. - Since we don't know in advance the origin of the frame where further transformation - will take place, it's not safe to drop the distances from those coords with them set. + return SkyCoord(lon_q, lat_q, frame=frame) + + if not any(targets_is_unitsphericalrep): + dist_unit = distances[0].unit + dist_vals = np.stack([d.to_value(dist_unit) for d in distances], axis=0) + dist_q = u.Quantity(dist_vals, unit=dist_unit) + return SkyCoord(lon_q, lat_q, dist_q, frame=frame) + + # Mixture of coords with distances and without. + # Assign large distances to UnitSphericalRepresentation objects. + filled_distances = [] + for dist, is_unitspherical in zip(distances, targets_is_unitsphericalrep): + if is_unitspherical: + fill_vals = np.broadcast_to(100.0, dist.shape if dist.shape else ()) + filled_distances.append(u.Quantity(fill_vals, u.kpc)) + else: + filled_distances.append(dist) - Instead, let's assign large distances to those objects with none. - """ - distances = [distance if distance != 1 else 100*u.kpc for distance in distances] - return SkyCoord(longitudes, latitudes, distances, frame=frame) + dist_unit = filled_distances[0].unit + dist_vals = np.stack([d.to_value(dist_unit) for d in filled_distances], axis=0) + dist_q = u.Quantity(dist_vals, unit=dist_unit) + return SkyCoord(lon_q, lat_q, dist_q, frame=frame) class SpecialObjectFlag: diff --git a/astroplan/tests/test_constraints.py b/astroplan/tests/test_constraints.py index cb599d1f..e4db748a 100644 --- a/astroplan/tests/test_constraints.py +++ b/astroplan/tests/test_constraints.py @@ -461,17 +461,7 @@ def test_caches_shapes(): targets = get_skycoord([m31, ippeg, htcas]) observer = Observer.at_site('lapalma') ac = AltitudeConstraint(min=30*u.deg) - - # When time and targets are the same size, - # they're broadcastable and grid_times_targets is ignored - assert ac(observer, targets, times, grid_times_targets=True).shape == (3,) - targets = get_skycoord([m31, ippeg]) - assert ac(observer, targets, times, grid_times_targets=True).shape == (2, 3) - - # When time and targets don't have the same size this fails with grid_times_targets=False - with pytest.raises(ValueError): - ac(observer, targets, times, grid_times_targets=False) - targets = get_skycoord([m31, ippeg, htcas]) + assert ac(observer, targets, times, grid_times_targets=True).shape == (3, 3) assert ac(observer, targets, times, grid_times_targets=False).shape == (3,) diff --git a/astroplan/tests/test_target.py b/astroplan/tests/test_target.py index 13d07145..2a9adbc4 100644 --- a/astroplan/tests/test_target.py +++ b/astroplan/tests/test_target.py @@ -80,11 +80,11 @@ def test_TLETarget(): assert isinstance(tle_target_no_observer.observer, Observer) assert abs(tle_target_no_observer.observer.location.lat) < 0.001*u.deg - tle_target_no_observer.coord(time) # Just needs to work + tle_target_no_observer.get_skycoord(time) # Just needs to work # Single time (Below Horizon) - ra_dec1 = tle_target1.coord(time) # '08h29m26.00003243s +07d31m36.65950907s' - ra_dec2 = tle_target2.coord(time) + ra_dec1 = tle_target1.get_skycoord(time) # '08h29m26.00003243s +07d31m36.65950907s' + ra_dec2 = tle_target2.get_skycoord(time) assert ra_dec1.to_string('hmsdms') == ra_dec2.to_string('hmsdms') # Comparison with the JPL Horizons System @@ -108,7 +108,7 @@ def test_TLETarget(): # Single time (Above Horizon) time_ah = Time("2023-08-02 07:20", scale='utc') - ra_dec_ah = tle_target1.coord(time_ah) # '11h19m48.53631001s +44d49m45.22194611s' + ra_dec_ah = tle_target1.get_skycoord(time_ah) # '11h19m48.53631001s +44d49m45.22194611s' ra_dec_ah_horizon_icrf = SkyCoord("11h19m49.660349s +44d49m34.65875s") assert ra_dec_ah.separation(ra_dec_ah_horizon_icrf) < 20*u.arcsec # 15.95″ @@ -125,8 +125,8 @@ def test_TLETarget(): # Should the accuracy be better than < 25*u.arcsec when compared to the JPL Horizons System? # Multiple times - ra_dec1 = tle_target1.coord(times) - ra_dec2 = tle_target2.coord(times) + ra_dec1 = tle_target1.get_skycoord(times) + ra_dec2 = tle_target2.get_skycoord(times) ra_dec_from_horizon = SkyCoord(["08h29m27.029117s +07d31m28.35610s", "06h25m46.672661s -54d32m16.77533s", "13h52m08.854291s +04d26m49.56432s", @@ -178,11 +178,11 @@ def test_TLETarget(): # InvalidTLEDataWarning and # ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" with pytest.warns(): - assert np.isnan(tle_target1.coord(time_invalid).ra) + assert np.isnan(tle_target1.get_skycoord(time_invalid).ra) # InvalidTLEDataWarning and # ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" with pytest.warns(): - assert np.isnan(tle_target1.coord(times_invalid)[2].ra) + assert np.isnan(tle_target1.get_skycoord(times_invalid)[2].ra) @pytest.mark.remote_data @@ -249,20 +249,3 @@ def test_get_skycoord_with_TLETarget(): assert mixed_output.shape == (3,) mixed_output = get_skycoord([skycoord_targed, fixed_target1, tle_target1], times) assert mixed_output.shape == (3, 4) - - # backwards_compatible - tle_output = get_skycoord(skycoord_targed, time, backwards_compatible=False) - assert tle_output.size == 1 - - tle_output = get_skycoord(skycoord_targed, times) - assert tle_output.size == 1 - tle_output = get_skycoord(skycoord_targed, times, backwards_compatible=False) - assert tle_output.shape == (4,) - - tle_output = get_skycoord([skycoord_targed, skycoord_targed], time, backwards_compatible=False) - assert tle_output.shape == (2,) - - tle_output = get_skycoord([skycoord_targed, skycoord_targed], times) - assert tle_output.shape == (2,) - tle_output = get_skycoord([skycoord_targed, skycoord_targed], times, backwards_compatible=False) - assert tle_output.shape == (2, 4) From 09dd5868ac3c53fbcde65547052215d3ceab95ab Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Thu, 12 Feb 2026 13:48:52 +0100 Subject: [PATCH 10/12] Fix code style checks IIIII --- astroplan/constraints.py | 4 ++-- astroplan/observer.py | 6 +++++- astroplan/target.py | 2 +- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 06d841ca..34604174 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -15,7 +15,7 @@ import numpy as np from astropy import table from astropy.time import Time -from astropy.coordinates import get_body, get_sun, Galactic, SkyCoord +from astropy.coordinates import get_body, get_sun, Galactic from numpy.lib.stride_tricks import as_strided # Package @@ -236,7 +236,7 @@ def __call__(self, observer, targets, times=None, The times to compute the constraint. time_range : `~astropy.time.Time` (length = 2) (optional) Lower and upper bounds on time sequence. - Only used when `times` is not provided. + Only used when ``times`` is not provided. time_grid_resolution : `~astropy.units.Quantity` Time-grid spacing. grid_times_targets : bool diff --git a/astroplan/observer.py b/astroplan/observer.py index 6cc4e9cd..c697115e 100644 --- a/astroplan/observer.py +++ b/astroplan/observer.py @@ -534,7 +534,11 @@ def _preprocess_inputs(self, time, target=None, grid_times_targets=False): if grid_times_targets: # Only ambiguous case: a single time-dependent target produces shape == time.shape # but grid mode requires a leading target axis (1, ...). - if is_single_time_dependent_target and (not target.isscalar) and (target.shape == time.shape): + if ( + is_single_time_dependent_target + and (not target.isscalar) + and (target.shape == time.shape) + ): target = target[np.newaxis, ...] # Ensure at least one targets axis for scalar targets diff --git a/astroplan/target.py b/astroplan/target.py index 0afc590c..f92f9c8d 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -498,7 +498,7 @@ def _broadcast_quantity(q, shape): dist = _broadcast_quantity(dist, times_shape) distances.append(dist) - # Convert all longitude/latitude quantities to a common unit + # Convert all longitude/latitude quantities to a common unit # and plain ndarrays before stacking (robust across units/Quantity subclasses). lon_unit = longitudes[0].unit lat_unit = latitudes[0].unit From a12c17aab875df13becf84b7db965a93489bfbc7 Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Tue, 17 Feb 2026 12:58:04 +0100 Subject: [PATCH 11/12] get_skycoord: improve performance by continuing to ignore non-scalar SkyCoord targets when checking for multiple targets --- astroplan/observer.py | 5 +---- astroplan/target.py | 7 +++---- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/astroplan/observer.py b/astroplan/observer.py index c697115e..30ba6de5 100644 --- a/astroplan/observer.py +++ b/astroplan/observer.py @@ -518,10 +518,7 @@ def _preprocess_inputs(self, time, target=None, grid_times_targets=False): return time, None # Remember whether target is a single time-dependent target - is_multiple_targets = ( - isinstance(target, (list, tuple)) or - (isinstance(target, SkyCoord) and not target.isscalar) - ) + is_multiple_targets = isinstance(target, (list, tuple)) is_target_time_dependent = ( callable(getattr(target, "get_skycoord", None)) and not hasattr(target, "coord") diff --git a/astroplan/target.py b/astroplan/target.py index f92f9c8d..11d53f75 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -422,10 +422,9 @@ def _as_coord(obj): return obj.get_skycoord(times) return obj - is_multiple_targets = ( - isinstance(targets, (list, tuple)) or - (isinstance(targets, SkyCoord) and not targets.isscalar) - ) + # Ignore non-scalar SkyCoords targets here + # e.g. from get_body/get_sun, because they represent a single target + is_multiple_targets = isinstance(targets, (list, tuple)) if not is_multiple_targets: return _as_coord(targets) From 6c4dfe58f926bdf02cfbb609c18099106e7500cc Mon Sep 17 00:00:00 2001 From: Michael Baisch Date: Mon, 16 Mar 2026 09:57:37 +0100 Subject: [PATCH 12/12] Target, Scheduling: implement review feedback * Schedule.to_table(): Re-add RA and Dec columns to avoid a breaking change, and replace broad `except Exception` handlers with narrower exception handling. * Target: Add an explicit `is_time_dependent` property to label targets that require evaluation at a specific time. --- astroplan/scheduling.py | 40 ++++++++++++++++++++++++++++++++++------ astroplan/target.py | 23 ++++++++++++++++++++++- 2 files changed, 56 insertions(+), 7 deletions(-) diff --git a/astroplan/scheduling.py b/astroplan/scheduling.py index 24220d68..d9b52796 100644 --- a/astroplan/scheduling.py +++ b/astroplan/scheduling.py @@ -11,6 +11,7 @@ from astropy import units as u from astropy.time import Time from astropy.table import Table +from astropy.coordinates import ConvertError from .utils import time_grid_from_range, stride_array from .constraints import AltitudeConstraint @@ -270,7 +271,7 @@ def _format_target_info(target): if hasattr(target, "coord"): try: return target.coord.icrs.to_string("hmsdms") - except Exception: + except ConvertError: return repr(target.coord) if hasattr(target, "alt") and hasattr(target, "az"): parts = [ @@ -279,11 +280,12 @@ def _format_target_info(target): ] if hasattr(target, "pressure"): try: - p = u.Quantity(target.pressure) - if p.to_value(u.hPa) != 0.0: - parts.append("pressure={:.3f} hPa".format(p.to_value(u.hPa))) - except Exception: + p_hpa = u.Quantity(target.pressure).to_value(u.hPa) + except (TypeError, ValueError, u.UnitConversionError): pass + else: + if p_hpa != 0.0: + parts.append("pressure={:.3f} hPa".format(p_hpa)) return ", ".join(parts) if hasattr(target, "satellite"): return ( @@ -298,6 +300,8 @@ def _format_target_info(target): durations = [] target_types = [] target_info = [] + ra = [] + dec = [] config = [] for slot in self.slots: @@ -307,6 +311,14 @@ def _format_target_info(target): durations.append(slot.duration.to(u.minute).value) target_names.append(slot.block.target.name) target_types.append(slot.block.target.__class__.__name__) + try: + ra.append(u.Quantity(slot.block.target.ra)) + except (AttributeError, NotImplementedError): + ra.append("") + try: + dec.append(u.Quantity(slot.block.target.dec)) + except (AttributeError, NotImplementedError): + dec.append("") target_info.append(_format_target_info(slot.block.target)) config.append(slot.block.configuration) elif show_transitions and slot.block: @@ -315,6 +327,8 @@ def _format_target_info(target): durations.append(slot.duration.to(u.minute).value) target_names.append("TransitionBlock") target_types.append("TransitionBlock") + ra.append("") + dec.append("") target_info.append("") changes = list(slot.block.components.keys()) if "slew_time" in changes: @@ -326,16 +340,30 @@ def _format_target_info(target): durations.append(slot.duration.to(u.minute).value) target_names.append("Unused Time") target_types.append("") + ra.append("") + dec.append("") target_info.append("") config.append("") return Table( - [target_names, start_times, end_times, durations, target_types, target_info, config], + [ + target_names, + start_times, + end_times, + durations, + ra, + dec, + target_types, + target_info, + config, + ], names=( "target", "start time (UTC)", "end time (UTC)", "duration (minutes)", + "ra", + "dec", "target type", "target info", "configuration", diff --git a/astroplan/target.py b/astroplan/target.py index 11d53f75..ba402510 100644 --- a/astroplan/target.py +++ b/astroplan/target.py @@ -60,6 +60,19 @@ def __init__(self, name=None, ra=None, dec=None, marker=None): """ raise NotImplementedError() + @property + def is_time_dependent(self): + """ + Whether this target requires evaluation at a specific time. + + Returns + ------- + is_time_dependent : bool + `True` for targets whose coordinates depend on ``obstime``, + otherwise `False`. + """ + return False + @property def ra(self): """ @@ -207,6 +220,14 @@ class TLETarget(Target): """ A target defined by TLE (Two-Line Element set) for satellites. """ + + @property + def is_time_dependent(self): + """ + Whether this target requires evaluation at a specific time. + """ + return True + def __init__(self, line1, line2, name=None, observer=None, skip_tle_check=False): """ Parameters @@ -413,7 +434,7 @@ def get_skycoord(targets, times=None): times = Time(times) def _is_time_dependent(obj): - return callable(getattr(obj, "get_skycoord", None)) and not hasattr(obj, "coord") + return isinstance(obj, Target) and obj.is_time_dependent def _as_coord(obj): if hasattr(obj, "coord"):