Select Git revision
conversions.py
TMSS-250: add tests for target rise and set
Jörn Künsemöller authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
conversions.py 13.04 KiB
from astropy.time import Time
import astropy.units
from datetime import datetime, timedelta, time as dtime
from astropy.coordinates.earth import EarthLocation
from astropy.coordinates import Angle, get_body
from astroplan.observer import Observer
import astropy.time
from functools import lru_cache
import logging
logger = logging.getLogger(__name__)
def create_astroplan_observer_for_station(station: str) -> Observer:
'''
returns an astroplan observer for object for a given station, located in the LBA center of the given station
:param station: a station name, e.g. "CS002"
:return: astroplan.observer.Observer object
'''
from lofar.lta.sip import station_coordinates
coords = station_coordinates.parse_station_coordinates()["%s_LBA" % station.upper()]
location = EarthLocation.from_geocentric(x=coords['x'], y=coords['y'], z=coords['z'], unit=astropy.units.m)
observer = Observer(location, name="LOFAR", timezone="UTC")
return observer
# default angle to the horizon at which the sunset/sunrise starts and ends, as per LOFAR definition.
SUN_SET_RISE_ANGLE_TO_HORIZON = Angle(10, unit=astropy.units.deg)
# default n_grid_points; higher is more precise but very costly; astropy defaults to 150, errors now can be in the minutes, increase if this is not good enough
SUN_SET_RISE_PRECISION = 30
@lru_cache(maxsize=256, typed=False) # does not like lists, so use tuples to allow caching
def timestamps_and_stations_to_sun_rise_and_set(timestamps: tuple, stations: tuple, angle_to_horizon: Angle=SUN_SET_RISE_ANGLE_TO_HORIZON) -> dict:
"""
Compute sunrise, sunset, day and night of the given stations at the given timestamps.
The day/sunrise/sunset is always on the date of the timestamp.
The night is usually the one _starting_ on the date of the time stamp, unless the given timestamp falls before sunrise, in which case it is the night _ending_ on the timestamp date.
:param timestamps: tuple of datetimes, e.g. (datetime(2020, 1, 1), datetime(2020, 1, 2))
:param stations: tuple of station names, e.g. ("CS002",)
:param angle_to_horizon: the angle between horizon and given coordinates for which rise and set times are returned
:return A dict that maps station names to a nested dict that contains lists of start and end times for sunrise, sunset, etc, on each requested date.
E.g.
{"CS002":
{ "sunrise": [{"start": datetime(2020, 1, 1, 6, 0, 0)), "end": datetime(2020, 1, 1, 6, 30, 0)},
{"start": datetime(2020, 1, 2, 6, 0, 0)), "end": datetime(2020, 1, 2, 6, 30, 0)}],
"sunset": [{"start": datetime(2020, 1, 1, 18, 0, 0)), "end": datetime(2020, 1, 1, 18, 30, 0)},
{"start": datetime(2020, 1, 2, 18, 0, 0)), "end": datetime(2020, 1, 2, 18, 30, 0)}],
"day": [{"start": datetime(2020, 1, 1, 6, 30, 0)), "end": datetime(2020, 1, 1, 18, 00, 0)},
{"start": datetime(2020, 1, 2, 6, 30, 0)), "end": datetime(2020, 1, 2, 18, 00, 0)}],
"night": [{"start": datetime(2020, 1, 1, 18, 30, 0)), "end": datetime(2020, 1, 2, 6, 0, 0)},
{"start": datetime(2020, 1, 2, 18,3 0, 0)), "end": datetime(2020, 1, 3, 6, 0, 0)}],
}
}
"""
return_dict = {}
for station in stations:
for timestamp in timestamps:
# todo: this can probably be made faster by moving the following logic to an own function with single station/timestamp as input and putting the lru_cache on there.
# This also means that we have to strip the time from the datetime. Can this be safely done?
observer = create_astroplan_observer_for_station(station)
sunrise_start = observer.sun_rise_time(time=Time(datetime.combine(timestamp.date(), dtime(12,0,0))), horizon=-angle_to_horizon, which='previous', n_grid_points=SUN_SET_RISE_PRECISION)
sunrise_end = observer.sun_rise_time(time=Time(sunrise_start), horizon=angle_to_horizon, which='next', n_grid_points=SUN_SET_RISE_PRECISION)
sunset_start = observer.sun_set_time(time=sunrise_end, horizon=angle_to_horizon, which='next', n_grid_points=SUN_SET_RISE_PRECISION)
sunset_end = observer.sun_set_time(time=sunset_start, horizon=-angle_to_horizon, which='next', n_grid_points=SUN_SET_RISE_PRECISION)
return_dict.setdefault(station, {}).setdefault("sunrise", []).append({"start": sunrise_start.to_datetime(), "end": sunrise_end.to_datetime()})
return_dict[station].setdefault("sunset", []).append({"start": sunset_start.to_datetime(), "end": sunset_end.to_datetime()})
return_dict[station].setdefault("day", []).append({"start": sunrise_end.to_datetime(), "end": sunset_start.to_datetime()})
if timestamp >= sunrise_start:
sunrise_next_start = observer.sun_rise_time(time=sunset_end, horizon=-angle_to_horizon, which='next', n_grid_points=SUN_SET_RISE_PRECISION)
return_dict[station].setdefault("night", []).append({"start": sunset_end.to_datetime(), "end": sunrise_next_start.to_datetime()})
else:
sunset_previous_end = observer.sun_set_time(time=sunrise_start, horizon=-angle_to_horizon, which='previous', n_grid_points=SUN_SET_RISE_PRECISION)
return_dict[station].setdefault("night", []).append({"start": sunset_previous_end.to_datetime(), "end": sunrise_start.to_datetime()})
return return_dict
# todo: Depending on usage patterns, we should consider refactoring this a little so that we cache on a function with a single timestamp as input. Requests with similar (but not identical) timestamps or bodies currently make no use of cached results for the subset computed in previous requests.
@lru_cache(maxsize=256, typed=False) # does not like lists, so use tuples to allow caching
def coordinates_and_timestamps_to_separation_from_bodies(angle1: float, angle2: float, direction_type: str, timestamps: tuple, bodies: tuple) -> dict:
"""
compute angular distances of the given sky coordinates from the given solar system bodies at the given timestamps (seen from LOFAR core)
:param angle1: first angle of celectial coordinates, e.g. RA
:param angle2: second angle of celectial coordinates, e.g. Dec
:param direction_type: direction_type of celectial coordinates, e.g. 'J2000'
:param timestamps: tuple of datetimes, e.g. (datetime(2020, 1, 1, 15, 0, 0), datetime(2020, 1, 1, 16, 0, 0))
:param bodies: tuple of solar system bodies, e.g. ('sun', 'moon', 'jupiter')
:return A dict that maps each body to a dict that maps the given timestamp to a separation angle from the given coordinate.
E.g.
{
"sun": {datetime(2020, 1, 1, 6, 0, 0): Angle("0.7rad"), datetime(2020, 1, 1, 7, 0, 0): Angle("0.7rad")},
"moon": {datetime(2020, 1, 1, 6, 0, 0): Angle("0.4rad"), datetime(2020, 1, 1, 7, 0, 0): Angle("0.4rad")},
"jupiter": {datetime(2020, 1, 1, 6, 0, 0): Angle("2.7rad"), datetime(2020, 1, 1, 7, 0, 0): Angle("2.7rad")}
}
"""
if direction_type == "J2000":
coord = astropy.coordinates.SkyCoord(ra=angle1, dec=angle2, unit=astropy.units.rad)
else:
raise ValueError("Do not know how to convert direction_type=%s to SkyCoord" % direction_type)
return_dict = {}
for body in bodies:
location = create_astroplan_observer_for_station("CS002").location
for timestamp in timestamps:
# get body coords at timestamp
body_coord = get_body(body=body, time=astropy.time.Time(timestamp), location=location)
angle = coord.separation(body_coord)
return_dict.setdefault(body, {})[timestamp] = angle
return return_dict
# default angle above horizon, above which the target it reporte as 'up'
TARGET_SET_RISE_ANGLE_TO_HORIZON = Angle(0, unit=astropy.units.deg) # if default should be non-zero, should we include it explicitly in response?
# default n_grid_points; higher is more precise but very costly; astropy defaults to 150, note that errors can be in the minutes with a lower values
TARGET_SET_RISE_PRECISION = 150
@lru_cache(maxsize=256, typed=False) # does not like lists, so use tuples to allow caching
def coordinates_timestamps_and_stations_to_target_rise_and_set(angle1: float, angle2: float, direction_type: str, timestamps: tuple, stations: tuple, angle_to_horizon: Angle=TARGET_SET_RISE_ANGLE_TO_HORIZON) -> dict:
"""
Compute rise and set times of the given coordinates above the provided horizon, for each given station and timestamp.
The set time is always the one following the provided timestamp.
This implies that if the target is up at a given timestamp, the surrounding rise and set times are returned.
Otherwise both rise and set times follow the timestamp.
:param angle1: first angle of celectial coordinates, e.g. RA
:param angle2: second angle of celectial coordinates, e.g. Dec
:param direction_type: direction_type of celectial coordinates, e.g. 'J2000'
:param timestamps: tuple of datetimes, e.g. (datetime(2020, 1, 1), datetime(2020, 1, 2))
:param stations: tuple of station names, e.g. ("CS002",)
:param angle_to_horizon: the angle between horizon and given coordinates for which rise and set times are returned
:return A dict that maps station names to a list of dicts with rise and set times for each requested date.
E.g.
{"CS002": [{"rise": datetime(2020, 1, 1, 4, 0, 0), "set": datetime(2020, 1, 1, 11, 0, 0)},
{"rise": datetime(2020, 1, 2, 4, 0, 0), "set": datetime(2020, 1, 2, 11, 0, 0)}]
}
"""
if direction_type == "J2000":
coord = astropy.coordinates.SkyCoord(ra=angle1, dec=angle2, unit=astropy.units.rad)
else:
raise ValueError("Do not know how to convert direction_type=%s to SkyCoord" % direction_type)
return_dict = {}
for station in stations:
for timestamp in timestamps:
# todo: this can probably be made faster by moving the following logic to an own function with single station/timestamp as input and putting the lru_cache on there.
observer = create_astroplan_observer_for_station(station)
target_set = observer.target_set_time(target=coord, time=Time(timestamp), horizon=angle_to_horizon, which='next', n_grid_points=TARGET_SET_RISE_PRECISION)
target_rise = observer.target_rise_time(target=coord, time=Time(target_set), horizon=angle_to_horizon, which='previous', n_grid_points=TARGET_SET_RISE_PRECISION)
return_dict.setdefault(station, []).append({"rise": target_rise.to_datetime(), "set": target_set.to_datetime()})
return return_dict
def local_sidereal_time_for_utc_and_station(timestamp: datetime = None,
station: str = 'CS002',
field: str = 'LBA',
kind: str = "apparent"):
"""
calculate local sidereal time for given utc time and station
:param timestamp: timestamp as datetime object
:param station: station name
:param field: antennafield, 'LBA' or 'HBA'
:param kind: 'mean' or 'apparent'
:return:
"""
from lofar.lta.sip import station_coordinates
if timestamp is None:
timestamp = datetime.utcnow()
station_coords = station_coordinates.parse_station_coordinates()
field_coords = station_coords["%s_%s" % (station, field)]
location = EarthLocation.from_geocentric(x=field_coords['x'], y=field_coords['y'], z=field_coords['z'], unit=astropy.units.m)
return local_sidereal_time_for_utc_and_longitude(timestamp=timestamp, longitude=location.lon.to_string(decimal=True), kind=kind)
def local_sidereal_time_for_utc_and_longitude(timestamp: datetime = None,
longitude: float = 6.8693028,
kind: str = "apparent"):
"""
:param timestamp: timestamp as datetime object
:param longitude: decimal longitude of observer location (defaults to CS002 LBA center)
:param kind: 'mean' or 'apparent'
:return:
"""
if timestamp is None:
timestamp = datetime.utcnow()
t = Time(timestamp, format='datetime', scale='utc')
return t.sidereal_time(kind=kind, longitude=longitude)
def antennafields_for_antennaset_and_station(antennaset:str, station:str) -> list:
"""
convert an antennaset to a list of antennafields
:param antennaset: A string identifier for an antennaset, like 'HBA_DUAL'
:param station: A string identifier for a station, like 'CS001'
:return: a list of antennafields that the station uses for the given antennaset ['HBA0', 'HBA1']
"""
if antennaset.startswith('LBA'):
fields = ['LBA']
elif antennaset.startswith('HBA') and not station.startswith('CS'):
fields = ['HBA']
elif antennaset.startswith('HBA_DUAL'):
fields = ['HBA0', 'HBA1']
elif antennaset.startswith('HBA_ZERO'):
fields = ['HBA0']
elif antennaset.startswith('HBA_ONE'):
fields = ['HBA1']
else:
raise ValueError('Cannot determine antennafields for station=%s antennaset=%s' % (station, antennaset))
return fields