Welcome to ska-fastimgproto’s documentation!

Version 0

Contents:

Jupyter Notebook Examples

Kernel Functions

Various analytic functions used for generating ‘kernels’ (see below). Some are very simple (e.g. pillbox, triangle), which makes them useful when testing other code routines. Some (e.g. gaussian-sinc, prolate spheroidal) have desirable properties relating to Fourier transforms.

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 6, 6

import fastimgproto.gridder.conv_funcs as conv_funcs

triangle3 = conv_funcs.Triangle(half_base_width=3.0)
pillbox = conv_funcs.Pillbox(half_base_width=2.5)
sinc = conv_funcs.Sinc(3)
gauss = conv_funcs.Gaussian(trunc=5.)
g_sinc = conv_funcs.GaussianSinc(trunc=5.)
narrow_g_sinc = conv_funcs.GaussianSinc(trunc=3.)

plot_radius = 10.
x=np.linspace(-plot_radius,plot_radius,501)

# %matplotlib notebook
fig = plt.figure()
ax = fig.gca()
# ax.plot(x, pillbox(x), color='r', label='pillbox')
# ax.plot(x, triangle3(x), color='g', label='triangle')

ax.set_xticks(np.arange(-5,6,1))
ax.set_yticks(np.arange(-1,1.1,0.1))
ax.plot(x, sinc(x), color='b', label='sinc')
ax.plot(x, gauss(x), color='m', label='Gaussian')
ax.plot(x, g_sinc(x), color='k', ls=':', label='Gaussian-sinc')
ax.plot(x, narrow_g_sinc(x), color='k', ls='--')
ax.grid()
ax.set_xlim(-5,5)
ax.set_ylim(-0.5,1.1)
ax.legend(loc='best')

<matplotlib.legend.Legend at 0x7fe6979106d0>

Convolution kernel generation

We use kernel functions to generate a small pixel-grid (‘kernel of convolution’):

from fastimgproto.gridder.kernel_generation import Kernel
gs_kernel = Kernel(kernel_func=narrow_g_sinc, support=3)

fig = plt.figure()
ax = fig.gca()
plt.imshow(gs_kernel.array)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fe6977c1350>

fig = plt.figure()
ax = fig.gca()
cross_section = gs_kernel.array[gs_kernel.centre_idx]
x_pixels_idx = np.arange(len(gs_kernel.array)) - gs_kernel.centre_idx
ax.grid()
ax.plot(x_pixels_idx, cross_section, ls='--', alpha=0.5, color='b', lw=3,
       label='cross-section')
ax.plot(x, narrow_g_sinc(x)*np.max(cross_section), ls=':', color='r', lw=3,
       label='reference curve')
ax.set_xlim(-4,4)
ax.legend(loc='upper right')

<matplotlib.legend.Legend at 0x7fe6976f3a50>

Anti-aliasing properties of kernel functions

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 6, 6

import fastimgproto.gridder.conv_funcs as conv_funcs
from fastimgproto.gridder.kernel_generation import Kernel
narrow_g_sinc = conv_funcs.GaussianSinc(trunc=3.)
oversampling=1
support = 32
gs_kernel = Kernel(kernel_func=narrow_g_sinc, support=support,oversampling=oversampling)
cross_section = gs_kernel.array[gs_kernel.centre_idx]

plt.plot(cross_section)

[<matplotlib.lines.Line2D at 0x7f610cde2ad0>]

shifted = np.fft.ifftshift(cross_section)

shifted[0]

0.41675735397798913

shifted[-1], shifted[1]

(0.15770968729136908, 0.15770968729136908)

amp = np.abs(np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(cross_section))))
fig = plt.figure()
ax = fig.gca()
# plot_section = np.log10(amp[:oversampling])
ax.plot(np.arange(len(amp))-support, amp)
ax.grid()

C_grid = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(gs_kernel.array)))

plt.imshow(np.abs(C_grid))

<matplotlib.image.AxesImage at 0x7f610cc2eb10>

plt.plot(np.abs(C_grid)[support])

[<matplotlib.lines.Line2D at 0x7f610cb70bd0>]

