:py:mod:`aces.holography.clean_holography` ========================================== .. py:module:: aces.holography.clean_holography .. autoapi-nested-parse:: Clean holography data Module Contents --------------- Classes ~~~~~~~ .. autoapisummary:: aces.holography.clean_holography.BeamModel Functions ~~~~~~~~~ .. autoapisummary:: aces.holography.clean_holography.tft aces.holography.clean_holography.report_flags aces.holography.clean_holography.robust_poly_fit aces.holography.clean_holography.fit_spec aces.holography.clean_holography.fit_all_spectra aces.holography.clean_holography.add_metadata_to_header aces.holography.clean_holography.write_to_fits aces.holography.clean_holography.write_tt_to_fits aces.holography.clean_holography.interp2D_cmplx aces.holography.clean_holography.normalise aces.holography.clean_holography.volt_to_power aces.holography.clean_holography.feed_to_stokes aces.holography.clean_holography.raw_to_stokes aces.holography.clean_holography.get_quality_flags aces.holography.clean_holography.form_mean aces.holography.clean_holography.save_mean aces.holography.clean_holography.centre_beams aces.holography.clean_holography.tt aces.holography.clean_holography.get_beampos_errs aces.holography.clean_holography.get_shift aces.holography.clean_holography.centre_mass aces.holography.clean_holography.get_beampos_errs_from_fits aces.holography.clean_holography.shifts_to_mapset aces.holography.clean_holography.save_shifts aces.holography.clean_holography.design_mat aces.holography.clean_holography.pfit aces.holography.clean_holography.get_posfit_flags aces.holography.clean_holography.set_resid_lim aces.holography.clean_holography.get_shift_flags aces.holography.clean_holography.main aces.holography.clean_holography.cli Attributes ~~~~~~~~~~ .. autoapisummary:: aces.holography.clean_holography.log .. py:data:: log .. py:function:: tft(t) .. py:function:: report_flags(flags, tag) .. py:function:: robust_poly_fit(spec, poly_ord=2, sigma_lim=0.8, max_iter=20, low_fact=10.0) Fit a polynomial to data with outliers The test for the presence of outliers is a comparison between standard deviation (std) and the inter-quartile range (IQR). The IQR is a robust measure of the distribution width, and for a gaussian distribution std/IQR = approx 0.6. The value of std/IQR at which to cease the search for outliers is sigma_lim. :param spec: Spectrum to fit :type spec: np.ndarray :param poly_ord: Order of polynomial, defaults to 5 :type poly_ord: int, optional :param sigma_lim: Sigma limit, defaults to 0.8 :type sigma_lim: float, optional :param max_iter: Maxium interations, defaults to 20 :type max_iter: int, optional :return: TODO :rtype: TODO .. py:function:: fit_spec(spec, poly_ord=5, sigma_lim=3.0, max_iter=50) Fit a polynomial to data with outliers :param spec: Spectrum to fit :type spec: np.ndarray :param poly_ord: Order of polynomial, defaults to 5 :type poly_ord: int, optional :param sigma_lim: Sigma limit, defaults to 3.0 :type sigma_lim: float, optional :param max_iter: Maxium interations, defaults to 50 :type max_iter: int, optional :return: TODO :rtype: TODO .. py:function:: fit_all_spectra(hm, poly_ord=3) Determines the best fitting models and generates that model for a given mean cube Args: hm ([array]): data cube (masked array) poly_ord (int): Returns: [array]: Returns the best fitting model over all pixels and the associated frequencies. .. py:function:: add_metadata_to_header(header: astropy.io.fits.header.Header, metadata: dict[str, Any]) -> None Common function to add metadata from a MapSet into the FITS header. Some items from the MapSet metadata structure are recorded as a list of values. Here these are handled by taking the first element of the list. :param header: The HDU header that is being added to :type header: fits.header.Header :param metadata: The key-value pair of items to add to header :type metadata: dict .. py:function:: write_to_fits(data_clean: numpy.ndarray, mapset: aces.beamset.mapset.MapSet, interp_factor: float, output_dir: str) Write data cube to FITS: :param data_clean: Cleaned array (with Stokes axis). :type data_clean: np.ndarray :param mapset: Input holography dataset :type mapset: MapSet :param interp_factor: Intperpolation grid used in deg :type interp_factor: float :param output_dir: Directory to save file to :type output_dir: str .. py:function:: write_tt_to_fits(beamtt: numpy.ndarray, mapset: aces.beamset.mapset.MapSet, fref: float, bw: float, dxy: float, output_dir: float) Write Taylor Terms to FITS: :param beamtt: Taylor term beams :type beamtt: np.ndarray :param mapset: Input holography dataset :type mapset: MapSet :param fref: Reference frequency for Taylor terms :type fref: float :param bw; The observing bandwidth that is recorded to the header :type bw: float :param dxy: Interpolation grid in deg :type dxy: float :param output_dir: Directory to save file to :type output_dir: str .. py:function:: interp2D_cmplx(z, pixpos, k=2) :param z: (ndarray) Complex cube of shape Nt,Na,Nb,Np,Nf,nx,ny :param pixpos: (ndarray) Array of shape (nbeams,2) giving beam positions in pixels :param k: (int) Order of spline interpolation used. :return interp (ndarray) Interpolated values at beam positions - shape Nt,Na,Nb,Np,Nf .. py:function:: normalise(v, bpos, refant) Given an array of raw holography visibilities v, normalise to give the response of each antenna relative to its response to a point source at beam centre. Each map value at beam centre is determined by interpolation: vxx0, vyy0. The normalised quantities are computed as Nxx = vxx/vxx0, Nxy = vxy/vxx0, Nyx = vyx/vyy0, Nyy = vyy/vyy0 :param v:(numpy.ndarray) visibilities :param bpos: (np.ndarray) beam position in pixels (nb,2) :param refant: (int) Reference antenna number (0 based) :return N: (numpy.ndarray) normalised visibilities NOTE: This scheme works provided the reference antenna has a lower index than all target antennas For cases where the correlation is Target * Reference we need: scale_xy = scale_yy scale_yx = scale_xx .. py:function:: volt_to_power(v) Given a Jones matrix v, return the corresponding brightness matrix b as b = v * v(h) where (h) is the hermitian transpose. In both v and b, the polarisation axis is the 4th of 7 axes: sh = (nt,ne,nb,4,nf,nx,ny) :param v: (numpy.ndarray) voltage vector of shape sh :return b: (numpy.ndarray) brightness vector of shape sh .. py:function:: feed_to_stokes(f) Given a matrix (voltage or power) in the feed frame, return the corresponding Stokes matrix. Do this for ASKAP, assuming the measurements were made with the antenna rotated +45 degrees. In both f and s, the polarisation axis is the 4th of 7 axes: sh = (nt,ne,nb,4,nf,nx,ny) :param f: (numpy.ndarray) vector in feed frame: XX,XY,YX,YY :return s: (numpy.ndarray) vector in Stokes frame: I,Q,U,V .. py:function:: raw_to_stokes(raw_file_name, refant) Reads holography data in hdf5 MapSet format. Performs operations on the data hyper-cube: * transforms to axes aligned with celestial coords * determines map values at beam positions in Vxx, Vyy * derives normalising factors for Vxx, Vxy, Vyx, Vyy visibilities * normalises the complex data * determines brightness vector XX, XY, YX, YY * converts to Stokes I Q U V * saves both Stokes vector and scale values in hdf5 MapSet format :param raw_file_name: (str) Name of hdf5 holography data file :param refant: (int) Reference antenna :return stokes_obj (MapSet) Holography maps, normalised and in Stokes brightness I Q U V :return centre_obj (MapSet) Normalising scales - visibilities at beam positions (interpolated). .. py:function:: get_quality_flags(mapset, amp_max=1.4) Take MapSet object holding Holography power patterns in the IQUV Stokes parameters; look for and flag image planes whose maximum stokes I value significantly exceeds the expected value of 2.0. Typically, the large value flags arise from planes that were effected by RFI and produce large values at random locations away from the beam peak: :param mapset: MapSet object holding cube of beam images :type mapset: aces.beamset.mapset.MapSet :param amp_max: Maximum allowed amplitude :type amp_max: float, optional :return: result - a dict holding flags and correspondingly masked data :rtype: dict .. py:function:: form_mean(mapset) Take MapSet object holding Holography power patterns in the IQUV Stokes parameters; Compute the mean over the antenna axis for all beams, polarisations and spectral points using the MapSet flags to include only image planes judged to be of high quality: :param mapset: Object holding cube of beam images :type mapset: aces.beamset.mapset.MapSet :return: mean_cube: Mean cube of shape [nB, nP, nF, nx, ny] :rtype: numpy.ndarray .. py:function:: save_mean(mean_cube, mapset, out_name) Takes the array-wide mean cube and writes it to an hdf5 file in the BeamSet format :param mean_cube: Data cube of shape (nb, npol, nf, nx, ny) :param mapset: BeamSet object from full Stokes cube; used to derive new metadata for output :param out_name: Name of output file :type out_name: str :return: .. py:function:: centre_beams(mean_cube, mapset, kxy=5, nxy=101) Takes the mean holography cube (averaged over antennas) and interpolates each beam onto a grid centred on its nominal position and with the given grid spacing :param mean_cube: shape Nb,Np,Nf,Nx,Ny :type mean_cube: np.ndarray :param mapset: parent object - for metadata :type mapset: [type] :param kxy: Degrees of bivariate spline used in RecBivariateSpline, defaults to 5 :type kxy: int, optional :param nxy: Size of output grid (must be odd integer)., defaults to 101 :type nxy: int, optional :return: TODO :rtype: TODO .. py:function:: tt(fstart, beam_cube, nterm) For correction of continuum images we turn the frequency axis into a taylor term axis. - Note this requires picking a centre frequency f0, for which the Taylor coefficients are computed :param fstart: Starting frequency in MHz :type fstart: float :param beam_cube: Beam maps, averaged over antennas :type beam_cube: np.ndarray :param nterm: Number of Taylor terms :type nterm: int :return: TODO :rtype: TODO .. py:class:: BeamModel(bs_obj) Bases: :py:obj:`object` Class - to be improved 1. hold data in hdf5 file, not pickle. .. py:method:: get(pixoff, freq) TODO :param pixoff: TODO :type pixoff: TODO :param freq: TODO :type freq: TODO :return: TODO :rtype: TODO .. py:method:: _radial_gau_gen() :staticmethod: Parameters used here were derived from a 'good set' of mid-band beams (SBID=16881), the Gaussian parameter (0.6945) minimises the difference between model gaussian and measured average beam shape at 1438.5 MHz. .. py:function:: get_beampos_errs(mapset, bmodel) TODO :param mapset: TODO :type mapset: TODO :param bmodel: TODO :type bmodel: TODO :return: TODO :rtype: TODO .. py:function:: get_shift(mapset, ref_im, target_im, frac=0.25) Compute the shift of target_im relative to ref_im. Use Fourier shift theorem The sense of the returned shifts: If the returned value shift = (dl, dm), and if the position of reference and target beams are (lr, mr) and (lt, mt) respectively, then dl > 0.0 indicates lt > lr and if dm > 0.0 then tm > rm. Coordinate $l$ increases towards the west; $m$ increases towards the north. :param mapset: Primary mapset :type mapset: [type] :param ref_im: Reference image ndim=2 :type ref_im: np.ndarray :param target_im: Target image :type target_im: np.ndarray :param frac: Mask transform values < ampMax * frac in phase fit, defaults to 0.25 :type frac: float, optional :return: TODO :rtype: TODO .. py:function:: centre_mass(im, xaxis) .. py:function:: get_beampos_errs_from_fits(fits_header, fits_data) .. py:function:: shifts_to_mapset(mapset, shifts, resids, amps) Package the beam shift data into a Mapset object :param mapset: Input mapset object, holdng the holography maps :param shifts: Stokes I beam position offsets :param resids: Residuals to offset fit :param amps: :return: Mapset object holdong shift data .. py:function:: save_shifts(beamset, holo_dir) Save a BeamSet object holding beam errors to an hdf5 file: :param beamset: BeamSet object holding beam pos errors :type beamset: BeamSet :param holo_dir: Directory to receive hdf5 file :type holo_dir: str .. py:function:: design_mat(x, op) .. py:function:: pfit(x, y, op) .. py:function:: get_posfit_flags(resid, amps, res_limits, amp_limits) TODO :param resid: TODO :type resid: TODO :param amps: TODO :type amps: TODO :return: n_flags New flags of shape to match obj.flags; True for bad data :rtype: TODO .. py:function:: set_resid_lim(resids, order=2, iqr_fac=2.0) .. py:function:: get_shift_flags(shifts, shift_fit_flags) TODO :param shifts: :param shift_fit_flags: :return: n_flags New flags of shape to match obj.flags; True for bad data :rtype: TODO .. py:function:: main(sbid, holo_dir='.', remake_stokes=False, max_order=5, param='BIC', max_iter=50, snr_limit=3) Main script :param sbid: SBID :type sbid: int :param holo_dir: Directory containing holo data, defaults to '.' :type holo_dir: str, optional :param remake_stokes: Remake Stokes parameters, defaults to False :type remake_stokes: bool, optional :param max_order: Maximum polynomial to fit, defaults to 5 :type max_order: int, optional :param param: Criterion to evaluate fit, defaults to 'BIC' :type param: str, optional :param max_iter: Maximum interations for fitting, defaults to 50 :type max_iter: int, optional :param snr_limit: SNR limit for fitting, defaults to 3 :type snr_limit: int, optional .. py:function:: cli()