Source code for xga.imagetools.profile

#  This code is a part of X-ray: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS).
#  Last modified by David J Turner ( 13/04/2023, 23:33. Copyright (c) The Contributors

from typing import Tuple

import numpy as np
from astropy.cosmology import Cosmology
from astropy.units import Quantity, UnitBase, pix, deg, arcsec, UnitConversionError

from .misc import pix_deg_scale, pix_rad_to_physical, physical_rad_to_pix, rad_to_ang
from .. import DEFAULT_COSMO
from ..products import Image, RateMap
from ..products.profile import SurfaceBrightness1D

[docs]def annular_mask(centre: Quantity, inn_rad: np.ndarray, out_rad: np.ndarray, shape: tuple, start_ang: Quantity = Quantity(0, 'deg'), stop_ang: Quantity = Quantity(360, 'deg')) -> np.ndarray: """ A handy little function to generate annular (or circular) masks in the form of numpy arrays. It produces the src_mask for a given shape of image, centered at supplied coordinates, and with inner and outer radii supplied by the user also. Angular limits can also be supplied to give the src_mask an annular dependence. This function should be properly vectorised, and accepts inner and outer radii in the form of arrays. The result will be an len_y, len_x, N dimensional array, where N is equal to the length of inn_rad. :param Quantity centre: Astropy pix quantity of the form Quantity([x, y], pix). :param np.ndarray inn_rad: Pixel radius for the inner part of the annular src_mask. :param np.ndarray out_rad: Pixel radius for the outer part of the annular src_mask. :param Quantity start_ang: Lower angular limit for the src_mask. :param Quantity stop_ang: Upper angular limit for the src_mask. :param tuple shape: The output from the shape property of the numpy array you are generating masks for. :return: The generated src_mask array. :rtype: np.ndarray """ # Split out the centre coordinates cen_x = centre[0].value cen_y = centre[1].value # Making use of the astropy units module, check that we are being pass actual angle values if start_ang.unit not in ['deg', 'rad']: raise ValueError("start_angle unit type {} is not an accepted angle unit, " "please use deg or rad.".format(start_ang.unit)) elif stop_ang.unit not in ['deg', 'rad']: raise ValueError("stop_angle unit type {} is not an accepted angle unit, " "please use deg or rad.".format(stop_ang.unit)) # Enforcing some common sense rules on the angles elif start_ang >= stop_ang: raise ValueError("start_ang cannot be greater than or equal to stop_ang.") elif start_ang > Quantity(360, 'deg') or stop_ang > Quantity(360, 'deg'): raise ValueError("start_ang and stop_ang cannot be greater than 360 degrees.") elif stop_ang < Quantity(0, 'deg'): raise ValueError("stop_ang cannot be less than 0 degrees.") else: # Don't want to pass astropy objects to numpy functions, but do need the angles in radians start_ang ='rad').value stop_ang ='rad').value # Check that if the inner and outer radii are arrays, they're the same length if isinstance(inn_rad, (np.ndarray, list)) and len(inn_rad) != len(out_rad): raise ValueError("inn_rad and out_rad are not the same length") elif isinstance(inn_rad, list) and len(inn_rad) == len(out_rad): # If it is a list, just quickly transform to a numpy array inn_rad = np.array(inn_rad) out_rad = np.array(out_rad) # This sets up the cartesian coordinate grid of x and y values arr_y, arr_x = np.ogrid[:shape[0], :shape[1]] # Go to polar coordinates rec_x = arr_x - cen_x rec_y = arr_y - cen_y # Leave this as r**2 to avoid square rooting and involving floats init_r_squared = rec_x**2 + rec_y**2 # arctan2 does just perform arctan on two values, but then uses the signs of those values to # decide the quadrant of the output init_arr_theta = (np.arctan2(rec_x, rec_y) - start_ang) % (2*np.pi) # Normalising to 0-2pi range # If the radius limits are an array, the arrays that describe the space we have constructed are copied # into a third dimension - This allows masks for different radii to be generated in a vectorised fashion if isinstance(inn_rad, np.ndarray): arr_r_squared = np.repeat(init_r_squared[:, :, np.newaxis], len(inn_rad), axis=2) arr_theta = np.repeat(init_arr_theta[:, :, np.newaxis], len(inn_rad), axis=2) else: arr_r_squared = init_r_squared arr_theta = init_arr_theta # This will deal properly with inn_rad and out_rads that are arrays if np.greater(inn_rad, out_rad).any(): raise ValueError("inn_rad value cannot be greater than out_rad") else: rad_mask = (arr_r_squared < out_rad ** 2) & (arr_r_squared >= inn_rad ** 2) # Finally, puts a cut on the allowed angle, and combined the radius and angular cuts into the final src_mask ang_mask = arr_theta <= (stop_ang - start_ang) ann_mask = rad_mask * ang_mask # Should ensure that the central pixel will be 0 for annular masks that are bounded by zero. # Sometimes they aren't because of custom angle choices if 0 in inn_rad: where_zeros = np.where(inn_rad == 0)[0] ann_mask[cen_y, cen_x, where_zeros] = 1 if ann_mask.shape[-1] == 1: ann_mask = np.squeeze(ann_mask) # Returns the annular src_mask(s), in the form of a len_y, len_x, N dimension np array return ann_mask
[docs]def ann_radii(im_prod: Image, centre: Quantity, rad: Quantity, z: float = None, pix_step: int = 1, rad_units: UnitBase = arcsec, cosmo: Cosmology = DEFAULT_COSMO, min_central_pix_rad: int = 3, start_pix_rad: int = 0) -> Tuple[np.ndarray, np.ndarray, Quantity]: """ Will probably only ever be called by an internal brightness calculation, but two different methods need it so it gets its own method. :param Image im_prod: An Image or RateMap product object that you wish to calculate annuli for. :param Quantity centre: The coordinates of the centre of the set of annuli. :param Quantity rad: The outer radius of the set of annuli. :param float z: The redshift of the source of interest, required if the output radius units are a proper radius. :param int pix_step: The width (in pixels) of each annular bin, default is 1. :param UnitBase rad_units: The output units for the centres of the annulli returned by this function. The inner and outer radii will always be in pixels. :param Cosmology cosmo: An instance of an astropy cosmology, the default is a concordance flat LambdaCDM model. :param int start_pix_rad: The pixel radius at which the innermost annulus starts, default is zero. :param int min_central_pix_rad: The minimum radius of the innermost circular annulus (will only be used if start_pix_rad is 0, otherwise the innermost annulus is not a circle), default is three. :return: Returns the inner and outer radii of the annuli (in pixels), and the centres of the annuli in cen_rad_units. :rtype: Tuple[np.ndarray, np.ndarray, Quantity] """ # Quickly get the central coordinates in degrees as well deg_cen = im_prod.coord_conv(centre, deg) pix_cen = im_prod.coord_conv(centre, pix) # If the radius is passed in a proper distance unit, this will convert it or raise an exception # if it can't if rad.unit.is_equivalent('deg') or rad.unit.is_equivalent('kpc'): rad = physical_rad_to_pix(im_prod, rad, deg_cen, z, cosmo) elif rad.unit == pix: # If the radius is already in pixels then all is well pass else: raise UnitConversionError('{} is not a recognised distance unit'.format(rad.unit.to_string())) rad = np.ceil(rad) # So I'm adding a safety feature here, and ensuring the the central circle is a minimum radius if pix_step < 3 and start_pix_rad == 0: central_circ = np.array([0, min_central_pix_rad], dtype=int) ann_rads = np.arange(min_central_pix_rad+pix_step, rad.value + 1, pix_step).astype(int) rads = np.concatenate([central_circ, ann_rads]) else: # By this point, the rad should be in pixels rads = np.arange(start_pix_rad, rad.value + 1, pix_step).astype(int) inn_rads = rads[:len(rads) - 1] out_rads = rads[1:len(rads)] pix_cen_rads = Quantity((inn_rads + out_rads)/2, pix) cen_rads = pix_rad_to_physical(im_prod, pix_cen_rads, rad_units, deg_cen, z, cosmo) # If the innermost radius is zero then the innermost annulus is actually a circle and we don't want # the central radius to be between it and the next radius, as that wouldn't be strictly accurate if rads[0] == 0: cen_rads[0] = Quantity(0, rad_units) return inn_rads, out_rads, cen_rads
[docs]def radial_brightness(rt: RateMap, centre: Quantity, outer_rad: Quantity, back_inn_rad_factor: float = 1.05, back_out_rad_factor: float = 1.5, interloper_mask: np.ndarray = None, z: float = None, pix_step: int = 1, rad_units: UnitBase = arcsec, cosmo: Cosmology = DEFAULT_COSMO, min_snr: float = 0.0, min_central_pix_rad: int = 3, start_pix_rad: int = 0) -> Tuple[SurfaceBrightness1D, bool]: """ A simple method to calculate the average brightness in circular annuli upto the radius of the chosen region. The annuli are one pixel in width, and as this uses the masks that were generated earlier, interloper sources should be removed. :param RateMap rt: A RateMap object to construct a brightness profile from. :param Quantity centre: The coordinates for the centre of the brightness profile. :param Quantity outer_rad: The outer radius of the brightness profile. :param float back_inn_rad_factor: This factor is multiplied by the outer pixel radius, which gives the inner radius for the background mask. :param float back_out_rad_factor: This factor is multiplied by the outer pixel radius, which gives the outer radius for the background mask. :param np.ndarray interloper_mask: A numpy array that masks out any interloper sources. :param float z: The redshift of the source of interest. :param int pix_step: The width (in pixels) of each annular bin, default is 1. :param BaseUnit rad_units: The desired output units for the central radii of the annuli. :param Cosmology cosmo: An astropy cosmology object for source coordinate conversions. :param float min_snr: The minimum signal to noise allowed for each bin in the profile. If any point is below this threshold the profile will be rebinned. Default is 0.0 :param int start_pix_rad: The pixel radius at which the innermost annulus starts, default is zero. :param int min_central_pix_rad: The minimum radius of the innermost circular annulus (will only be used if start_pix_rad is 0, otherwise the innermost annulus is not a circle), default is three. :return: The brightness is returned in a flat numpy array, then the radii at the centre of the bins are returned in units of kpc, the width of the bins, and finally the average brightness in the background region is returned. :rtype: Tuple[SurfaceBrightness1D, bool] """ def _iterative_profile(annulus_masks: np.ndarray, inner_rads: np.ndarray, outer_rads: np.ndarray) \ -> Tuple[bool, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Internal function to calculate and re-bin (ONCE) a surface brightness profile. The profile, along with modified (if rebinned) masks and radii arrays are returned to the user. This can be called once, or iteratively by a loop. :param np.ndarray annulus_masks: 512x512xN numpy array of annular masks, where N is the number of annuli :param np.ndarray inner_rads: The inner radii of the annuli. :param np.ndarray outer_rads: The outer radii of the annuli. :return: Boolean variable that describes whether another re-binning iteration is required, the brightness profile and uncertainties, the modified annular masks, inner radii, outer radii, and annulus areas. :rtype: """ # These are annular masks with interloper sources removed, sensor and edge masks applied corr_ann_masks = annulus_masks * interloper_mask[..., None] * rt.sensor_mask[..., None] \ * rt.edge_mask[..., None] # This calculates the area of each annulus mask num_pix = np.sum(corr_ann_masks, axis=(0, 1)) ann_areas = num_pix * to_arcmin ** 2 # Applying the annular masks, with interlopers removed, and chip edges masked_countrate_data = corr_ann_masks *[..., None] masked_countrate_error_data = corr_ann_masks * count_rate_err_map[..., None] # Defining the brightness profile with the current annuli bright_profile = np.sum(masked_countrate_data, axis=(0, 1)) / ann_areas # And the uncertainties on the profiles. Adding pixels in quadrature bright_profile_errors = np.sqrt(np.sum(masked_countrate_error_data**2, axis=(0, 1))) / ann_areas # We want to check if all parts of the profile are above the defined minimum signal to noise snr_prof = bright_profile / countrate_bg_per_area # Checking if we are below the minimum signal to noise anywhere below = np.argwhere(snr_prof < min_snr).flatten() if below.shape != (0,) and below[0] != (bright_profile.shape[0] - 1): annulus_masks[:, :, below[0]] = annulus_masks[:, :, below[0]] + annulus_masks[:, :, below[0] + 1] annulus_masks = np.delete(annulus_masks, below[0] + 1, axis=2) outer_rads[below[0]] = outer_rads[below[0] + 1] outer_rads = np.delete(outer_rads, below[0] + 1) inner_rads = np.delete(inner_rads, below[0] + 1) another_pass = True elif below.shape != (0,) and below[0] == (bright_profile.shape[0] - 1): annulus_masks[:, :, below[0] - 1] = annulus_masks[:, :, below[0] - 1] + annulus_masks[:, :, below[0]] annulus_masks = np.delete(annulus_masks, below[0], axis=2) outer_rads[below[0] - 1] = outer_rads[below[0]] outer_rads = np.delete(outer_rads, below[0]) inner_rads = np.delete(inner_rads, below[0]) another_pass = True else: another_pass = False return another_pass, bright_profile, bright_profile_errors, annulus_masks, inner_rads, outer_rads, ann_areas if interloper_mask is not None and rt.shape != interloper_mask.shape: raise ValueError("The shape of the src_mask array {0} must be the same as that of im_prod " "{1}.".format(interloper_mask.shape, rt.shape)) elif interloper_mask is None: interloper_mask = np.ones(rt.shape) # Returns conversion factor to degrees, so multiplying by 60 goes to arcminutes # Getting this because we want to be able to convert pixel distance into arcminutes for dividing by the area to_arcmin = pix_deg_scale(centre, rt.radec_wcs).value * 60 # Just making sure we have the centre in pixel coordinates pix_cen = rt.coord_conv(centre, pix) # This sets up the initial annular bin radii, as well as finding the central radii of the bins in the chosen units. inn_rads, out_rads, cen_rads = ann_radii(rt, centre, outer_rad, z, pix_step, rad_units, cosmo, min_central_pix_rad, start_pix_rad) # These calculate the inner and out pixel radii for the background mask - placed in arrays because the # annular mask function expects iterable radii. As pixel radii have to be integer for generating a mask, # I've used ceil. back_inn_rad = np.array([np.ceil(out_rads[-1]*back_inn_rad_factor).astype(int)]) back_out_rad = np.array([np.ceil(out_rads[-1]*back_out_rad_factor).astype(int)]) # Using my annular mask function to make a nice background region, which will be corrected for instrumental # stuff and interlopers in a second back_mask = annular_mask(pix_cen, back_inn_rad, back_out_rad, rt.shape) # Includes chip gaps, interloper removal, and edge removal (to try and avoid artificially bright pixels) # in the background mask corr_back_mask = back_mask * rt.sensor_mask * rt.edge_mask * interloper_mask # Calculates the area of the background region in arcmin^2 back_pix = np.sum(corr_back_mask, axis=(0, 1)) back_area = back_pix * to_arcmin**2 if back_area == 0: raise ValueError("The background mask combined with the sensor mask is is all zeros, this is probably" " because you're looking at a large cluster at a low redshift.") # Per pix is the background which has been divided by the number of pixels count_bg_per_pix = np.sum( * corr_back_mask, axis=(0, 1)) / back_pix # This is where the pixels have been converted into an actual area measurement in arcmin units countrate_bg_per_area = np.sum( * corr_back_mask, axis=(0, 1)) / back_area # Defined here so we can use the where argument in np.sqrt err_calc = ( - count_bg_per_pix) + 2*count_bg_per_pix # Errors on the raw photon counts count_err_map = np.zeros(rt.shape) # Turning those into count rate errors np.sqrt(err_calc, where=err_calc > 0, out=count_err_map) # I divide by a copy of the expmap data array here, otherwise it breaks, and if I'm honest I don't know why... count_rate_err_map = np.divide(count_err_map,, != 0) * rt.edge_mask # Using the ellipse adds enough : to get all the dimensions in the array, then the None adds an empty # dimension. Helps with broadcasting the annular masks with the region src_mask that gets rid of interlopers ann_masks = annular_mask(pix_cen, inn_rads, out_rads, rt.shape) # Make copies of originals in case re-binning fails init_inn = inn_rads.copy() init_out = out_rads.copy() # This gets switched to false if the signal to noise conditions are not fulfilled calculate_profile = True while calculate_profile: calculate_profile, br_prof, br_errs, ann_masks, inn_rads, \ out_rads, areas = _iterative_profile(ann_masks, inn_rads, out_rads) if len(inn_rads) == 0: # This means that rebinning failed, and the profile couldn't be made greater than the minimum signal to noise inn_rads = init_inn out_rads = init_out ann_masks = annular_mask(pix_cen, inn_rads, out_rads, rt.shape) # Set the min_snr to 0 to get an non-rebinned profile min_snr = 0 prof_results = _iterative_profile(ann_masks, inn_rads, out_rads) br_prof, br_errs = prof_results[1:3] areas = prof_results[-1] succeeded = False else: succeeded = True final_inn_rads = pix_rad_to_physical(rt, Quantity(inn_rads, pix), rad_units, centre, z, cosmo) final_out_rads = pix_rad_to_physical(rt, Quantity(out_rads, pix), rad_units, centre, z, cosmo) # Need these simply because the brightness profile must always be aware of the radii in degrees in order to # to assemble a storage key properly deg_inn_rads = pix_rad_to_physical(rt, Quantity(inn_rads, pix), deg, centre, z, cosmo) deg_out_rads = pix_rad_to_physical(rt, Quantity(out_rads, pix), deg, centre, z, cosmo) cen_rads = (final_inn_rads + final_out_rads) / 2 deg_cen_rads = (deg_inn_rads + deg_out_rads) / 2 if final_inn_rads[0].value == 0: cen_rads[0] = Quantity(0, cen_rads.unit) deg_cen_rads[0] = Quantity(0, 'deg') rad_err = (final_out_rads-final_inn_rads) / 2 if not outer_rad.unit.is_equivalent('deg'): deg_outer_rad = rad_to_ang(outer_rad, z, cosmo).to('deg') else: deg_outer_rad ="deg") # Now I've finally implemented some profile product classes I can just smoosh everything into a convenient product br_prof = SurfaceBrightness1D(rt, cen_rads, Quantity(br_prof, 'ct/(s*arcmin**2)'), centre, pix_step, min_snr, deg_outer_rad, rad_err, Quantity(br_errs, 'ct/(s*arcmin**2)'), Quantity(countrate_bg_per_area, 'ct/(s*arcmin**2)'), np.insert(out_rads, 0, inn_rads[0]), np.concatenate([back_inn_rad, back_out_rad]), Quantity(areas, 'arcmin**2'), deg_cen_rads, succeeded) return br_prof, succeeded
[docs]def pizza_brightness(im_prod: Image, src_mask: np.ndarray, back_mask: np.ndarray, centre: Quantity, rad: Quantity, num_slices: int = 4, z: float = None, pix_step: int = 1, cen_rad_units: UnitBase = arcsec, cosmo=DEFAULT_COSMO) -> Tuple[np.ndarray, Quantity, Quantity, np.float64, np.ndarray, np.ndarray]: raise NotImplementedError("The supporting infrastructure to allow pizza profile product objects hasn't been" " written yet sorry!")