Data handling#

All image and visibility cubes in ps_eor.datacube store brightness temperature in Kelvin. Frequencies are stored in Hz, angular quantities in radians, and Cartesian UV coordinates in wavelengths unless a method says otherwise.

Cartesian visibility cubes#

CartDataCube stores a complex array shaped (frequency, spatial mode) together with its frequency and UV coordinates, image metadata, and optional weights. Saving a cube to HDF5 preserves this information:

from ps_eor import datacube

cube = datacube.CartDataCube.load("visibilities.h5")
cube.save("visibilities-copy.h5")

Selecting frequencies and UV range#

Use get_slice to select an inclusive frequency interval, or get_slice_from_idx for a boolean mask, index array, or Python slice. Both return a new cube. filter_uvrange selects a baseline interval in place:

selected = cube.get_slice(120e6, 130e6)  # frequencies in Hz
selected = selected.get_slice_from_idx(selected.freqs != 125e6)
selected.filter_uvrange(60, 200)          # UV distance in wavelengths

Apply the same selection to related data, noise, and model cubes so their coordinates remain aligned. Methods named filter_* generally modify their cube; use cube.copy() first when the original sampling must be retained.

Creating a cube from FITS images#

For frequency-dependent images in Jy/PSF, the preferred path supplies a matching PSF FITS file for every image. The measured PSF gives a more accurate Jy/PSF-to-Kelvin normalization and visibility weight than inferring the PSF area from image headers alone. Each PSF header must also contain a visibility weight/normalization key: WSCENVIS (preferred), WSCNVIS, or WEIGHT. Without one of these keys a weight cube cannot be constructed.

from pathlib import Path

import numpy as np

from ps_eor import datacube, psutil

image_files = psutil.sort_by_fits_key(
    [str(path) for path in Path("fits/images").glob("*.fits")],
    "CRVAL3",
)
psf_files = psutil.sort_by_fits_key(
    [str(path) for path in Path("fits/psfs").glob("*.fits")],
    "CRVAL3",
)

cube = datacube.CartDataCube.load_from_fits_image_and_psf(
    image_files,
    psf_files,
    umin=50,
    umax=250,
    theta_fov=np.deg2rad(4),
    int_time=10,
    total_time=10 * 3600,
)

Here UV limits are in wavelengths, theta_fov is in radians, and observing times are in seconds. If separate PSF images are unavailable, load_from_fits_image() can infer the PSF normalization from image metadata, but does not provide the same measured PSF weight information.

If the FITS images already contain Kelvin values, disable the Jy/PSF conversion explicitly:

cube = datacube.CartDataCube.load_from_fits_image(
    image_files,
    umin=50,
    umax=250,
    theta_fov=np.deg2rad(4),
    convert_jy2k=False,
)

Creating a cube from an in-memory Kelvin image cube#

An array shaped (frequency, x, y) can be wrapped as a CartImageCube and Fourier transformed to a visibility cube. The input array is already expected to be in Kelvin:

from ps_eor.datacube import CartImageCube, ImageMetaData

# image_data_k has shape (n_freq, nx, ny); freqs_hz has shape (n_freq,).
meta = ImageMetaData.from_res(
    res=np.deg2rad(1.5 / 60),
    shape=image_data_k.shape[-2:],
)
meta.wcs.wcs.cdelt[2] = freqs_hz[1] - freqs_hz[0]

image_cube = CartImageCube(image_data_k, freqs_hz, meta)
cube = image_cube.ft(umin=50, umax=250)

Transforming, changing the field of view, and windowing#

image(), regrid(), and ft() connect the sparse visibility and image representations:

gridded = cube.regrid()       # sparse UV samples -> regular UV grid
image_cube = gridded.image()  # regular UV grid -> Kelvin sky images
round_trip = image_cube.ft(umin=50, umax=250)

For convenience, cube.image() performs the regridding and imaging in one call. To restrict the image-plane field of view or apply a taper while remaining in visibility space:

smaller = cube.reduce_fov(np.deg2rad(3), umin=50, umax=250)
tapered = smaller.apply_window_function(
    datacube.WindowFunction("blackmanharris"),
    umin=50,
    umax=250,
)

Both methods return new cubes and propagate the weights when present. reduce_fov updates the WCS and image shape. apply_window_function records the window name and parameters in the output metadata by default, so the information survives save() and can be used in later power-spectrum normalization. Pass add_to_meta=False only when that bookkeeping is intentionally unwanted.

