spacecruft
/
lofarimaging
1
0
Fork 0
Code Issues Releases Activity

Merge 3c177be359 into e183746b8f

pull/3/merge
Mattia Mancini 2020-06-17 08:14:12 +00:00 committed by GitHub
parent e183746b8f 3c177be359
commit 1313681c22
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 273 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
lofarimaging.ipynb
View File

@ -170,8 +170,7 @@
"output_type": "stream",
"text": [
"Searching for available files in ./test\n",
"0: ./test/20200525_084652_mode_5_xst_sb100.dat\n",
"1: ./test/20170720_095816_mode_3_xst_sb297.dat\n"
"0: ./test/20170720_095816_mode_3_xst_sb297.dat\n"
]
}
],
@ -198,7 +197,7 @@
],
"source": [
"# Select a file\n",
"xst_filename = files[1]\n",
"xst_filename = files[0]\n",
"\n",
"print(\"File selected:\", xst_filename)"
]

148
lofarimaging/lofarimaging.py
View File

@ -1,15 +1,18 @@
"""Functions for working with LOFAR single station data"""
import logging
from typing import Dict, List
import numpy as np
from numpy.linalg import norm, lstsq
import numexpr as ne
import numba
from astropy.coordinates import SkyCoord, SkyOffsetFrame, CartesianRepresentation
import numba
import numexpr as ne
import numpy as np
import scipy.optimize
from astropy.coordinates import SkyCoord, SkyOffsetFrame, CartesianRepresentation
from numpy.linalg import norm, lstsq
__all__ = ["nearfield_imager", "sky_imager", "ground_imager", "skycoord_to_lmn", "calibrate", "simulate_sky_source",
"subtract_sources"]
"subtract_sources", "compute_calibrated_model", "compute_pointing_matrix", "estimate_model_visibilities",
"estimate_sources_flux", "self_cal"]
__version__ = "1.5.0"
SPEED_OF_LIGHT = 299792458.0
@ -170,13 +173,144 @@ def calibrate(vis, modelvis, maxiter=30, amplitudeonly=True):
return residual, gains
def compute_calibrated_model(vis, model_vis, maxiter=30):
"""
Calibrate the gains phases and apply to the model
Args:
vis: visibility matrix, shape [n_st, n_st]
model_vis: model visibility matrix, shape [n_st, n_st]
maxiter: max iterations (default 30)
Returns:
calibrated: model visibilities, shape [n_st, n_st]
gains: the antenna gains, shape [n_st]
residual: amplitude difference between model and actual visibilities, float
"""
n_ant = vis.shape[0]
gain_phase = np.zeros((n_ant), dtype=complex)
def gains_difference(gain_phase, vis, model_vis):
gains = np.exp(1.j * gain_phase)
gains_mat = np.diag(gains)
gains_mat_star = np.diag(np.conj(gains))
predicted = gains_mat_star @ model_vis @ gains_mat
residual = np.sum(np.abs(vis - predicted))
return residual
result = scipy.optimize.minimize(gains_difference, gain_phase, args=(vis, model_vis), options={'maxiter': maxiter})
gain_phase = result.x
gains = np.exp(1.j * gain_phase)
gains_mat = np.diag(gains)
calibrated = np.conj(gains_mat) @ model_vis @ gains_mat
return calibrated, gains, gains_difference(gain_phase, vis, model_vis)
def compute_pointing_matrix(sources_positions, baselines, frequency):
"""
Compute the pointing matrix used to compute the visibilities from
the sources position
Args:
sources_positions: matrix of the lmn sources position, shape [n_dir, 3]
baselines: baseline matrix, shape [n_st, n_st]
frequency: frequency
Returns:
pointing_matrix: pointing matrix, shape [n_st, n_st, n_dir]
"""
n_baselines = baselines.shape[0]
n_sources = sources_positions.shape[0]
pointing_matrix = np.zeros(shape=(n_baselines, n_baselines, n_sources), dtype=complex)
for q in range(n_sources):
pointing_matrix[:, :, q] = simulate_sky_source(sources_positions[q, :], baselines, frequency)
return pointing_matrix
def estimate_sources_flux(visibilities, pointing_matrix):
"""
Estimate the flux of the sources given the pointing matrix and the measured visibilities
Args:
visibilities: visibilities, shape [n_st, n_st]
pointing_matrix: pointing matrix, shape [n_st, n_st, n_dir]
Returns:
fluxes: amplitude of the source at the given direction, shape [n_st, n_st, n_dir]
"""
def residual_flux(signal, vis, point_matrix):
residual = vis - point_matrix @ signal
residual = np.sum(np.abs(residual))
return residual
n_antennas = visibilities.shape[0]
n_sources = pointing_matrix.shape[2]
lin_visibilities = visibilities.reshape((n_antennas * n_antennas))
lin_pointing_matrix = pointing_matrix.reshape((n_antennas * n_antennas, n_sources))
fluxes, *_ = np.linalg.lstsq(lin_pointing_matrix, lin_visibilities)
result = scipy.optimize.minimize(residual_flux, fluxes, args=(visibilities, pointing_matrix))
fluxes = result.x
return fluxes
def estimate_model_visibilities(sources_positions, visibilities, baselines, frequency):
"""
Estimate the visibilities of the models given the sources position and the measured visibilities
Args:
sources_positions: lmn coords of the sources, shape [n_dir, 3]
visibilities: visibilities, shape [n_st, n_st]
baselines: baselines matrix in m, shape [n_st, n_st, 3]
frequency: frequency
Returns:
model: the model visibilities, shape [n_st, n_st]
"""
pointing_matrix = compute_pointing_matrix(sources_positions, baselines, frequency)
fluxes = estimate_sources_flux(visibilities, pointing_matrix)
model = pointing_matrix @ fluxes
return model
def self_cal(visibilities, expected_sources, baselines, frequency, max_iterations=10, tolerance=1.e-4):
"""
Compute the gain phase comparing the model and the actual visibilities returning
the model with the gain phase correction applied.
Args:
visibilities: visibilities, shape [n_st, n_st]
expected_sources: lmn coords of the sources, shape [n_dir, 3]
baselines: baselines matrix in m, shape [n_st, n_st, 3]
frequency: frequency
max_iterations: maximum number of iterations to perform
Returns:
model: the model visibilities with the applied gains, shape [n_st, n_st]
gains: the antenna gains, shape [n_st]
"""
model = estimate_model_visibilities(expected_sources, visibilities, baselines, frequency)
model, gains, residual = compute_calibrated_model(visibilities, model, maxiter=50)
for i in range(1, max_iterations):
model, gains, next_residual = compute_calibrated_model(visibilities, model, maxiter=50)
if next_residual != 0 and abs(1 - (residual / next_residual)) >= tolerance:
residual = next_residual
else:
break
return model, gains
def simulate_sky_source(lmn_coord: np.array, baselines: np.array, freq: float):
"""
Simulate visibilities for a sky source
Args:
lmn_coord (np.array): l, m, n coordinate
baselines (np.array): baseline distances in metres, shape (n_ant, n_ant)
baselines (np.array): baseline distances in metres, shape (n_ant, n_ant, 3)
freq (float): Frequency in Hz
"""
return np.exp(2j * np.pi * freq * baselines.dot(np.array(lmn_coord)) / SPEED_OF_LIGHT)

