from IPython.display import Image
from IPython.terminal import embed
from IPython.core.display import display, HTML, Latex
from astropy.coordinates import EarthLocation, Distance, SkyCoord
from astropy import units as u
nch = EarthLocation(lat="37.8732", lon="-122.2573", height=123.1*u.m)
%matplotlib notebook
'''load laundry list of modules from <libradiomods.py>'''
import libradiomods
modules = libradiomods.get_modules()
for module in modules:
try:
exec(module)
except ImportError:
print("Error ", module)
continue
'''load functions from <libradiofns.py>'''
from libradiofns import *
First all the of raw data will be imported. This is often called a "measurement set". There is often weather information, seeing conditions, and space weather notes included. These will be appended in an ad-hoc way to reconstruct a more robust assessment of the observations.
crab_nebula = np.load('crab_6hrs_3-31-2020.npz')['arr_0']
orion_nebula = np.load('orion_3hrs_4-3-2020.npz')['arr_0']
virgo = np.load('virgoA_9.5hrs_4-2-2020.npz')['arr_0']
cassiopeia = np.load('casA_24hrs_4-5-2020.npz')['arr_0']
cygnus = np.load('cygA_17hrs_4-5-2020.npz')['arr_0']
targets = ['crab_nebula', 'orion_nebula', 'virgo', 'cassiopeia', 'cygnus']
observations = { targets[i]:globals()[targets[i]] for i in range(len(targets)) }
The FITs format is probably the most common and well supported ways to package observations. It allows for arbitrarily large images in both size (pixels) and in terms of per pixel encoding. Most computers are not able to show more than 24 bit colors with R, G, and B channels each having $2^8 = 256$ bits. Images in industry and the science can require "stretching" before being rendered. This is basically the same as using
$$ \log{P_\nu} $$to see a power spectrum except with other options. Having an image calibrated in terms of location can also be important if it is going to be composed with an other observation. Plus if you want to show some scale or astronomical coordinates then it will be useful to have the image wrapped with the important parameters needed to transform from one coordinate system to another.
The specifications of FITs files are available here: {https://fits.gsfc.nasa.gov/fits_standard.html}
# use astropy to get object coordinates
nch = EarthLocation(lat="37.8732", lon="-122.2573", height=123.1*u.m)
crab_nebula_wcs = SkyCoord.from_name("crab nebula")
orion_nebula_wcs = SkyCoord.from_name("orion nebula")
virgo_wcs = SkyCoord.from_name("virgo A")
cassiopeia_wcs = SkyCoord.from_name("cassiopeia A")
cygnus_wcs = SkyCoord.from_name("cygnus A")
wcsnames = [ "".join([targets[i], "_wcs"]) for i in range(len(targets)) ]
skycoordinates = { targets[i]:globals()[wcsnames[i]] for i in range(len(targets)) }
The methods above will get coordinates on the celestial sphere but will not provide any information about distance. These objects have similar radio fluxes, with varying spectral energy distributions, however two of them are extragalactic. The best way to calibrate the telescope is with a point source because it effectively normalizes the range of beams and sensitivities. If there is an established measurement of the absolute flux from a source then that value can be used as a fixed variable in solving for the gain of a telescope. The more often this is done the more checks you will have to adjust the pipeline to get accurate measures of error, accuracy, and precision.*
*A good reference for this is J. Baars, “History of Flux-Density Calibration in Radio Astronomy,” The Radio Science Bulletin, vol. March, no. 348, pp. 47-66, (2014). These sources can vary in time so just as the coordinates of the telescope are updated so should the calibration factors.
# import standard catalog distances manually
crab_nebula_z = (6500*u.lyr).to(u.pc); crab_nebula_z = Distance(crab_nebula_z)
orion_nebula_z = (1344*u.lyr).to(u.pc); orion_nebula_z = Distance(orion_nebula_z)
virgo_z = (5.35e7*u.lyr).to(u.pc); virgo_z = Distance(virgo_z)
cassiopeia_z = (11090*u.lyr).to(u.pc); cassiopeia_z = Distance(cassiopeia_z)
cygnus_z = (6e8*u.lyr).to(u.pc); cygnus_z = Distance(cygnus_z)
targets_z = [ "".join([targets[i], "_z"]) for i in range(len(targets))]
z_lookup = { targets[i]:globals()[f'{targets_z[i]}'] for i in range(len(targets)) }
One of the important measures is the beam size of the telescope. For single dish radio telescopes the "full width half maximum" is typically used as a proxy for this. The Rayleigh criterion is used
$$ \theta = \frac{\lambda}{2R} $$to approximate this, with $\lambda$ being the wavelength of the band of interest and $R$ being the radius of the antenna/aperture. For two element interferometers $2R \rightarrow B$ and is the typically much larger length of the baseline. In this special case the result is more like the spacing of a fringe pattern than a beam. The resolution is increased, however the sensitivity is not. For larger arrays it is the same although the result is a synthesized beam usually shown as an elliptical stereographic projection. There is an astropy-related package for radio beams which can be used.
# constants (priors)
beam_fwhm = 3*u.degree # the larger beam of the antenae sets the field of view (FOV)
baseline = (20)*u.m # estimate the baseline
xband = 10.5e9*u.hertz # set the observation frequency/wavelength
xband_wavelength = (xband).to(u.cm, equivalencies=u.spectral()) # get the right units, see astropy units
print(f'The X-band at {xband:1.2e} is equivalently {xband_wavelength:1.3f}.')
# get field of view
fwhm_to_sigma = 1. / (8 * np.log(2))**0.5 # this is a Gaussian-like factor to shape the beam
beam_sigma = beam_fwhm * fwhm_to_sigma # make the shaped beam
omega_B = 2 * np.pi * beam_sigma**2 # get the beam solid angle
# get resolution
interfres = (xband_wavelength/baseline).to(u.dimensionless_unscaled) # calculate interf resolution (radians)
interfres_degrees = (interfres*u.radian).to(u.degree) # convert to coordinate-like units
interfres_arcmin = (interfres*u.radian).to(u.arcmin) # convert to coordinate-like units
interfres_arcsec = (interfres*u.radian).to(u.arcsec) # convert to coordinate-like units
'''instantiate a Beam obejct'''
interfbeam = Beam(interfres_degrees)
print(f'The interf dishes have a beam size approximated by a FWHM of {beam_fwhm} and solid angle {omega_B:2.1f} [sr]')
print(f'The interf baseline corresponds to a synthetic beam size of {interfres_degrees:3.2e}')
print(f'The interf baseline corresponds to a synthetic beam size of {interfres_arcmin:3.2e}')
print(f'The interf baseline corresponds to a synthetic beam size of {interfres_arcsec:3.2e}')
The beam object now has functionality. It can be projected to sources, observations can be simulated, and archival files can be re-calibrated. The FITs format is good for packaging, however there are better data structures when using python like pandas dataframes. This provides a more flexible interface with other modules geared towards quantitative needs.
# print out the beam pattern area; what is the best way to size up extended radio sources?
# if the beam area is less than the object area how many beams fill the area?
for target in targets:
print(f'For {target} @ {z_lookup[target]:1.2e} the interf beam covers {interfbeam.beam_projected_area(distance=z_lookup[target]):1.3e}')
This step should include checking for other potential sources around the field of view of the beam. The reason is any additional flux entering into the main beam or the side lobes could reduce the fidelity of the observation so filtering or at least accounting for this is a good idea.
One of the most crucial parts of the FITs file is the header. This is were all the metadata is stored. The coordinates of the source, the name of the source, the telescope and so on. There is a standardized format for this detailed at the link given above and in the included PDF. A blank template for radio observations is also included and will be filled in below. astropy treats the header like a python dictionary. The FITs keywords are the keys and their values are the values in the dictionary. The header allows for interchange between formats hence "Flexible Image Transport System".
# import template file and get the header
template = fits.open('template.fits')
template_header = template[0].header # pick the first HDU in the HDU list
template_data = template[0].data
print(repr(template_header)) # repr() grabs the __repr__ of an object, the print representation
To get the observations organized correctly, the data is first loaded into the primary header/data unit (HDU). For these observations the FITs files will be simply the primary HDU. There are other more complicated types of FITs files including "cubes" with a third axis. This is common with spectroscopic observations over range of frequencies.
# all the observation data is used to generate HDU objects
hdus = dict()
for target in targets:
hdus[target] = fits.PrimaryHDU(observations[target])
hdus[target].header = template_header
# the primary HDU is supposed to be in an HDU list so the last step is repeated
hdu_lists = dict()
for target in targets:
hdu_lists[target] = fits.HDUList(hdus[target])
With the data ready all that is needed is to modify the headers. Generating headers is possible with astropy and makes sense to do as a script-generated object while observing. However, for building it up by hand like this it will be a little easier to use a template and change the values to match the observations. Editing FITs headers can be tedious and something like fv FITs viewer is probably the way to go.*
This template example shows three axes: NAXIS1 is a pixel image x-axis, NAXIS2 is a pixel image y-axis, and NAXIS3 is a pixel image z-axis. The full file (this is actually an integrated version) is like a movie with the frames being the elements of the z-axis. The data from the interf observations should be similarly summarized. The units of the x-axis will be time series measurements, the y-axis a voltage reading from the measurement, and the z-axis can be removed to reduce clutter. The value of CDELT1 and CDELT2 should also be changed to reflect the increments along each axis. For the x-axis it could be the time between measurements and the y-axis would be the range divided by the bit depth of the sampler.
Before making the headers it would be good to grab all the weather information that is missing. There are tutorials on how to get this type of information with an API. There's an example implemented in libradiofns but you will have to get your own API Key if you want to use it.
*https://www.visualcrossing.com/resources/documentation/weather-api/weather-api-documentation/
After sending in a query the database will return a good amount of information. Astropy wants the pressure, temperature, and relative_humidity to make corrections. Noting cloud cover, wind speed, and precipitation might be good too so they can be included as notes.
# start by updating all the SkyCoord objects to get derivied values
# at the same time make a dictionary for the weather information
weather_info = dict()
for target in targets:
observation_time = Time(hdus[target].data[1][0], format='unix')
skycoordinates[target].obstime = observation_time
skycoordinates[target].location = nch
weather_info[target] = get_weather(observation_time.value, observation_time.value)
pressure = float(weather_info[target]['Sea Level Pressure'][0])
skycoordinates[target].pressure = (pressure*u.mbar).to(u.hPa)
temperature = float(weather_info[target]['Temperature'])
skycoordinates[target].temperature = f2C(temperature)*u.deg_C
skycoordinates[target].relative_humidity = float(weather_info[target]['Relative Humidity'])
skycoordinates[target].obswl = xband_wavelength.to(u.micron)
# the SkyCoord objects are all set and now provide alt/az info with atmospheric corrections
# these could be compared with the uncorrected values to get a more robust error
for target in targets:
print(f'{target} is located in the constellation {skycoordinates[target].get_constellation()}')
alt = skycoordinates[target].altaz.alt; az = skycoordinates[target].altaz.az;
print(f'{target} had the following corrected alt/az coordinates {alt:3.1f} {az:3.1f}\n')
# the weather data now stored as pandas data frames which can be displayed
for target in targets:
display(weather_info[target])
The weather information can be written into a second HDU (with its own header) of the FITs files and the rest can go into the more-or-less typical locations. The Beam object from above will be the same for all of them although the projection would vary the geometry. If there is a way to convert the data into an image then new FITs files can be made with WCS axes to allow accurate positions when plotting.*
The HDU lists need to be updated and the weather information can be appended before writing each HDU list to separate FITs files. Thankfully there are a few built-in methods that make going from pandas to FITs not too tortuous.
# the number of time series measurements is the length of the second array
# in the .data extension of each HDU the voltage readings are uncalibrated
# so the value can be set to the number of mesurements but the increment
# should be set to uncalibrated; comments are changed by passing tuples
for target in targets:
header = hdu_lists[target][0].header
header['OBJECT'] = (target, 'Point source calibrator for interf')
header['OBSERVER'] = 'spacewaves'
header['TELESCOP'] = 'interf'
header['INTRUME'] = 'HP Multimeter'
header['BTYPE'] = ('Uncalibrated Intensity', 'These should be set once calibrated')
header['BSCALE'] = ('Uncalibrated')
header['BUNIT'] = ('Uncalibated', 'Needs to be calibrated into Jy/beam')
header['BMAJ'] = (interfbeam.major.value, 'Major axis of beam ellipse [degrees]')
header['BMIN'] = (interfbeam.minor.value, 'Minor axis of beam ellipse [degrees]')
header['BPA'] = (interfbeam.beam_projected_area(z_lookup[target]).value, 'Beam projected area [pc^2]')
header['RESTFRQ'] = (xband.value/1e9, 'GHz +- 20 MHz')
header['DATE-OBS'] = (ugradio.timing.utc(hdus[target].data[1][0], '%Y-%m-%d'), 'UTC')
header['NAXIS'] = (3, 'Number of axes in the data')
header['NAXIS1'] = (hdus[target].data[1].shape[0], 'Unix time of time series measurements')
header['NAXIS2'] = (hdus[target].data[0].shape[0], 'Voltage readings (time series)')
header['CTYPE1'] = ('s', 'Seconds with unix epoch offset')
header['CTYPE2'] = ('voltage', 'Uncalibrated')
header['CUNIT1'] = ('s', 'Seconds with unix epoch offset')
header['CUNIT2'] = ('voltage', 'Needs to be calibrated into Jy/beam')
header['CDELT1'] = (.2, 'HP mm @ dt = 0.1 -> sampling slower/varied')
header['CDELT2'] = (1, 'Needs to be calibrated into Jy/beam')
header['CRVAL1'] = (hdus[target].data[1][0], 'Unix time at start')
header['CRVAL2'] = (0, 'voltage should be zero without source')
header['CRPIX1'] = (0, 'Reference pixel coordinates')
header['CRPIX2'] = (0, 'Reference pixel coordinates')
header['OBSRA'] = (skycoordinates[target].ra.value, 'Right ascension [degrees]')
header['OBSDEC'] = (skycoordinates[target].dec.value, 'Declination [degrees]')
header['OBSALT'] = (skycoordinates[target].altaz.alt.value, 'Altitude [degrees]')
header['OBSAZ'] = (skycoordinates[target].altaz.az.value, 'Azimuth [degrees]')
v_radial_correction = skycoordinates[targets[0]].radial_velocity_correction().to(u.km/u.s).value
header['VRAD'] = (v_radial_correction, 'km/s add this value to observered v_radial')
header['OBSLAT'] = nch.lat.value
header['OBSLON'] = nch.lon.value
header['NAXIS3'] = (1, 'Could add proper motion of object')
header['CTYPE3'] = ('parsecs', f'Distance to {target} in pc')
header['CRVAL3'] = (z_lookup[target].value, 'Distance to {target} in pc')
header['CRPIX3'] = (0, 'Reference pixel coordinates')
header['CONSTEL'] = (skycoordinates[target].get_constellation(), 'Constellation')
header['CDELT3'] = (1, 'Could add proper motion of object')
header['CUNIT3'] = ('parsecs', f'Distance to {target} in pc')
header['DATE'] = (ugradio.timing.local_time(), 'UTC')
header['ORIGIN'] = ('RADIOLAB', 'generated with astropy')
start = hdus[target].data[1][0]; stop = hdus[target].data[1][-1]
header['ELAPTIME'] = f'total observation time was {(stop-start)/3600} hours'
header['HISTORY'] = 'weather information added with table in second HDU'
table = astropy.table.Table.from_pandas(weather_info[target])
hdu_lists[target].append(fits.table_to_hdu(table))
hdu_lists[target].writeto(f'{target}.fits')
crab_nebula_fits = fits.open('crab_nebula.fits')
print(repr(crab_nebula_fits[0].header))