Browse Source

Improved documentation and stability of the code

pull/1/head
Benito Marcote 1 year ago
parent
commit
ab0e5fcde2
  1. 5
      CHANGES.txt
  2. 38
      app.py
  3. 20
      doc/notes.md
  4. 5
      setup.py
  5. 28
      tests/test.py
  6. 4
      vlbiplanobs/graphical_elements.py
  7. 544
      vlbiplanobs/observation.py
  8. 486
      vlbiplanobs/stations.py

5
CHANGES.txt

@ -1 +1,4 @@
v1.0, 2020-10-23 -- First stable release.
v1.0.1 2020-10-26 -- Source coordinates from server processed seamlessly with extra empty spaces between RA and DEC.
-- Improved documentation in functions.
-- Improved stability of the code with additional checks and explicit type hints.
v1.0 2020-10-23 -- First stable release.

38
app.py

@ -35,15 +35,15 @@ from astropy import units as u
## THIS WILL NEED TO GO AWAY IN THE NEW VERSION OF ASTROPY, WHICH IS STILL NOT
## SUPPORTED BY THE CURRENT VERSION OF ASTROPLAN
# Tweak to not let astroplan crashing...
from astropy.utils.data import clear_download_cache
from astropy.utils import iers
clear_download_cache() # to be sure it is really working
iers.conf.auto_download = False
iers.conf.iers_auto_url = None
iers.conf.auto_max_age = None
iers.conf.remote_timeout = 100.0
iers.conf.download_cache_lock_attempts = 10
# from astropy.utils.data import clear_download_cache
# from astropy.utils import iers
# clear_download_cache() # to be sure it is really working
#
# iers.conf.auto_download = False
# iers.conf.iers_auto_url = None
# iers.conf.auto_max_age = None
# iers.conf.remote_timeout = 100.0
# iers.conf.download_cache_lock_attempts = 10
from astroplan import FixedTarget
@ -90,7 +90,7 @@ default_datarates = {'EVN': 2048, 'e-EVN': 2048, 'eMERLIN': 4096, 'LBA': 1024, '
# Safety check that all these antennas are available in the file
for a_array in default_arrays:
for a_station in default_arrays[a_array]:
assert a_station in all_antennas.keys()
assert a_station in all_antennas.codenames
doc_files = {'About this tool': 'doc-contact.md',
'About the antennas': 'doc-antennas.md',
@ -176,7 +176,7 @@ def convert_colon_coord(colon_coord):
def alert_message(message, title="Warning!"):
"""Produces an alert-warning message.
message can be either a string or a list with different string/dash components.
'message' can be either a string or a list with different string/dash components.
"""
if type(message) == str:
return [html.Br(), \
@ -189,7 +189,8 @@ def alert_message(message, title="Warning!"):
def update_sensitivity(obs):
"""Given the observation, it sets the text about all properties of the observation.
"""Given the observation, it sets the text for all summary cards
with information about the observation.
"""
cards = []
# The time card
@ -520,6 +521,8 @@ app.layout = html.Div([
@app.callback(Output('onsourcetime-label', 'children'),
[Input('onsourcetime', 'value')])
def update_onsourcetime_label(onsourcetime):
"""Keeps the on-source time label updated with the value selected by the user.
"""
return f"% of on-target time ({onsourcetime}%)"
@ -557,6 +560,12 @@ def select_antennas(selected_band, selected_networks, is_eEVN):
[Input('starttime', 'date'), Input('starthour', 'value'),
Input('endtime', 'date'), Input('endhour', 'value')])
def check_obstime(starttime, starthour, endtime, endhour):
"""Verify the introduced times/dates for correct values.
Once the user has introduced all values for the start and end of the observation,
it guarantees that they have the correct shape:
- the end of the observation is afte the start of the observation.
- The total observing length is less than five days (value chosen for computational reasons).
"""
times = [None, None]
if None not in (starttime, starthour):
times[0] = Time(dt.strptime(f"{starttime} {starthour}", '%Y-%m-%d %H:%M'))
@ -578,12 +587,15 @@ def check_obstime(starttime, starthour, endtime, endhour):
@app.callback(Output('error_source', 'children'),
[Input('source', 'value')])
def get_source(source_coord):
"""Verifies that the introduced source coordinates have a right format.
If they are correct, it does nothing. If they are incorrect, it shows an error label.
"""
if source_coord != 'hh:mm:ss dd:mm:ss' and source_coord != None:
try:
dummy_target = observation.Source(convert_colon_coord(source_coord), 'Source')
return ''
except ValueError as e:
return "Use 'hh:mm:ss dd:mm:ss' format"
return "Use 'hh:mm:ss dd:mm:ss' or 'XXhXXmXXs XXdXXmXXs' format"
else:
return dash.no_update

20
doc/notes.md

@ -7,7 +7,23 @@
[ ] Update the 'install_requires' in setup.py (line 54).
[ ] Move the app.py into bin/ folder, so it is directly called from terminal.
https://python-packaging.readthedocs.io/en/latest/command-line-scripts.html
[ ] Add __all__ = ['', ...] line in all modules (functions/classes to be imported with *).
[ ] Add __all__ = ['', ...] line in all modules (functions/classes to be imported with \*).
[ ] Target field allowing source name instead of only coordinates.
[ ] Resolution with different robust weighting.
## Dependencies
App.py
numpy
dash
dash_core_components
dash_html_components
dash_bootstrap_components
plotly
astropy
astroplan
@ -108,7 +124,7 @@ Source
+ __init__(coordinates : coord.SkyCoord, name=None)
Observation
- target : FixedTarget
- target : Source
- times : Time
- gstimes : Longitude (hourangle)
- duration : Time

5
setup.py

@ -10,7 +10,7 @@ setup(name='vlbiplanobs',
# Versions should comply with PEP 440:
# https://www.python.org/dev/peps/pep-0440/
# [N!]N(.N)*[{a|b|rc}N][.postN][.devN]
version='1.0',
version='1.0.1',
# one-line description or tagline of what your project does
description='Planner for VLBI observations',
# "Description" metadata field
@ -54,8 +54,7 @@ setup(name='vlbiplanobs',
#
# For an analysis of "install_requires" vs pip's requirements files see:
# https://packaging.python.org/en/latest/requirements.html
# TODO
# install_requires=[''],
install_requires=['numpy', 'astropy>=4.0.2', 'astroplan>=0.7'],
# If there are data files included in your packages that need to be
# installed, specify them here.

28
tests/test.py

@ -7,16 +7,15 @@ Note that this is not a full test...
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
# import matplotlib.pyplot as plt
# import matplotlib.dates as mdates
from astropy import coordinates as coord
from astropy import units as u
from astropy.time import Time
from astropy.io import ascii
from src import freqsetups as fs
from src import stations
from src import functions as fx
from src import observation
from vlbiplanobs import freqsetups as fs
from vlbiplanobs import stations
from vlbiplanobs import observation
@ -32,7 +31,7 @@ obs.channels = 32
obs.polarizations = 2
obs.inttime = 2
all_stations = fx.get_stations_from_configfile(f"data/stations_catalog.inp")
all_stations = stations.Stations.get_stations_from_configfile(f"data/stations_catalog.inp")
def get_selected_antennas(list_of_selected_antennas):
"""Given a list of antenna codenames, it returns a Stations object containing
@ -53,15 +52,15 @@ obs.stations = get_selected_antennas(evn6)
elevs = obs.elevations()
srcup = obs.is_visible()
fig, ax = plt.subplots()
# fig, ax = plt.subplots()
for ant in elevs:
ax.plot(obs.times.datetime[srcup[ant]], elevs[ant][srcup[ant]].deg, '-', label=ant)
# for ant in elevs:
# ax.plot(obs.times.datetime[srcup[ant]], elevs[ant][srcup[ant]].deg, '-', label=ant)
ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
ax.set_xlabel(f"Time (UTC) - {obs.times.datetime[0].strftime('%d-%m-%Y')}")
# ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
# ax.set_xlabel(f"Time (UTC) - {obs.times.datetime[0].strftime('%d-%m-%Y')}")
@ -92,9 +91,8 @@ def gst2t(t):
ax_gst = ax.secondary_xaxis("top", functions=(t2gst, gst2t))
ax_gst.cla()
# ax_gst = ax.secondary_xaxis("top", functions=(t2gst, gst2t))
# ax_gst.cla()

4
vlbiplanobs/graphical_elements.py

@ -108,9 +108,9 @@ def summary_card_antennas(app, obs):
if len(ants_up[an_ant][0]) == 0:
ant_no_obs.append(an_ant)
ant_text = ', '.join([ant for ant in obs.stations.keys() if ant not in ant_no_obs]) + '.'
ant_text = ', '.join([ant for ant in obs.stations.codenames if ant not in ant_no_obs]) + '.'
ant_text = []
for ant in obs.stations.keys():
for ant in obs.stations.codenames:
if ant not in ant_no_obs:
ant_text += tooltip_card(antenna_card(app, obs.stations[ant]),
idname=f"basel-{ant}", trigger=ant, placement='top')

544
vlbiplanobs/observation.py

@ -1,56 +1,124 @@
import numpy as np
from astropy import units as u
from astropy import coordinates as coord
from astropy.time import Time
from astropy.time import Time, TimeDelta
from astroplan import FixedTarget
import vlbiplanobs
from vlbiplanobs.stations import Stations
"""Defines an observation with given network(s), a target source,
time range and the observing band.
"""Defines an observation, which basically consist of a given network of stations,
observing a target source for a given time range and at an observing band.
"""
class SourceNotVisible(Exception):
"""Exception produced when a given target source cannot be observed for any
antenna in the network.
"""
pass
class Source(FixedTarget):
"""Defines a source with some coordinates and a name.
Inputs
- coordinates : str
In a format recognized by astropy.coordinates.SkyCoord
(e.g. XXhXXmXXs XXdXXmXXs)
- name : str (optional)
Name associated to the source
"""Defines a target source located at some coordinates and with a given name.
"""
def __init__(self, coordinates, name=None):
def __init__(self, coordinates: str, name: str = None):
"""Initializes a Source object.
Inputs
- coordinates : str
Coordinates of the target source in a str format recognized by
astropy.coordinates.SkyCoord (e.g. XXhXXmXXs XXdXXmXXs).
J2000 coordinates are assumed.
- name : str [OPTIONAL]
Name associated to the source. By default is None.
"""
super().__init__(coord.SkyCoord(coordinates), name)
class Observation(object):
def __init__(self, target=None, times=None, band=None, datarate=None,subbands=None,
channels=None, polarizations=None, inttime=None, ontarget=1.0,
stations=None, bits=2):
if target is not None:
self.target = target
if times is not None:
"""Defines an observation of a single target source with a given network of stations.
The observation can be set up in incremental steps (i.e. not all inputs are necessary
at initialization time). Depending on the functions that you want to execute not all
attributes may be required.
An Observation allows the user to determine the elevation of the target source for
each antenna in the network, the estimated rms noise level that would reach in such
observation, or the resolution of the resulting images, assuming a standard neutral
robust weighting.
"""
def __init__(self, target: Source = None, times: Time = None, band: str = None,
datarate=None, subbands: int = None, channels: int = None,
polarizations: int = None, inttime=None, ontarget: float = 1.0,
stations: Stations = None, bits: int = 2):
"""Initializes an observation.
Note that you can initialize an empty observation at this stage and add the
information for the different attributes later. However, you may raise exception
if running methods that require some of the unset attributes.
Inputs
- target : vlbiplanobs.observation.Source
Target source to be observed.
- times : astropy.time.Time
An array of times defining the duration of the observation. That is,
the first time will define the start of the observation and the last one
will be the end of the observation. Note that the higher time resolution
in `times` the more precise values you will obtain for antenna source
visibility or determination of the rms noise levels. However, that will
also imply a longer computing time. Steps of ~5-15 min seem appropriate
for typical VLBI observations.
- band : str
Observing band to conduct the observation. Note that while this is a
free-format string, it should match the type used when defining the bands
at which each station can observe, and the default is the str format `XXcm`
where XX represents the wavelength to observe in cm units.
You can always check the available bands in `{Stations object}.observing_bands`.
- datarate : int or astropy.units.Quantity
Data rate for each antenna. It assumes that all antennas will run at the same
data rate, which may not be true. If an int is introduce, it will be assumed to
be given in Mbit/s. Otherwise, you can specify a Quantity with compatible units.
- subbands : int
Number of subbands in which the total bandwidth of the observation will be divided
during correlation.
- channels : int
Number of channels for each subband to be created during correlation.
- polarizations : int (1, 2, 4)
Number of polarizations to record in the observation. Three values are only possible
0: single polarization.
2: dual polarization (only direct terms: RR and LL, or XX and YY, are kept).
4: full polarization (crossed terms are kept: RR, LL, RL, LR; or XX, YY, XY, YX).
- inttime : float/int or astropy.units.Quantity
Integration time (time resolution) that the final, correlated, data will show.
If no units provided, seconds are assumed.
- ontarget : float
Fraction of the total observing time (end time minus start time) that will be
spent on the target source. Note that in a typical VLBI observation only a fraction
of the total observing time will end up in on-target source, commonly between 0.4-0.8
(~40-80% of the time). This has an effect on the determination of the final rms noise level.
- stations : vlbiplanobs.stations.Stations
Network of stations that will participate in the given observation.
- bits : int
Number of bits at which the data have been recorded (sampled). A typical VLBI observation is
almost always recorded with 2-bit data.
"""
self.target = target
# Because otherwise gstimes is not initialized
if times is None:
self._times = None
self._gstimes = None
else:
self.times = times
if band is not None:
self.band = band
if datarate is not None:
self.datarate = datarate
if subbands is not None:
self.subbands = subbands
if channels is not None:
self.channels = channels
if polarizations is not None:
self.polarizations = polarizations
if inttime is not None:
self.inttime = inttime
self.band = band
self.datarate = datarate
self.subbands = subbands
self.channels = channels
self.polarizations = polarizations
self.inttime = inttime
if stations is not None:
self.stations = stations
else:
self.stations = vlbiplanobs.stations.Stations('empty', [])
self.stations = Stations('empty', [])
self.bitsampling = bits
self.ontarget_fraction = ontarget
@ -61,25 +129,39 @@ class Observation(object):
@property
def target(self):
def target(self) -> Source:
"""Returns the target source to be observed during the current observation.
It can return None if the target has not been set yet, showing a warning.
"""
if self._target is None:
print("WARNING: 'target' source not set yet but used in 'Observation'.")
return self._target
@target.setter
def target(self, new_target):
assert isinstance(new_target, Source)
def target(self, new_target: Source):
assert isinstance(new_target, Source) or new_target is None
self._target = new_target
# Resets all parameters than depend on the source coordinates
self._uv_baseline = None
self._uv_array = None
self._rms = None
self._synth_beam = None
@property
def times(self):
def times(self) -> Time:
"""Returns the times when the observation runs as an astropy.time.Time object.
It can return None if the times have not been set yet, showing a warning.
"""
if self._times is None:
print("WARNING: 'times' not set yet but used in 'Observation'.")
return self._times
@times.setter
def times(self, new_times):
assert isinstance(new_times, Time)
assert isinstance(new_times, Time) and len(new_times) >= 2
self._times = new_times
self._gstimes = self._times.sidereal_time('mean', 'greenwich')
self._uv_baseline = None
@ -87,129 +169,234 @@ class Observation(object):
self._rms = None
self._synth_beam = None
@property
def gstimes(self):
def gstimes(self) -> coord.angles.Longitude:
"""Returns the GST times when the observation runs as an astropy.coordinates.angles.Longitude
object (meaning in hourangle units).
It can return None if the times have not been set yet, showing a warning.
"""
if self._gstimes is None:
print("WARNING: 'times' not set yet but used in 'Observation'.")
return self._gstimes
@property
def duration(self):
return self.times[-1]-self.times[0]
def duration(self) -> u.Quantity:
"""Returns the total duration of the observation.
It raises AttributeError if the attribute 'times' has not been set yet.
"""
if self.times is None:
raise AttributeError("'times' in 'Observation' has not been set yet.")
return (self.times[-1] - self.times[0]).to(u.h)
@property
def band(self):
def band(self) -> str:
"""Returns the observing band at which the observation will be conducted.
It can return None if the band has not been set yet, showing a warning.
"""
if self._band is None:
print("WARNING: 'band' not set yet but used in 'Observation'.")
return self._band
@band.setter
def band(self, new_band):
assert isinstance(new_band, str) and 'cm' in new_band
assert (isinstance(new_band, str) and 'cm' in new_band) or (new_band is None)
self._band = new_band
self._uv_baseline = None
self._uv_array = None
self._rms = None
self._synth_beam = None
@property
def wavelength(self):
"""Returns the central wavelength of the observations with units
def wavelength(self) -> u.Quantity:
"""Returns the central wavelength of the observation.
"""
assert self.band is not None
return float(self.band.replace('cm',''))*u.cm
@property
def frequency(self):
"""Returns the central frequency of the observations with units
def frequency(self) -> u.Quantity:
"""Returns the central frequency of the observations.
"""
assert self.band is not None
return 30*u.GHz/self.wavelength.to(u.cm).value
@property
def datarate(self):
"""Datarate in Mbps
def datarate(self) -> u.Quantity:
"""Retuns the data rate (per station) used at which the observation is conducted.
It can return None if the data rate has not been set yet, showing a warning.
"""
if self._datarate is None:
print("WARNING: 'datarate' not set yet but used in 'Observation'.")
return self._datarate
@datarate.setter
def datarate(self, new_datarate):
"""If no units provided, Mbps assumed.
"""Sets the data rate used at each station during the observation.
Inputs
- new_datarate : int or astropy.units.Quantity
If no units provided, Mbps assumed.
"""
if isinstance(new_datarate, int):
if new_datarate is None:
self._datarate = None
elif isinstance(new_datarate, int):
self._datarate = new_datarate*u.Mbit/u.s
elif isinstance(new_datarate, u.Quantity):
self._datarate = new_datarate.to(u.Mbit/u.s)
else:
raise ValueError(f"Unknown type for {new_datarate} (int/Quantity(bit/s) expected)")
raise ValueError(f"Unknown type for datarate {new_datarate}" \
"(int or astropy.units.Quantity (~bit/s) expected)")
self._rms = None
@property
def subbands(self):
def subbands(self) -> int:
"""Returns the number of subbands (also known as intermediate frequencies for AIPS users)
in which the total bandwidth of the observation will be divided during correlation.
It can return None if the number of subbands has not been set yet, showing a warning.
"""
if self._subbands is None:
print("WARNING: 'subbands' not set yet but used in 'Observation'.")
return self._subbands
@subbands.setter
def subbands(self, n_subbands):
assert isinstance(n_subbands, int)
def subbands(self, n_subbands: int):
assert (isinstance(n_subbands, int) and n_subbands > 0) or n_subbands is None
self._subbands = n_subbands
@property
def channels(self):
def channels(self) -> int:
"""Returns the number of channels in which each subband will be divided during correlation.
It can return None if the number of channels has not been set yet, showing a warning.
"""
if self._channels is None:
print("WARNING: 'channels' not set yet but used in 'Observation'.")
return self._channels
@channels.setter
def channels(self, n_channels):
assert isinstance(n_channels, int)
def channels(self, n_channels: int):
assert (isinstance(n_channels, int) and n_channels > 0) or (n_channels is None)
self._channels = n_channels
@property
def polarizations(self):
def polarizations(self) -> int:
"""Returns the number of polarizations that will be stored in the final data.
It can return None if the number of polarizations has not been set yet, showing a warning.
"""
if self._polarizations is None:
print("WARNING: 'polarizations' not set yet but used in 'Observation'.")
return self._polarizations
@polarizations.setter
def polarizations(self, n_pols):
assert n_pols in (1, 2, 4)
def polarizations(self, n_pols: int):
assert (n_pols in (1, 2, 4)) or (n_pols is None)
self._polarizations = n_pols
@property
def inttime(self):
"""Integration time during the observation.
def inttime(self) -> u.Quantity:
"""Returns the integration time used when correlating the observation as an astropy.units.Quantity.
It can return None if the integration time has not been set yet, showing a warning.
"""
if self._inttime is None:
print("WARNING: 'inttime' not set yet but used in 'Observation'.")
return self._inttime
@inttime.setter
def inttime(self, new_inttime):
if isinstance(new_inttime, float) or isinstance(new_inttime, int):
"""Sets the integration time of the observation.
Inputs
- new_inttime float/int or astropy.units.Quantity.
If no units provided, seconds are assumed.
"""
if new_inttime is None:
self._inttime = None
elif isinstance(new_inttime, float) or isinstance(new_inttime, int):
self._inttime = new_inttime*u.s
elif isinstance(new_inttime, u.Quantity):
self._inttime = new_inttime.to(u.s)
else:
raise ValueError(f"Unknown type for {new_inttime} (float/int/Quantity(s) expected)")
raise ValueError(f"Unknown type for 'inttime' {new_inttime} (float/int/Quantity(~seconds) expected)")
@property
def ontarget_fraction(self):
"""Fraction of the total observing time spent on the target source
def ontarget_fraction(self) -> float:
"""Fraction of the total observing time spent on the target source.
It can return None if the ontarget_fraction has not been set yet, showing a warning.
"""
if self._ontarget is None:
print("WARNING: 'ontarget_fraction' not set yet but used in 'Observation'.")
return self._ontarget
@ontarget_fraction.setter
def ontarget_fraction(self, ontarget):
assert 0.0 < ontarget <= 1.0
def ontarget_fraction(self, ontarget: float):
assert (0.0 < ontarget <= 1.0) or (ontarget is None)
self._ontarget = ontarget
self._rms = None
@property
def ontarget_time(self):
def ontarget_time(self) -> u.Quantity:
"""Total time spent on the target source during the observation.
It can return None if the ontarget_fraction and duration have not been set yet, showing a warning.
"""
if (self._ontarget is None):
raise AttributeError("'ontarget_time' in 'Observation' has not been set yet.")
return self.duration*self.ontarget_fraction
@property
def bandwidth(self):
def bandwidth(self) -> u.Quantity:
"""Returns the total bandwidth of the observation.
It raises AttributeError if the attributes 'polarizations', 'datarate', or 'bitsampling'
have not been set yet.
"""
if None in (self.polarizations, self.datarate, self.bitsampling):
raise AttributeError("'polarizations', 'datarate', and 'bitsampling' in 'Observation'" \
"have not been set yet.")
pols = self.polarizations % 3 + self.polarizations // 3 # Either 1 or 2
return (self.datarate/(pols*self.bitsampling*2)).to(u.MHz)
@property
def bitsampling(self):
def bitsampling(self) -> u.Quantity:
"""Returns the bit sampling at which the data are recorded during the observation.
"""
return self._bitsampling
@bitsampling.setter
def bitsampling(self, new_bitsampling):
"""Sets the bit sampling of the observation.
Inputs
- new_bitsampling : int/float or astropy.units.Quantity
If no units provided, bits are assumed.
"""
if isinstance(new_bitsampling, float) or isinstance(new_bitsampling, int):
self._bitsampling = new_bitsampling*u.bit
elif isinstance(new_inttime, u.Quantity):
@ -217,47 +404,84 @@ class Observation(object):
else:
raise ValueError(f"Unknown type for {new_bitsampling} (float/int/Quantity(bit) expected)")
@property
def stations(self):
def stations(self) -> Stations:
"""Returns the network of stations 'Stations' that will participate in this observation
observing the target source.
"""
return self._stations
@stations.setter
def stations(self, new_stations):
assert isinstance(new_stations, vlbiplanobs.stations.Stations)
def stations(self, new_stations: Stations):
assert isinstance(new_stations, Stations)
self._stations = new_stations
self._uv_baseline = None
self._uv_array = None
self._rms = None
self._synth_beam = None
def elevations(self):
"""Returns the target elevations for all stations along the observation.
def elevations(self) -> dict:
"""Returns the elevation of the target source for each stations participating in the observation
for all the given observing times.
Returns
elevations : dict
Dictionary where they keys are the station code names, and the values will be
an astropy.coordinates.angles.Latitude object with the elevation at each observing time.
"""
elevations = {}
for a_station in self.stations:
elevations[a_station.codename] = a_station.elevation(self.times, self.target)
return elevations
def altaz(self):
"""Returns the target altaz for all stations along the observation.
def altaz(self) -> dict:
"""Returns the altitude/azimuth of the target source for each stations participating
in the observation for all the given observing times.
Returns
altaz : dict
Dictionary where they keys are the station code names, and the values will be
an astropy.coordinates.SkyCoord object with the altitude and azimuth
of the source at each observing time.
"""
aa = {}
for a_station in self.stations:
aa[a_station.codename] = a_station.altaz(self.times, self.target)
return aa
def is_visible(self):
"""Returns when the target is visible for all stations for each time in the observation.
def is_visible(self) -> dict:
"""Returns whenever the target source is visible for each station for each time
of the observation.
Returns
is_visible : dict
Dictionary where they keys are the station code names, and the values will be
a tuple containing a numpy array with the indexes in the `Observation.times`
array with the times where the target source can be observed by the station.
In this sense, you can e.g. call obs.times[obs.is_visible[a_station_codename]]
to get such times.
"""
iv = {}
for a_station in self.stations:
iv[a_station.codename] = a_station.is_visible(self.times, self.target)
return iv
def longest_baseline(self):
def longest_baseline(self) -> tuple:
"""Returns the longest baseline in the observation.
It retuns the tuple ((ant1,ant2), length)
where `ant1,ant2` are the antennas in such baseline, and `length` its length.
Returns
- ('{ant1}-{ant2}', length) : tuple
- '{ant1}-{ant2}' : str
Composed by the codenames of the two antennas (ant1, ant2) conforming the longest baseline.
- length : astropy.units.Quantity
The projected length of the baseline as seen from the target source position.
"""
uv = self.get_uv_baseline()
longest_bl = {'bl': '', 'value': None}
@ -269,10 +493,16 @@ class Observation(object):
return longest_bl['bl'], longest_bl['value']*self.wavelength
def shortest_baseline(self):
def shortest_baseline(self) -> tuple:
"""Returns the shortest baseline in the observation.
It retuns the tuple ((ant1,ant2), length)
where `ant1,ant2` are the antennas in such baseline, and `length` its length.
Returns
- ('{ant1}-{ant2}', length) : tuple
- '{ant1}-{ant2}' : str
Composed by the codenames of the two antennas (ant1, ant2) conforming the shortest baseline.
- length : astropy.units.Quantity
The projected length of the baseline as seen from the target source position.
"""
uv = self.get_uv_baseline()
shortest_bl = {'bl': '', 'value': None}
@ -285,32 +515,66 @@ class Observation(object):
return shortest_bl['bl'], shortest_bl['value']*self.wavelength
def bandwidth_smearing(self):
def bandwidth_smearing(self) -> u.Quantity:
"""Returns the bandwidth smearing expected for the given observation.
The peak response to a point target source decreases at positions farther away from the
pointing (correlated) sky position due to the frequency averaging performed in the data.
This function returns the angular separation at which the bandwidth smearing produces
a reduction of a 10% in the response of the telescope. The field of view should then
be limited to this range to avoid significant loses.
"""
# TODO: Check if with units it works
return ((49500*u.arcsec*u.MHz*u.km)*self.channels/ \
(self.longest_baseline()[1]*self.bandwidth/self.subbands)).to(u.arcsec)
def time_smearing(self):
def time_smearing(self) -> u.Quantity:
"""Returns the time smearing expected for the given observation.
The peak response to a point target source decreases at positions farther away from the
pointing (correlated) sky position due to the time averaging performed in the data.
This function returns the angular separation at which the time smearing produces
a reduction of a 10% in the response of the telescope. The field of view should then
be limited to this range to avoid significant loses.
"""
# TODO: Check if with units it works
return ((18560*u.arcsec*u.km*u.s/u.cm)* \
(self.wavelength/(self.longest_baseline()[1]*self.inttime))).to(u.arcsec)
def datasize(self):
"""Returns the expected size for the output FITS files.
def datasize(self) -> u.Quantity:
"""Returns the expected size for the output FITS IDI files.
A regular observation with the European VLBI Network is stored in FITS IDI files,
typically several 2-GB files. This function provides an estimation of the total
size for these stored files.
Note that this function does not take into account down times for the different
stations. The provided value will thus always be un upper-limit for the real, final,
value.
"""
temp = len(self.stations)**2*((self.times[-1]-self.times[0])/self.inttime).decompose()
temp = len(self.stations)**2*(self.duration/self.inttime).decompose()
temp *= self.polarizations*self.subbands*self.channels
temp *= 1.75*u.GB/(131072*3600)
return temp.to(u.GB)
return temp*1.75*u.GB/(131072*3600)
def thermal_noise(self) -> u.Quantity:
"""Returns the expected rms thermal noise for the given observation.
def thermal_noise(self):
"""Returns the expected thermal noise for the given observation
Each antenna has a different sensitivity for the observing band (established from
the SEFD value). This function computes the available baselines at each timestamp
and estimates the final thermal noise reached when integrating the whole observation.
Note that this function takes into account when each station can observe the source,
but does not take into account sensitivity drops doe to external factors like e.g.
low elevations of the source. The provided thermal noise is also assumed when no
weighting is applied to the data when imaging. The thermal noise can thus be a bit
lower if a natural robust weighting is used, but slightly higher if an uniform
roubst is used.
"""
# As the determination is computationally heavy, if no parameters have been updated
# the nit returns the stored value from a previous run.
if self._rms is not None:
return self._rms
@ -318,7 +582,6 @@ class Observation(object):
visible = self.is_visible()
for i,stat in enumerate(self.stations):
main_matrix[:,i][visible[stat.codename]] = stat.sefd(self.band)
# For each timestamp
# Determines the noise level for each time stamp.
temp = 0.0
for i,ti in enumerate(main_matrix[:-1]):
@ -332,19 +595,27 @@ class Observation(object):
raise SourceNotVisible('No single baseline can observe the source.')
temp = 1.0/np.sqrt(temp*self.ontarget_fraction)
# TODO: fix units problem.
self._rms = ((1.0/0.7)*temp/np.sqrt(self.datarate.to(u.bit/u.s).value/2))*u.Jy
return self._rms
def get_uv_baseline(self):
"""Returns the uv values for each baseline and each timestamp when the source
def get_uv_baseline(self) -> dict:
"""Returns the (u, v) values for each baseline and each timestamp for which the source
is visible.
It returns a dictionary containing the uv values in lambda units
for each baseline as key.
Complex conjugates are not provided.
It may raise the exception SourceNotVisible if no antennas can observe the source
at all during the observation.
Returns
- {'{ant1}-{ant2}': uv_data} : dict
- '{ant1}-{ant2}' : str
The keys of the dictionary are the baselines forming the array, where ant1 and ant2
are the code names of each station.
- uv_data : astropy.units.Quantity
A 2-d array with all (u, v) values for each timestamp in `times` when the source is
visible for the given baseline. (u,v) are given in lambda units.
Note that complex conjugate values are not provided.
Exceptions
- It may raise the exception SourceNotVisible if no baselines can observe the source
at all during the observation.
"""
if self._uv_baseline is not None:
return self._uv_baseline
@ -365,7 +636,7 @@ class Observation(object):
bl_names.append("{}-{}".format(self.stations[i].codename,
self.stations[j].codename))
# Matrix to convert xyz to uvw for each timestamp (w is not considered)
# Matrix to convert xyz to uvw for each timestamp (but w is not considered)
m = np.array([[np.sin(hourangle), np.cos(hourangle), np.zeros(len(hourangle))],
[-np.sin(self.target.dec)*np.cos(hourangle),
np.sin(self.target.dec)*np.sin(hourangle),
@ -386,10 +657,20 @@ class Observation(object):
return bl_uv_up
def get_uv_array(self):
"""Returns a 2-D array with the (u,v) points from this observation.
Note that complex conjugates are not provided, thus one should compute
V(-u, -v), V(u, v) if the full uv coverage is desired.
def get_uv_array(self) -> np.ndarray:
"""Returns the (u, v) values for each baseline and each timestamp for which the source
is visible.
The difference with `get_uv_baseline` is that `get_uv_array` only returns the (u,v)
values, dropping the information of baselines and times to which these values belong to.
Returns a (N, 2)-dimensional numpy.ndarray containing all N (u,v) points resulting for
each timestamp and baseline. The (u,v) values are given in lambda units.
Note that complex conjugate values are not provided.
Exceptions
- It may raise the exception SourceNotVisible if no baselines can observe the source
at all during the observation.
"""
if self._uv_array is not None:
return self._uv_array
@ -408,15 +689,23 @@ class Observation(object):
return self._uv_array
def synthesized_beam(self):
def synthesized_beam(self) -> dict:
"""Estimates the resulting synthesized beam of the observations based on
the expected uv coverage.
This is just an estimation made by a ellipse fitting to the uv coverage,
the expected (u,v) coverage.
This is just an estimation made by a ellipse fitting to the (u, v) coverage,
from which we obtain the resolution on the two axes following
https://science.nrao.edu/facilities/vlba/docs/manuals/oss/ang-res
theta_HPBW (mas) \sim 2063 x lambda(cm)/b_max^km
returns a dict with the following keys: bmaj, bmin, pa.
Note that different robust weighting maps during imaging would provide slightly
different synthesized beams. The provided value here does not assumed any weighting
in the data. A natural weighting would thus likely provide a slightly larger beam,
while an uniform weighting would provide a slightly smaller beam.
Returns a dict with the following keys: 'bmaj', 'bmin', 'pa' (major axis, minor axis,
and position angle).
The three values are astropy.units.Quantity objects with units of angles.
"""
if self._synth_beam is not None:
return self._synth_beam
@ -454,21 +743,26 @@ class Observation(object):
def print_obs_times(self, date_format='%d %B %Y'):
"""Given an observation, it returns the time range (starttime-endtime) in a smart
way. If the observation lasts for less than one day it omits the end date:
20 Jan 1971 10:00-20:00UT
It also adds the GST range after that.
"""Returns the time range (starttime-to-endtime) of the observation in a smart way.
If the observation lasts for less than one day it omits the end date:
20 January 1971 10:00-20:00 UTC
GST: 05:00-15:00
If the observation ends the day after, then it returns:
20 January 1971 10:00-20:00 UTC (+1d)
GST: 05:00-15:00
If the observation is longer, then it returns
20 January 1971 10:00 to 24 January 1971 20:00 UTC
GST: 05:00-15:00
Input:
- obs : observation.Observation
It must already have set the .times part with an array of astropy.Time times.
- date_format : str [optional]
Format for the date part (only the date part) of the string to represent
the time range.
- date_format : str [OPTIONAL]
Format for the date part (only the date part, not the times) of
the string to represent the time range.
Output:
- printed_time : str
- strtime : str
A string showing the time-range of the observation.
"""
gsttext = "{:02n}:{:02.2n}-{:02n}:{:02.2n}".format((self.gstimes[0].hour*60) // 60,
(self.gstimes[0].hour*60) % 60,

486
vlbiplanobs/stations.py

@ -1,3 +1,6 @@
# -*- coding: utf-8 -*-
# Licensed under GPLv3+ - see LICENSE
from __future__ import annotations
import configparser
from importlib import resources
import numpy as np
@ -5,40 +8,63 @@ from astropy import units as u
from astropy import coordinates as coord
from astropy.io import ascii
from astropy.time import Time
from astroplan import Observer
from astroplan import Observer, FixedTarget
class Station(object):
"""Module that defines the `Station` and `Stations` objects, which represent a station (antenna)
or a network composed of antennas.
"""
__all__ = ['Station', 'SelectedStation', 'Stations']
def __init__(self, name, codename, network, location, freqs_sefds, min_elevation=20*u.deg,
fullname=None, all_networks=None, country='', diameter='', real_time=False):
"""Initializes a station. The given name must be the name of the station that
observes, with the typical 2-letter format used in the EVN (with exceptions).
class Station(object):
"""Defines an astronomical station (antenna).
A station is defined by some names and codenames, coordinates, and its sensitivity for the
radio bands (wavelengths) it can observe.
Apart of the metadata related to the station, it allows to compute the altitude/azimuth, elevation,
or simply when a source is visible from the station for a given time range.
"""
def __init__(self, name: str, codename: str, network: str, location: coord.EarthLocation,
freqs_sefds: dict, min_elevation=20*u.deg, fullname: str = None,
all_networks: str = None, country: str = '', diameter: str = '', real_time: bool = False):
"""Initializes a station.
Inputs
- name : str
Name of the observer (the station that is going to observe).
If it contains undercores (_), they will be converted to blank spaces.
- codename : str
A code for the name of the station. It can be the same as name.
A short code (accronym) for the name of the station. It is meant to follow the standard approach
from the EVN: an (often) two-letter code unique for each station.
- network : str
Name of the network to which the station belongs.
- location : EarthLocation
Position of the observer on Earth.
Name of the network to which the station belongs (e.g. EVN).
- location : astropy.coordinates.EarthLocation
Position of the station on Earth in (x,y,z) gecentric coordinates.
- freqs_sefds : dict
Dictionary with all frequencies the station can observe, and as values
the SEFD at each frequency.
- min_elevation : Quantity
Minimum elevation that the station can observe a source. If no units
provided, degrees are assumed. By default it 20 degrees.
Dictionary with all frequencies the station can observe as keys of the dictionary, and the
values representing the system equivalent flux density (SEFD; in Jansky units)
at each frequency.
Although the key format is in principle free, we recommend to use the syntax 'XXcm' (str type).
This will be then consistent with the default station catalog.
- min_elevation : Quantity [OPTIONAL]
Minimum elevation that the station can reach to observe a source. If no units (astropy.units)
provided, degrees are assumed. By default it 20 degrees. It does not support an azimuth-dependent
elevation limits.
- fullname : str [OPTIONAL]
Full name of the station. If not given, `name` is assumed.
Full name of the station. If not given, same as `name` is assumed.
It can be used to expand the full name if an abbreviation is typically used for the name.
For example, name: VLA, fullname: Karl G. Jansky Very Large Array.
- all_networks : str [OPTIONAL]
Networks where the station can participate (free style).
Networks where the station can participate (free style string).
- country : str [OPTIONAL]
Country where the station is placed.
Country where the station is located.
- diameter : str [OPTIONAL]
Diameter of the station (free format).
- real_time : bool [OPTIONAL, False by default]
If the station can participate in real-time observations (e.g. e-EVN).
Diameter of the station (free format string). We recommend a syntax of e.g. '30 m' for normal
single-dish antennas, and in case of interferometers it can have a form like '25 x 20 m',
meaning that the station is composed of 25 antennas of 20 m each.
- real_time : bool [OPTIONAL]
If the station can participate in real-time observations (e.g. e-EVN), False by default.
"""
self.observer = Observer(name=name.replace('_', ' '), location=location)
self._codename = codename
@ -65,201 +91,370 @@ class Station(object):
@property
def name(self):
def name(self) -> str:
"""Name of the station.
"""
return self.observer.name
@property
def codename(self):
def codename(self) -> str:
"""Codename of the station (typically a two-letter accronym).
"""
return self._codename
@property
def fullname(self):
def fullname(self) -> str:
"""Full name of the station. If not specified, it can be the same as 'name'.
"""
return self._fullname
@property
def network(self):
def network(self) -> str:
"""Name of the network to which the station belongs.
"""
return self._network
@property
def all_networks(self):
def all_networks(self) -> str:
"""Name of all networks to which the station belongs.
If not specified it can be the same as 'network'.
"""
return self._all_networks
@property
def country(self):
def country(self) -> str:
"""Country where this station is located.
It can be an empty string if not specified.
"""
return self._country
@property
def diameter(self):
def diameter(self) -> str:
"""String representing the diameter of the station, and/or how many antennas compose
the station in case of connected-interferometers.
"""
return self._diameter
@property
def real_time(self):
def real_time(self) -> bool:
"""If the station can participate in real-time observations (e.g. e-EVN).
"""
return self._real_time
@property
def location(self):
"""Location of the station in EarthLocation type.
def location(self) -> coord.EarthLocation:
"""Location of the station as an astropy.coordinates.EarthLocation object.
"""
return self.observer.location
@property
def bands(self):
"""Observing bands the station can observe.
Returns a dict_keys object with all bands in a string format as introduced in the freqs_sefd
attribute when the Station was created.
"""
return self._freqs_sefds.keys()
@property
def sefds(self):
"""Returns a dictionary with the SEFDs for each of the frequencies
the station can observe (given as keys).
def sefds(self) -> dict:
"""Returns a dictionary with the system equivalent flux density (SEFDs) for each
of the frequencies the station can observe (given as keys).
"""
return self._freqs_sefds
@property
def min_elevation(self):
def min_elevation(self) -> u.Quantity:
"""Minimum elevation the station can observe a source.
Returns an astropy.units.Quantity (i.e. number with units).
"""
return self._min_elev
def elevation(self, obs_times, target):
"""Returns the elevation of the source as seen by the Station during obs_times.
def elevation(self, obs_times: Time, target: FixedTarget) -> coord.angles.Latitude:
"""Returns the elevation of the target source as seen by the Station during obs_times.
Inputs
- obs_times : astropy.time.Time
Time to compute the elevation of the source (either single time or a list of times).
Time to compute the elevation of the source
(either a single timestamp or an array of times).
- target : astroplan.FixedTarget
Target to observe.
Output
- elevations : ndarray
Elevation of the source at the given obs_times
- elevations : astropy.coordinates.angles.Latitute
Elevation of the source at the given obs_times.
"""
# source_altaz = source_coord.transform_to(coord.AltAz(obstime=obs_times,
# location=self.location))
# return source_altaz.alt
return self.observer.altaz(obs_times, target).alt
def altaz(self, obs_times, target):
"""Returns the altaz coordinates of the target for the given observing times.
def altaz(self, obs_times: Time, target: FixedTarget) -> coord.sky_coordinate.SkyCoord:
"""Returns the altaz coordinates of the target source for the given observing times.
Inputs
- obs_times : astropy.time.Time
Time to compute the elevation of the source
(either a single timestamp or an array of times).
- target : astroplan.FixedTarget
Target coordinates to observe.
Output
- altaz : astropy.coordinates.sky_coordinate.SkyCoord
Altitude and azimuth of the source for each given time.
"""
return self.observer.altaz(obs_times, target)
def is_visible(self, obs_times, target):
"""Return if the source is visible for this station at the given times.
def is_visible(self, obs_times: Time, target: FixedTarget) -> tuple:
"""Returns when the target source is visible for this station at the given times.
Inputs
- obs_times : astropy.time.Time
Time to compute the elevation of the source
(either a single timestamp or an array of times).
- target : astroplan.FixedTarget
Target coordinates to observe.
Output
- visible : tuple
Tuple containing the indexes of obs_times when the target source is visible
from the station. Therefore obs_times[visible] would return only those times.
"""
elevations = self.elevation(obs_times, target)
return np.where(elevations >= self.min_elevation)
def has_band(self, band):
def has_band(self, band: str) -> bool:
"""Returns if the Station can observed the given band `the_band`.
Inputs
- band : str
A string representing an observing band, following the same syntax as used
when the station was initialized and the bands where defined in the keys of
the freqs_sefds attribute.
Output
- bool whenever the station has the given observing band.
"""
return band in self.bands
def sefd(self, band):
"""Returns the SEFD of the Station at the given band.
def sefd(self, band: str) -> float:
"""Returns the system equivalent flux density (SEFD) of the Station at the given band,
in Jansky (Jy) units.
Input
- band : str
A string representing an observing band, following the same syntax as used
when the station was initialized and the bands where defined in the keys of
the freqs_sefds attribute.
Output
- SEFD : float
The SEFD at the given band, in Jy units.
Exception
- It may raise KeyError if the given band is not available for this station.
"""
return self._freqs_sefds[band]
def __str__(self):
return f"<{self.codename}>"
def __repr__(self):
return f"<stations.Station: {self.codename}>"
return f"<Station: {self.codename}>"
class SelectedStation(Station):
def __init__(self, name, codename, network, location, freqs_sefds, min_elevation=20*u.deg,
fullname=None, all_networks=None, country='', diameter='', real_time=False, selected=True):
"""Extends the Station class with an additional attribute: `selected` (bool).
This allows the sub-seting of lists of stations by considering all of them to observe
or just disabling some of them.
"""
def __init__(self, name: str, codename: str, network: str, location: coord.EarthLocation,
freqs_sefds: dict, min_elevation=20*u.deg, fullname: str = None,
all_networks: str = None, country: str = '', diameter: str = '',
real_time: bool = False, selected: bool = True):
"""Initializes a SelectedStation.
This class extends Station by adding one additional attribute: `selected` (bool).
Inputs
- name : str
Name of the observer (the station that is going to observe).
If it contains undercores (_), they will be converted to blank spaces.
- codename : str
A short code (accronym) for the name of the station. It is meant to follow the standard approach
from the EVN: an (often) two-letter code unique for each station.
- network : str
Name of the network to which the station belongs (e.g. EVN).
- location : astropy.coordinates.EarthLocation
Position of the station on Earth in (x,y,z) gecentric coordinates.
- freqs_sefds : dict
Dictionary with all frequencies the station can observe as keys of the dictionary, and the
values representing the system equivalent flux density (SEFD; in Jansky units)
at each frequency.
Although the key format is in principle free, we recommend to use the syntax 'XXcm' (str type).
This will be then consistent with the default station catalog.
- min_elevation : Quantity [OPTIONAL]
Minimum elevation that the station can reach to observe a source. If no units (astropy.units)
provided, degrees are assumed. By default it 20 degrees. It does not support an azimuth-dependent
elevation limits.
- fullname : str [OPTIONAL]
Full name of the station. If not given, same as `name` is assumed.
It can be used to expand the full name if an abbreviation