Example visibility-simulation and imaging

import numpy as np
import logging
import os

import astropy.units as u
import astropy.constants as const
import scipy.stats

import fastimgproto.imager as imager
import fastimgproto.visibility as visibility

from astropy.coordinates import Angle, SkyCoord, AltAz, EarthLocation
from astropy.time import Time

from fastimgproto.gridder.conv_funcs import GaussianSinc
from fastimgproto.skymodel.helpers import SkyRegion, SkySource
from fastimgproto.sourcefind.image import SourceFindImage
from fastimgproto.telescope.readymade import Meerkat

# from tqdm import tqdm_notebook as Tqdm
from tqdm import tqdm as Tqdm

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 10, 10

telescope = Meerkat()
print("Telescope with {} antennae == {} baselines".format(
    len(telescope.ant_local_xyz), len(telescope.baseline_local_xyz)))
print("Centre: {!r}, {!r}".format(telescope.lon, telescope.lat))

Telescope with 64 antennae == 2016 baselines
Centre: <Longitude 21.443270030936162 deg>, <Latitude -30.712485463619643 deg>

pointing_centre = SkyCoord(0 * u.deg, -30 * u.deg)
obs_central_frequency = 3. * u.GHz
wavelength = const.c / obs_central_frequency
transit_time = telescope.next_transit(pointing_centre.ra,
                                      start_time=Time('2017-01-01'))

/home/docs/checkouts/readthedocs.org/user_builds/ska-fastimaging-python/envs/latest/lib/python2.7/site-packages/astropy-2.0.1-py2.7-linux-x86_64.egg/astropy/time/core.py:749: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif longitude == 'greenwich':

Let’s check that the transit-time calculation is approximately correct. We chose the target SkyCoord to have the same celestial latitude as Meerkat’s geographical latitude, so it should transit near zenith:

altaz = pointing_centre.transform_to(
    AltAz(obstime=transit_time,
         location=telescope.centre))
altaz.alt.deg

89.17591321744462

nstep=10
obs_times = transit_time + np.linspace(-1, 1, nstep) * u.hr
print("Generating UVW-baselines for {} timesteps".format(nstep))
uvw_m = telescope.uvw_tracking_skycoord(pointing_centre, obs_times)
# From here on we use UVW as multiples of wavelength, lambda:
uvw_lambda = (uvw_m / wavelength).to(u.dimensionless_unscaled).value



# Additional source to North-East of pointing centre
extra_src_position = SkyCoord(ra=pointing_centre.ra + 0.01 * u.deg,
                              dec=pointing_centre.dec + 0.01 * u.deg, )

steady_sources = [
    SkySource(pointing_centre, flux=1 * u.Jy),
    SkySource(extra_src_position, flux=0.4 * u.Jy),
]

# Simulate incoming data; includes transient sources, noise:
print("Simulating visibilities")
data_vis = visibility.visibilities_for_source_list(
    pointing_centre,
    source_list = steady_sources,
    uvw = uvw_lambda)

vis_noise_level = 0.1 * u.Jy
data_vis = visibility.add_gaussian_noise(vis_noise_level, data_vis)
vis_weights =  np.minimum(1.0, scipy.stats.lognorm.rvs(s=0.25, size=len(data_vis)))

Generating UVW-baselines for 10 timesteps
Simulating visibilities

print("Simulated {} visibilities".format(len(data_vis)))

Simulated 20160 visibilities

image_size=1024 * u.pixel
cell_size=1 * u.arcsecond

kernel_support = 3
kernel_func = GaussianSinc(trunc=kernel_support)
with Tqdm() as progress_bar:
    image, beam = imager.image_visibilities(vis=data_vis,
                                            vis_weights=vis_weights,
                                            uvw_lambda=uvw_lambda,
                                            image_size=image_size,
                                            cell_size=cell_size,
                                            kernel_func=kernel_func,
                                            kernel_support=kernel_support,
                                            kernel_exact=True,
                                            kernel_oversampling=None,
                                            progress_bar=progress_bar
                                            )
image = np.real(image)

Gridding visibilities100%|██████████| 20160/20160 [00:05<00:00, 4013.87it/s]