Stored metadata#

Every image or visibility cube has an ImageMetaData object available as cube.meta. Its properties provide the main observing and image information without requiring direct access to the FITS header:

Property	Information	Obtained from
`res`	Angular pixel size in radians.	Image WCS
`theta_fov`	Image field of view in radians.	Image shape and pixel size
`ra_dec_center_deg`	Phase-centre right ascension and declination in degrees.	Image WCS
`freq_width`	Spacing between frequency samples in Hz.	Frequency WCS
`chan_width`	Physical channel width in Hz.	Channel metadata, or frequency spacing if it is unavailable
`int_time` and `total_time`	Integration time and accumulated observing time in seconds.	FITS metadata or values supplied when creating the cube
`obs_mjd`	Observation reference time as an MJD.	Image WCS
`win_fct`, `win_fct_power`, and `win_fct_area`	Applied image window and its normalizations.	Window metadata; a boxcar window is assumed if none was recorded

This information is copied when new cubes are created and saved with the cube. Operations such as trimming the field of view or applying a window update the corresponding metadata automatically.

Plotting and cube arithmetic#

cube.plot() displays an image-plane view of one channel (or a reduction over channels), while cube.plot_uv() scatters the real visibility values on the UV plane:

cube.plot(fmhz=125, theta_lines=[1, 2])
cube.plot_uv(fmhz=125, uv_lines=[50, 100, 200])
cube.plot_uv(action_fct=np.mean)

Compatible cubes support addition and subtraction, and a cube can be scaled by a number. These operations return a new cube and preserve its geometry:

residual = stokes_i - foreground_model
average = 0.5 * (split_a + split_b)

The operands must have the same frequency channels (and, in practice, the same spatial grid). Cube weights are propagated according to the relevant cube type.

Weights and noise#

A cube’s weights describe sampling or relative statistical weight; they are not a noise realization. Missing or rejected samples should be represented consistently by filtering them from all associated cubes or by assigning zero weight.

When a measured noise realization or noise proxy has weights, it can estimate the system-equivalent flux density (SEFD):

noise_proxy = stokes_v.make_diff_cube()

sefd = noise_proxy.estimate_sefd()
sefd_by_frequency = noise_proxy.estimate_freqs_sefd(
    sefd_poly_fit_deg=3,
)
sefd_by_uv = noise_proxy.estimate_uv_sefd()

estimate_sefd returns one overall value, estimate_freqs_sefd an array over frequency, and estimate_uv_sefd a cube over spatial modes. Values are in Jy by default; pass sefd_jansky=False for Kelvin. The same estimators are available from the weight cube as weights.estimate_*sefd(noise_proxy) when the weights and noise cube are held separately.

Given an SEFD, a CartWeightCube can calculate the expected noise level or draw a realization with the radiometer equation:

total_time = 100 * 3600  # seconds

noise_std = data_cube.weights.get_noise_std_cube(
    sefd_by_frequency,
    time=total_time,
)
simulated_noise = data_cube.weights.simulate_noise(
    sefd_by_frequency,
    time=total_time,
)

NoiseStdCube stores the full complex thermal-noise standard deviation, sqrt(E[abs(noise)**2]). In the example above, noise_std is a NoiseStdCube obtained directly from the visibility weights, a measured SEFD, and the requested accumulated time. Its values are real and non-negative. Use generate_noise_cube() only when an analysis needs a random realization rather than the expected noise level. It is equivalent in purpose to drawing with weights.simulate_noise after the standard-deviation cube has already been constructed.

Combining observations#

DataCubeCombiner incrementally combines observations on their common coordinates and records the accumulated weights and observing time:

from ps_eor.datacube import DataCubeCombiner

combiner = DataCubeCombiner(
    umin=50,
    umax=250,
    weighting_mode="full",
)

for night_id, cube in observations.items():
    selected = cube.copy()
    selected.filter_uvrange(50, 250)
    combiner.add(selected, night_id)

combined = combiner.get(min_n_nights=3)
combined.save("combined_I.h5")

weighting_mode="full" uses the frequency- and UV-dependent sample weights. The alternatives are "uv" for one weight per spatial mode, "global" for one weight per cube, and "none" for an unweighted mean. Set inhomogeneous=True when observations do not all contain the same frequency and UV samples; min_n_nights can then discard channels with insufficient coverage.

Use separate combiners for Stokes I, Stokes V, and time-difference cubes, but feed them the same observation order and quality selections so their output grids remain compatible.