Data handling
=============

All image and visibility cubes in :mod:`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
--------------------------

:class:`~ps_eor.datacube.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:

.. code-block:: python

   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:

.. code-block:: python

   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.

.. code-block:: python

   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,
:meth:`~ps_eor.datacube.CartDataCube.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:

.. code-block:: python

   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
:class:`~ps_eor.datacube.CartImageCube` and Fourier transformed to a
visibility cube. The input array is already expected to be in Kelvin:

.. code-block:: python

   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:

.. code-block:: python

   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:

.. code-block:: python

   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 :meth:`~ps_eor.datacube.CartDataCube.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
:class:`~ps_eor.datacube.ImageMetaData` object available as ``cube.meta``.
Its properties provide the main observing and image information without
requiring direct access to the FITS header:

.. list-table::
   :header-rows: 1
   :widths: 28 42 30

   * - 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:

.. code-block:: python

   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:

.. code-block:: python

   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:

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):

.. code-block:: python

   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 :class:`~ps_eor.datacube.CartWeightCube` can calculate the
expected noise level or draw a realization with the radiometer equation:

.. code-block:: python

   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,
   )

:class:`~ps_eor.datacube.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
:meth:`~ps_eor.datacube.NoiseStdCube.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:

Combining observations
----------------------

:class:`~ps_eor.datacube.DataCubeCombiner` incrementally combines observations
on their common coordinates and records the accumulated weights and observing
time:

.. code-block:: python

   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.