Finally, let’s run our rudimentary sourcefinder on the image, and plot the results:

from fastimgproto.sourcefind.image import SourceFindImage

detection_n_sigma=20
analysis_n_sigma=15
sfimage = SourceFindImage(data=np.real(image),
                          detection_n_sigma=detection_n_sigma,
                          analysis_n_sigma=analysis_n_sigma,
                          )

sfimage.islands

[IslandParams(parent=<fastimgproto.sourcefind.image.SourceFindImage object at 0x7fa04ba69790>, label_idx=1, sign=1, extremum_val=0.40170229680359693, extremum_x_idx=481, extremum_y_idx=476, xbar=481.33041435161402, ybar=477.02546415144235),
 IslandParams(parent=<fastimgproto.sourcefind.image.SourceFindImage object at 0x7fa04ba69790>, label_idx=3, sign=1, extremum_val=1.0023925008556798, extremum_x_idx=512, extremum_y_idx=512, xbar=511.79100217442266, ybar=511.73170038622021)]

src = sfimage.islands[0]

fig, ax1 = plt.subplots(1,1)
ax1.imshow(image)
axlims = 400, 600
ax1.set_xlim(axlims)
ax1.set_ylim(axlims)
for src in sfimage.islands:
    ax1.axvline(src.xbar, ls=':')
    ax1.axhline(src.ybar, ls=':')

Sourcefinding example

from astropy.io import fits
from astropy.stats import sigma_clip
import numpy as np
from scipy import ndimage

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 10, 10

Load the data. A mask can be applied if necessary - this may be useful e.g. for excluding the region around a bright source, to avoid false detections due to sidelobes.

# fits_path = '../casapy-simulation-scripts/simulation_output/vla.image.fits'
# hdu0 = fits.open(fits_path)[0]
# img_data = hdu0.data.squeeze()
# imgdata = np.ma.MaskedArray(imgdata, mask=np.zeros_like(imgdata))
# imgdata.mask[900:1100,900:1100] = True
# imgdata.mask.any()

Alternatively, we can simulate a rudimentary image by adding a Gaussian model to background noise, although note that the noise will be uncorrelated, unlike a radio-synthesis image:

import fastimgproto.fixtures.image as fixture
img_shape = (128,192)
img_data = fixture.uncorrelated_gaussian_noise_background(img_shape,sigma=0.1)
srcs = []
srcs.append(fixture.gaussian_point_source(x_centre=32.5, y_centre=32.66, amplitude=0.5))
srcs.append(fixture.gaussian_point_source(x_centre=64.12, y_centre=48.88, amplitude=1.0))
srcs.append(fixture.gaussian_point_source(x_centre=128.43, y_centre=94.5, amplitude=1.5))
for s in srcs:
    img_data += fixture.evaluate_model_on_pixel_grid(img_data.shape, s)

imgmax = np.max(img_data)
plt.imshow(img_data,vmax=imgmax*0.5)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fafc84cb510>

from fastimgproto.sourcefind.image import SourceFindImage
sfimage = SourceFindImage(data = img_data, detection_n_sigma=5, analysis_n_sigma=3)

Background level and RMS are crudely estimated via median and sigma-clipped std. dev., respectively:

sfimage.bg_level

0.00059247585779820266

sfimage.rms_est

0.099449361884591136

We use two thresholds when identifying our source ‘islands’ (connected pixel regions). The high threshold is our detection level, and should be set high enough to avoid false detections due to noise spikes. The lower threshold expands each island, such that it is more likely to contain enough pixels to reasonably fit a Gaussian profile (otherwise the island may consist of only a single pixel over the detection threshold).

Note that this thresholding approach may result in multi-peaked regions (e.g. two distinct but adjacent sources) being assigned to a single island / label. This can be tackled with ‘deblending’ algorithms if desired, but is not covered in this notebook.

The thresholded data is then run through scipy.ndimage.label which numbers the connected regions:

plt.imshow(sfimage.label_map,cmap='Paired')

<matplotlib.image.AxesImage at 0x7fafc84595d0>