88
lofarimaging/singlestationutil.py
View File

@ -1,36 +1,32 @@
"""Functions for working with LOFAR single station data"""
import os
import datetime
import os
from typing import List, Dict, Tuple, Union
import numpy as np
from packaging import version
import tqdm
import h5py
import matplotlib.pyplot as plt
import matplotlib.animation
from matplotlib.ticker import FormatStrFormatter
from matplotlib import cm
from matplotlib.figure import Figure
from matplotlib.colors import ListedColormap, Normalize
from matplotlib.patches import Circle
import matplotlib.axes as maxes
from mpl_toolkits.axes_grid1 import make_axes_locatable
from astropy.coordinates import SkyCoord, GCRS, EarthLocation, AltAz, get_sun, get_moon
import astropy.units as u
from astropy.time import Time
import lofargeotiff
from lofarantpos.db import LofarAntennaDatabase
import h5py
import lofarantpos
import lofargeotiff
import matplotlib.animation
import matplotlib.axes as maxes
import matplotlib.pyplot as plt
import numpy as np
import tqdm
from astropy.coordinates import SkyCoord, GCRS, EarthLocation, AltAz, get_sun, get_moon
from astropy.time import Time
from lofarantpos.db import LofarAntennaDatabase
from matplotlib import cm
from matplotlib.colors import ListedColormap, Normalize
from matplotlib.figure import Figure
from matplotlib.patches import Circle
from matplotlib.ticker import FormatStrFormatter
from mpl_toolkits.axes_grid1 import make_axes_locatable
from packaging import version
from .maputil import get_map, make_leaflet_map
from .lofarimaging import nearfield_imager, sky_imager, skycoord_to_lmn, subtract_sources
from .hdf5util import write_hdf5
from .lofarimaging import nearfield_imager, sky_imager, skycoord_to_lmn, subtract_sources
from .maputil import get_map, make_leaflet_map
__all__ = ["sb_from_freq", "freq_from_sb", "find_caltable", "read_caltable",
"rcus_in_station", "read_acm_cube", "get_station_pqr", "get_station_type",
@ -331,7 +327,7 @@ def get_station_pqr(station_name: str, rcu_mode: Union[str, int], db):
'intl': GENERIC_INT_201512, 'remote': GENERIC_REMOTE_201512, 'core': GENERIC_CORE_201512
}
selected_dipoles = selected_dipole_config[station_type] + \
np.arange(len(selected_dipole_config[station_type])) * 16
np.arange(len(selected_dipole_config[station_type])) * 16
station_pqr = db.hba_dipole_pqr(full_station_name)[selected_dipoles]
else:
raise RuntimeError("Station name did not contain LBA or HBA, could not load antenna positions")
@ -340,7 +336,7 @@ def get_station_pqr(station_name: str, rcu_mode: Union[str, int], db):
def make_ground_plot(image: np.ndarray, background_map: np.ndarray, extent: List[float], title: str = "Ground plot",
subtitle: str = "", opacity: float = 0.6, fig: Figure = None, draw_contours: bool = True, **kwargs) \
subtitle: str = "", opacity: float = 0.6, fig: Figure = None, draw_contours: bool = True, **kwargs) \
-> Tuple[Figure, np.ndarray]:
"""
Make a ground plot of an array with data
@ -585,7 +581,6 @@ def make_xst_plots(xst_data: np.ndarray,
opacity: Opacity for map overlay. Defaults to 0.6.
hdf5_filename: Filename where hdf5 results can be written. Defaults to outputpath + '/results.h5'
outputpath: Directory where results can be saved. Defaults to 'results'
subtract: List of sources to subtract. Defaults to None
Returns:
@ -662,14 +657,14 @@ def make_xst_plots(xst_data: np.ndarray,
marked_bodies = {
'Cas A': SkyCoord(ra=350.85 * u.deg, dec=58.815 * u.deg),
'Cyg A': SkyCoord(ra=299.86815191 * u.deg, dec=40.73391574 * u.deg),
'Per A': SkyCoord(ra=49.95066567*u.deg, dec=41.51169838 * u.deg),
'Her A': SkyCoord(ra=252.78343333*u.deg, dec=4.99303056*u.deg),
'Cen A': SkyCoord(ra=201.36506288*u.deg, dec=-43.01911267*u.deg),
'Vir A': SkyCoord(ra=187.70593076*u.deg, dec=12.39112329*u.deg),
'3C295': SkyCoord(ra=212.83527917*u.deg, dec=52.20264444*u.deg),
'Per A': SkyCoord(ra=49.95066567 * u.deg, dec=41.51169838 * u.deg),
'Her A': SkyCoord(ra=252.78343333 * u.deg, dec=4.99303056 * u.deg),
'Cen A': SkyCoord(ra=201.36506288 * u.deg, dec=-43.01911267 * u.deg),
'Vir A': SkyCoord(ra=187.70593076 * u.deg, dec=12.39112329 * u.deg),
'3C295': SkyCoord(ra=212.83527917 * u.deg, dec=52.20264444 * u.deg),
'Moon': get_moon(obstime_astropy, location=station_earthlocation).transform_to(GCRS),
'Sun': get_sun(obstime_astropy),
'3C196': SkyCoord(ra=123.40023371*u.deg, dec=48.21739888*u.deg)
'3C196': SkyCoord(ra=123.40023371 * u.deg, dec=48.21739888 * u.deg)
}
marked_bodies_lmn = {}
@ -752,7 +747,7 @@ def make_xst_plots(xst_data: np.ndarray,
"pixels_per_metre": pixels_per_metre}
tags.update(calibration_info)
lon_min, lon_max, lat_min, lat_max = extent_lonlat
lofargeotiff.write_geotiff(ground_img[::-1,:], os.path.join(outputpath, f"{fname}_nearfield_calibrated.tiff"),
lofargeotiff.write_geotiff(ground_img[::-1, :], os.path.join(outputpath, f"{fname}_nearfield_calibrated.tiff"),
(lon_min, lat_max), (lon_max, lat_min), as_pqr=False,
stationname=station_name, obsdate=obstime, tags=tags)
@ -769,7 +764,7 @@ def make_sky_movie(moviefilename: str, h5file: h5py.File, obsnums: List[str], vm
"""
Make movie of a list of observations
"""
fig = plt.figure(figsize=(10,10))
fig = plt.figure(figsize=(10, 10))
for obsnum in tqdm.tqdm(obsnums):
obs_h5 = h5file[obsnum]
skydata_h5 = obs_h5["sky_img"]
@ -787,7 +782,28 @@ def make_sky_movie(moviefilename: str, h5file: h5py.File, obsnums: List[str], vm
# Thanks to Maaijke Mevius for making this animation work!
ims = fig.get_children()[1:]
ims = [ims[i:i+2] for i in range(0, len(ims), 2)]
ims = [ims[i:i + 2] for i in range(0, len(ims), 2)]
ani = matplotlib.animation.ArtistAnimation(fig, ims, interval=30, blit=False, repeat_delay=1000)
writer = matplotlib.animation.writers['ffmpeg'](fps=5, bitrate=800)
ani.save(moviefilename, writer=writer, dpi=fig.dpi)
def compute_baselines(station_name, rcu_mode):
# Setup the database
db = LofarAntennaDatabase()
station_pqr = get_station_pqr(station_name, rcu_mode, db)
station_name = get_full_station_name(station_name, rcu_mode)
# Rotate station_pqr to a north-oriented xyz frame, where y points North, in a plane through the station.
rotation = db.rotation_from_north(station_name)
pqr_to_xyz = np.array([[np.cos(-rotation), -np.sin(-rotation), 0],
[np.sin(-rotation), np.cos(-rotation), 0],
[0, 0, 1]])
station_xyz = (pqr_to_xyz @ station_pqr.T).T
baselines = station_xyz[:, np.newaxis, :] - station_xyz[np.newaxis, :, :]
return baselines

1
requirements.txt
View File

@ -13,3 +13,4 @@ Pillow
opcua
h5py
tqdm
scipy

77
test/test_single_station.py 100644
View File

@ -0,0 +1,77 @@
import numpy
from numpy.testing import assert_almost_equal
from lofarimaging.lofarimaging import compute_calibrated_model, \
compute_pointing_matrix, simulate_sky_source, \
estimate_sources_flux, \
estimate_model_visibilities, \
self_cal
def test_compute_calibrated_model():
# ----TEST DATA
vis = numpy.diag(numpy.ones((10)))
model_vis = numpy.diag(numpy.ones(10))
# ----TEST
calibrated, gains, residual = compute_calibrated_model(vis, model_vis)
assert_almost_equal(calibrated, vis)
assert residual < 1.e-5
def test_simulate_sky_source():
# ----TEST DATA
source_position = [0, 0, 0]
baselines = numpy.zeros((2, 2, 3))
baselines[0, 1, :] = [1, 0, 0]
baselines[1, 0, :] = [1, 0, 0]
# ----TEST
visibilities = simulate_sky_source(source_position, baselines, 1)
expected = numpy.ones((2, 2))
assert_almost_equal(visibilities, expected)
def test_compute_pointing_matrix():
# ----TEST DATA
source_position = numpy.zeros((1, 3))
baselines = numpy.zeros((2, 2, 3))
baselines[0, 1, :] = [1, 0, 0]
baselines[1, 0, :] = [1, 0, 0]
expected = numpy.ones((2, 2, 1))
# ----TEST
pointing_matrix = compute_pointing_matrix(source_position, baselines, 1)
assert_almost_equal(pointing_matrix, expected)
def test_estimate_sources_flux():
# ----TEST DATA
source_position = numpy.zeros((1, 3))
pointing_matrix = numpy.ones((2, 2, 1))
visibilities = numpy.ones((2, 2))
# ----TEST
source_fluxes = estimate_sources_flux(visibilities, pointing_matrix)
expected = numpy.ones((1))
assert_almost_equal(source_fluxes, expected)
def test_estimate_model_visibilities():
# ----TEST DATA
source_position = numpy.zeros((1, 3))
visibilities = numpy.ones((2, 2))
baselines = numpy.zeros((2, 2, 3))
baselines[0, 1, :] = [1, 0, 0]
baselines[1, 0, :] = [1, 0, 0]
# ----TEST
model_visibilities = estimate_model_visibilities(source_position, visibilities, baselines, 1)
assert_almost_equal(visibilities, model_visibilities)
def test_self_cal_zero_residual():
# ----TEST DATA
source_position = numpy.zeros((1, 3))
visibilities = numpy.ones((2, 2))
baselines = numpy.zeros((2, 2, 3))
baselines[0, 1, :] = [1, 0, 0]
baselines[1, 0, :] = [1, 0, 0]
# ----TEST
calibrated_model, _ = self_cal(visibilities, source_position, baselines, 1)
assert_almost_equal(calibrated_model, visibilities)