Plotting all data which is merely above the analysis threshold shows the importance of a usefully high detection threshold - there are many noise spikes above the analysis threshold.

plt.imshow(sfimage.data > sfimage.analysis_n_sigma*sfimage.rms_est, cmap='binary')

<matplotlib.image.AxesImage at 0x7fafc835cc50>

Each island is then analysed for peak value, barycentre, etc (and in may be model-fitted in future):

island = sfimage.islands[1]
island

IslandParams(parent=<fastimgproto.sourcefind.image.SourceFindImage object at 0x7fafc8253810>, label_idx=13, sign=1, extremum_val=0.94112170230174208, extremum_x_idx=64, extremum_y_idx=49, xbar=64.067026526873718, ybar=48.92987056569585)

plt.imshow(island.data)
plt.xlim(island.extremum_x_idx-10,island.extremum_x_idx+10,)
plt.ylim(island.extremum_y_idx-10,island.extremum_y_idx+10,)
plt.scatter(island.xbar,island.ybar, marker='*', s=200, c='y',)

<matplotlib.collections.PathCollection at 0x7fafc69e4d50>

print("Bright source model:")
print(srcs[-1])
print()
print("Island barycentres:")
for i in sfimage.islands:
    print(i.xbar, i.ybar)

Bright source model:
Model: Gaussian2D
Inputs: (u'x', u'y')
Outputs: (u'z',)
Model set size: 1
Parameters:
    amplitude x_mean y_mean x_stddev y_stddev theta
    --------- ------ ------ -------- -------- -----
          1.5 128.43   94.5      1.5      1.2   1.0
()
Island barycentres:
(128.50568016955827, 94.631554654595774)
(64.067026526873718, 48.92987056569585)

API Reference

fastimgproto.coords module

fastimgproto.coords.time_of_next_transit(observer_longitude, target_ra, start_time, tolerance=<Quantity 0.02 arcsec>)

Calculate time, t when the local-hour-angle of target will be 0.

(Such that start_time <= t < start_time + 24*u.hr.)

NB: This gives no guarantees about visibility / horizons, as we are ignoring observer-declination / target-declination here. Just provides a useful way of setting up an observation at a good hour-angle, for a given time-period.

Parameters:

• observer_longitude (astropy.coordinates.Longitude) – Longitude of position on Earth
• target_ra (astropy.coordinates.Longitude) – Right ascension of sky-target.
• start_time (astropy.time.Time) – Time to start from.
• tolerance (astropy.units.Quantity) – If target’s LHA is within tolerance of zero at start_time, simply return start_time. Otherwise look for the next transit. Dimensions of angular width (arc).

Returns:

Approximate time of the next transit

(good to around 0.02 arcsec)

Return type:

astropy.time.Time

fastimgproto.coords.xyz_to_uvw_rotation_matrix(hour_angle, declination)

Generates rotation matrix from XYZ -> UVW co-ordinates, for given HA, Dec.

For an explanation and derivation, see Notebook 4.1, Sections 4.1.1d and 4.1.2 of https://github.com/griffinfoster/fundamentals_of_interferometry (“The (u,v,w) Space”)

The definition of XYZ used here is that:

• X axis points towards the crossing of the celestial equator and meridian at hour-angle \(H=0h\)
• Y axis points towards hour-angle \(H=-6h\) (Due East from telescope, in the local-XYZ / local hour-angle case.)
• Z axis points towards the North Celestial Pole (NCP).

The definition of UVW used here is that:

• the u-axis lies in the celestial equatorial plane, and points toward the hour angle \(H_0-6h\) (where \(H_0\) is the hour-angle of the source) (i.e. East of source).
• the v-axis lies in the plane of the great circle with hour angle \(H_0\), and points toward the declination and points toward the declination \(\frac{\pi}{2}-\delta_0\) (North of source).
• the w-axis points in the direction of the source.

Can be used to generate a rotation matrix for either local XYZ co-ordinates, or for Global XYZ (also known as ECEF or ITRF) co-ordinates. In the case of local-XYZ co-ords the hour-angle \(H_0\) supplied should be Local-hour-angle (relative to transit, \(LHA = LST - RA\)). In the case of Global XYZ co-ords it should be Greenwich hour angle (\(GHA = GST - RA\), where GST is Greenwich Sidereal Time).

NB LST = GST + Longitude, equivalently LHA = GHA + Longitude (with the standard convention of East-is-positive)

Parameters:

• hour_angle (astropy.units.Quantity) – Source hour-angle. Dimensions of angular width (arc).
• dec (astropy.units.Quantity) – Source declination. Dimensions of angular width (arc).

Returns:

Rotation matrix for converting from XYZ to UVW. [Dtype numpy.float_, shape (3, 3)]

Return type:

numpy.ndarray

fastimgproto.coords.z_rotation_matrix(rotation_angle)

Rotation matrix for a passive transformation counter-clockwise about Z-axis

(i.e. rotation of co-ordinate system)

Parameters:rotation_angle (astropy.units.Quantity) – Rotation-angle. Dimensions of angular width (arc). Returns:Rotation matrix Return type:numpy.ndarray

fastimgproto.gridder

Modules

fastimgproto.gridder.conv_funcs module

Convolution functions.

Used for generating the kernel used in convolutional gridding.

We actually paramerize the functions at initialization and return a simple callable with one parameter, the distance in pixels.

This allows us to pass the convolution routine the minimum of extra parameters.

class fastimgproto.gridder.conv_funcs.ConvFuncBase(trunc)

Implements truncation (via __call__), numpy array reshaping.

Always returns 0 outside truncation radius, i.e.:

if np.fabs(x) > trunc:
    conv_func(x)==0 # True

Parameters:trunc – truncation radius.

f(radius)

The convolution function to be evaluated and truncated

class fastimgproto.gridder.conv_funcs.Gaussian(trunc, w=1.0)

Gaussian function (with truncation).

evaluates the function:

exp(-(x/w)**2)

(Using the notation of Taylor 1998, p143, where x = u/delta_u and alpha==2. Default value of w=1).

Parameters:

• trunc – truncation radius.
• w (float) – Width normalization of the Gaussian. Default = 1.0

class fastimgproto.gridder.conv_funcs.GaussianSinc(trunc, w1=2.52, w2=1.55)

Gaussian times sinc function (with truncation).

evaluates the function:

exp(-(x/w1)**2) * sinc(x/w2)

(Using the notation of Taylor 1998, p143, where x = u/delta_u and alpha==2. Default values for w1,w2 are chosen according to recommendation therein).

Parameters:

• trunc – truncation radius.
• w1 (float) – Width normalization of the Gaussian. Default = 2.52
• w2 (float) – Width normalization of the sinc. Default = 1.55

class fastimgproto.gridder.conv_funcs.Pillbox(half_base_width)

Valued 1.0 from origin to half_base_width, zero thereafter.

AKA ‘TopHat’ function.

Symmetric about the origin.

Makes a terrible anti-aliasing function. But, because it’s so simple, it’s easy to verify and therefore a useful tool in verifying convolution codes.

Parameters:half_base_width (float) – Half-base width pillbox.

class fastimgproto.gridder.conv_funcs.Sinc(trunc)

Sinc function (with truncation).

class fastimgproto.gridder.conv_funcs.Triangle(half_base_width)

Linearly declines from 1.0 at origin to 0.0 at half_base_width, zero thereafter. “

Symmetric about the origin.

Makes a terrible anti-aliasing function. But, because it’s so simple, it’s easy to verify and therefore a useful tool in verifying convolution codes.

Parameters:half_base_width (float) – Half-base width of the triangle.

fastimgproto.gridder.gridder module

Convolutional gridding of visibilities.

fastimgproto.gridder.gridder.calculate_oversampled_kernel_indices(subpixel_coord, oversampling)

Find the nearest oversampled gridpoint for given sub-pixel offset.

Effectively we are mapping the range [-0.5, 0.5] to the integer range [-oversampling//2, ..., oversampling//2].

Inputs will be between -0.5 and 0.5 inclusive. This is an issue, because inputs at the extreme (e.g. 0.5) might round UP, taking them outside the desired integer output range. We simply correct this edge-case by replacing outlier values before returning.

Parameters:

• subpixel_coord (numpy.ndarray) – Array of ‘fractional’ co-ords, that is the subpixel offsets from nearest pixel on the regular grid. dtype: float, shape: (n_vis, 2).
• oversampling (int) – How many oversampled pixels to one regular pixel.

Returns:

Corresponding oversampled pixel indexes. These are in oversampled pixel widths from the kernel centre pixel, to a maximum of half a regular pixel, so they have integer values ranging from -oversampling/2 to oversampling/2. [Dtype: int, shape: (n_vis, 2)].

Return type:

numpy.ndarray

fastimgproto.gridder.gridder.convolve_to_grid(kernel_func, support, image_size, uv, vis, vis_weights, exact=True, oversampling=0, raise_bounds=True, progress_bar=None)

Grid visibilities using convolutional gridding.

Returns the un-normalized weighted visibilities; the weights-renormalization factor can be calculated by summing the sample grid.

If exact == True then exact gridding is used, i.e. the kernel is recalculated for each visibility, with precise sub-pixel offset according to that visibility’s UV co-ordinates. Otherwise, instead of recalculating the kernel for each sub-pixel location, we pre-generate an oversampled kernel ahead of time - so e.g. for an oversampling of 5, the kernel is pre-generated at 0.2 pixel-width offsets. We then pick the pre-generated kernel corresponding to the sub-pixel offset nearest to that of the visibility.

Kernel pre-generation results in improved performance, particularly with large numbers of visibilities and complex kernel functions, at the cost of introducing minor aliasing effects due to the ‘step-like’ nature of the oversampled kernel. This in turn can be minimised (at the cost of longer start-up times and larger memory usage) by pre-generating kernels with a larger oversampling ratio, to give finer interpolation.

Parameters:

• kernel_func (callable) – Callable object, (e.g. conv_funcs.Pillbox,) that returns a convolution co-efficient for a given distance in pixel-widths.
• support (int) – Defines the ‘radius’ of the bounding box within which convolution takes place. Box width in pixels = 2*support+1. (The central pixel is the one nearest to the UV co-ordinates.) (This is sometimes known as the ‘half-support’)
• image_size (int) – Width of the image in pixels. NB we assume the pixel [image_size//2,image_size//2] corresponds to the origin in UV-space.
• uv (numpy.ndarray) – UV-coordinates of visibilities. 2d array of float_, shape: (n_vis, 2). assumed ordering is u-then-v, i.e. u, v = uv[idx]
• vis (numpy.ndarray) – Complex visibilities. 1d array, shape: (n_vis,).
• vis_weights (numpy.ndarray) – Visibility weights. 1d array, shape: (n_vis,).
• exact (bool) – Calculate exact kernel-values for every UV-sample.
• oversampling (int) – Controls kernel-generation if exact==False. Larger values give a finer-sampled set of pre-cached kernels.
• raise_bounds (bool) – Raise an exception if any of the UV samples lie outside (or too close to the edge) of the grid.
• progress_bar (tqdm.tqdm) – [Optional] progressbar to update.

Returns:

(vis_grid, sampling_grid)

Tuple of ndarrays representing the gridded visibilities and the sampling weights. These are 2d arrays of same dtype as vis, shape (image_size, image_size). Note numpy style index-order, i.e. access like vis_grid[v,u].

Return type:

tuple

fastimgproto.gridder.gridder.populate_kernel_cache(kernel_func, support, oversampling)

Generate a cache of normalised kernels at oversampled-pixel offsets.

We need kernels for offsets of up to oversampling//2 oversampling-pixels in any direction, in steps of one oversampling-pixel (i.e. steps of width 1/oversampling in the original co-ordinate system).

Parameters:

• kernel_func (callable) – Callable object, (e.g. conv_funcs.Pillbox,) that returns a convolution co-efficient for a given distance in pixel-widths.
• support (int) – See kernel generation routine.
• oversampling (int) – Oversampling ratio. cache_size = ((oversampling // 2 * 2) + 1)**2

Returns:

Dictionary mapping oversampling-pixel offsets to normalised kernels.

Return type:

dict

fastimgproto.gridder.kernel_generation module

Generation of convolution kernel arrays, with optional sub-pixel origin offset and and oversampling.

class fastimgproto.gridder.kernel_generation.Kernel(kernel_func, support, offset=(0.0, 0.0), oversampling=None, pad=False, normalize=True)

Generates a 2D array representing a sampled kernel function.

Parameters:

• kernel_func (callable) – Callable object, (e.g. conv_funcs.Pillbox,) that returns a convolution co-efficient for a given distance in pixel-widths.
• support (int) – Defines the ‘radius’ of the bounding box within which convolution takes place. Box width in pixels = 2*support+1. (The central pixel is the one nearest to the UV co-ordinates.) For a kernel_func with truncation radius trunc, the support should be set to ceil(trunc+0.5) to ensure that the kernel function is fully supported for all valid subpixel offsets.
• offset (tuple) – 2-vector subpixel offset from the sampling position of the central pixel to the origin of the kernel function. Ordering is (x_offset,y_offset). Should have values such that fabs(offset) <= 0.5 otherwise the nearest integer grid-point would be different!
• oversampling (int) – Oversampling ratio, how many kernel pixels to each UV-grid pixel. Defaults to 1 if not given or oversampling=None is passed.
• pad (bool) – Whether to pad the array by an extra pixel-width. This is used when generating an oversampled kernel that will be used for interpolation.

array

numpy.ndarray – The sampled kernel function.

centre_idx

int – Index of the central pixel

kernel_func, support, offset, oversampling

See params.

fastimgproto.telescope

Modules

fastimgproto.telescope.base module

class fastimgproto.telescope.base.Telescope(centre, ant_labels=None, ant_itrf_xyz=None, ant_local_xyz=None)

Data-structure for holding together array layout and central position.

Also takes care of baseline calculation and antenna/baseline label-tracking.

Instantiate using the named constructors, e.g.:

meerkat = Telescope.from_itrf(
    *parse_itrf_ascii_to_xyz_and_labels(meerkat_itrf_filepath)
)

centre

astropy.coordinates.EarthLocation – Centre of the array

ant_labels

list – Antennae labels

ant_itrf_xyz

numpy.ndarray – Antennae xyz_positions in ITRF frame. [dtype: np.float_, shape (n_antennae,3,).

ant_local_xyz

numpy.ndarray – Antennae xyz_positions in local-XYZ frame. [dtype: np.float_, shape (n_antennae,3,).

baseline_local_xyz

numpy.ndarray – Baseline xyz-vectors in local-XYZ frame. [dtype: np.float_, shape (n_baselines,3,).

baseline_labels

list – Baseline labels. (See also generate_baselines_and_labels())

static from_itrf(ant_itrf_xyz, ant_labels)

Instantiate telescope model from ITRF co-ordinates of antennae

(aka ECEF-coords, i.e. Earth-Centered-Earth-Fixed reference frame.)

Takes care of calculating central Latitude and Longitude from the mean antenna position, and converts the antenna positions into local-XYZ frame.

Parameters:

• ant_itrf_xyz (numpy.ndarray) – Array co-ordinatates in the ITRF frame
• ant_labels (list[str]) – Antennae labels

Returns:

A telescope class with the given array co-ords.

Return type:

Telescope

lha(ra, time)

Calculate the local hour angle of a target-RA at a given time

Parameters:

• ra (astropy.coordinates.Longitude) – Right ascension
• time (astropy.time.Time) – Timestamp

Returns:

Local Hour Angle of target-RA.

Return type:

astropy.coordinates.Longitude

lst(time)

Calculate the local sidereal time at the telescope

Parameters:time (astropy.time.Time) – Global timestamp Returns:Local sidereal time expressed as an angle-Quantity. Return type:astropy.coordinates.Longitude

next_transit(target_ra, start_time)

Wrapper around time_of_next_transit()

See time_of_next_transit() for details.

Parameters:

• target_ra (astropy.coordinates.Longitude) – Right ascension of sky-target.
• start_time (astropy.time.Time) – Time to start from.

Returns:

Approximate time of the next transit

Return type:

astropy.time.Time

uvw_at_local_hour_angle(local_hour_angle, dec)

Transform baselines to UVW for source at given local-hour-angle and Dec.

Parameters:

• local_hour_angle (astropy.units.Quantity) – Local hour angle. LHA=0 is transit, LHA=-6h is rising, LHA=+6h is setting. [Dimensions of angular width / arc, e.g. u.deg, u.rad, or u.hourangle]
• dec (astropy.units.Quantity) – Declination. [Dimensions of angular width / arc]

Returns:

UVW-array, with units of metres. [Shape: (3, n_baselines), dtype: numpy.float_].

Return type:

astropy.units.Quantity

uvw_tracking_skycoord(pointing_centre, obs_times, progress_bar=None)

Calculate the UVW-array towards pointing centre for all obs_times.

This function calculates an ndarray of uvw-coords corresponding to an instantaneous observation at each of obs_times, concatenating the results to produce an ndarray of shape (n_baselines * len(obs_times), 3)

Parameters:

• pointing_centre (astropy.coordinates.SkyCoord) – Pointing centre co-ords for UVW calculation.
• obs_times (list) – List of astropy.time.Time, the instants of observation.
• progress_bar (tqdm.tqdm) – [Optional] progressbar to update.

Returns:

UVW-array, with units of metres.

Return type:

astropy.units.Quantity

fastimgproto.telescope.base.generate_baselines_and_labels(antenna_positions, antenna_labels)

Given an array of antenna positions, generate an array of baseline-vectors.

Baseline-vectors are generated as follows:

• It is assumed that the antenna-positions supplied are already in the desired order, with matching labels in the antenna_labels list. (Typically something like ['ANT1','ANT2','ANT3'].

• Pairings are then generated according to position in the supplied list, as follows:

  pairings = list(itertools.combinations(range(len(antenna_positions)), 2))
  

• Baseline labels are concatenated with a comma, e.g. pairing of 'ANT1' and 'ANT2' has label 'ANT1,ANT2'.

• Baseline vectors represent translation from first antenna to second antenna, therefore:

  baseline_1_2 = ant_2_pos - ant_1_pos
  

Parameters:

• antenna_positions (astropy.units.Quantity) – Antennae xyz_positions in local-XYZ frame. [dtype: np.float_, shape (n_antennae,3,).
• antenna_labels (list) – Antennae labels

Returns:

(baseline_vecs, baseline_labels). Both of length \(n_{baselines} = {n_{antennae} \choose 2}\), where baseline_labels is simply a list of strings, and baseline_vecs is a (numpy.ndarray) [dtype: np.float_, shape (n_baselines,3,).

Return type:

tuple

fastimgproto.telescope.base.hstack_table_columns_as_ndarray(cols)

When passed a subset of astropy Table columns, turn them into a 2-d ndarray.

(Note, columns should have same dtype)

Example

>>> import numpy as np
>>> import astropy.table
>>> x = np.arange(5)
>>> y = x**2
>>> t = astropy.table.Table(data=(x,y), names=('x','y'))
>>> hstack_table_columns_as_ndarray(t.columns[:2])
array([[ 0,  0],
       [ 1,  1],
       [ 2,  4],
       [ 3,  9],
       [ 4, 16]])

Parameters:cols (list) – List (or iterable) of astropy.table.TableColumns. Typically a table or a slice of a table’s columns. Returns:H-stacked columns-array. Return type:numpy.ndarray

fastimgproto.telescope.base.parse_itrf_ascii_to_xyz_and_labels(tabledata)

Read ITRF data in ascii format and return (antenna positions, labels).

ASCII table should be simple whitespace-delimited, column order: x y z diameter label mount

Parameters:tabledata – (str, file-like, or list). File-object or path-to-file. Returns:Tuple of (ant_itrf_xyz, ant_labels), where ant_itrf_xyz is a numpy.ndarray of shape (n,3), dtype np.float_ and ant_labels is a numpy.ndarray of shape (n,), dtype np.str_. Return type:tuple

Indices and tables

• Index
• Module Index
• Search Page