12. API Reference

You can find automatically generated documentation for each module in the nemo package below. Note that some functions listed here may not be used by any of the scripts found in the bin/ directory.

12.1. catalogs

This module contains tools for handling catalogs, which are usually astropy.table.Table objects.

catalog2DS9(catalog, outFileName, constraintsList=[], addInfo=[], idKeyToUse='name', RAKeyToUse='RADeg', decKeyToUse='decDeg', color='cyan', showNames=True, writeNemoInfo=True, coordSys='fk5')

Writes a DS9 region file corresponding to the given catalog.

Parameters
  • catalog (astropy.table.Table) – An astropy Table where each row represents an object.

  • outFileName (str) – A file name for the output DS9 region file.

  • constraintsList (list, optional) – A list of constraints in the same format as used by selectFromCatalog().

  • addInfo (list, optional) – A list of dictionaries with keys named key and fmt (e.g., {'key': "SNR", 'fmt': "%.3f"}). These will be added to the object label shown in DS9.

  • idKeyToUse (str, optional) – The name of the key in each object dictionary that defines the object’s name. Used to label objects in the DS9 region file.

  • RAKeyToUse (str, optional) – The name of the key in each object dictionary that contains the RA of the object in decimal degrees.

  • decKeyToUse (str, optional) – The name of the key in each object dictionary that contains the declination of the object in decimal degrees.

  • color (str, optional) – The color of the plot symbol used by DS9.

  • writeNemoInfo (bool, optional) – If True, writes a line with the nemo version and date generated at the top of the DS9 .reg file.

  • coordSys (str, optional) – A string defining the coordinate system used for RA, dec, as understood by DS9.

Returns

None

catalogListToTab(catalogList, keysToWrite=['name', 'RADeg', 'decDeg', 'SNR', 'numSigPix', 'template', 'tileName', 'galacticLatDeg', 'deltaT_c', 'err_deltaT_c', 'y_c', 'err_y_c', 'Y500_sr', 'err_Y500_sr', 'fluxJy', 'err_fluxJy', 'redshift', 'redshiftErr', 'ellipse_PA', 'ellipse_A', 'ellipse_B', 'ellipse_x0', 'ellipse_y0', 'ellipse_e', 'fixed_deltaT_c', 'fixed_err_deltaT_c', 'fixed_y_c', 'fixed_err_y_c'])

Converts a catalog in the form of a list of dictionaries (where each dictionary holds the object properties) into an astropy Table.

Parameters
  • catalogList (list) – Catalog in the form of a list of dictionaries.

  • keysToWrite (list, optional) – Keys to convert into columns in the output table.

Returns

An astropy Table object.

crossMatch(refCatalog, matchCatalog, radiusArcmin=2.5)

Cross matches matchCatalog onto refCatalog for objects found within some angular radius (specified in arcmin).

Parameters
  • refCatalog (astropy.table.Table) – The reference catalog.

  • matchCatalog (astropy.table.Table) – The catalog to match onto the reference catalog.

  • radiusArcmin (float, optional) – Cross-match radius in arcmin.

Returns

Cross-matched reference catalog, matchCatalog, and array of angular separation in degrees, for objects in common within the matching radius. The cross matched columns are sorted such that rows in each correspond to the matched objects.

generateRandomSourcesCatalog(mapData, wcs, numSources)

Generate a random source catalog (with amplitudes in deltaT uK), with random positions within the footprint of the given map (areas where pixel values == 0 are ignored). The distribution of source amplitudes is roughly similar to that seen in the 148 GHz ACT maps, but this routine should only be used for tests - it is not a serious attempt at simulating the real extragalactic source population.

Parameters
  • mapData (numpy.ndarray) – Map pixel-data, only used for determining valid area in which sources may be randomly placed (pixel values == 0 are ignored).

  • wcs (astWCS.WCS) – WCS corresponding to the map.

  • numSources (int) – Number of random sources to put into the output catalog.

Returns

An astropy Table object containing the catalog.

generateTestCatalog(config, numSourcesPerTile, amplitudeColumnName='fixed_y_c', amplitudeRange=[0.001, 1], amplitudeDistribution='linear', selFn=None, avoidanceRadiusArcmin=20.0)

Generate a catalog of objects with random positions and amplitudes. This is for testing purposes - see, e.g., nemo.maps.sourceInjectionTest().

Parameters
  • config (nemo.startup.NemoConfig) – Nemo configuration object.

  • numSourcesPerTile (int) – The maximum number of sources to insert into each tile. The number of sources actually inserted may be less than this depending on the value of avoidanceRadiusArcmin.

  • amplitudeColumnName (str) – Name of the column in the output catalog in which source (or cluster) amplitudes will be stored. Typically this should be “deltaT_c” for sources, and “fixed_y_c” for clusters.

  • amplitudeRange (list) – Range for the random amplitudes, in the form [minimum, maximum].

  • amplitudeDistribution (str) – Either ‘linear’ or ‘log’.

  • selFn (nemo.completeness.SelFn, optional) – Nemo selection function object, used to access area masks and coordinates info. If not given, a selFn object will be created using the info in the config. Providing this saves time, as the area mask files don’t have to be read from disk.

  • avoidanceRadiusArcmin (float) – Minimum separation between two objects in the output catalog. This should be set large enough to avoid crowding and spurious cross-matching problems.

Returns

An astropy Table object containing the catalog.

getCatalogWithinImage(tab, shape, wcs)

Returns the subset of the catalog with coordinates within the image defined by the given shape and wcs.

Parameters
  • tab (astropy.table.Table) – Catalog, as an astropy Table object. Must have columns called ‘RADeg’, ‘decDeg’ that contain object coordinates in decimal degrees.

  • shape (list) – Shape of the array corresponding to the image / map.

  • wcs (astWCS.WCS) – WCS of the image.

Returns

An astropy Table containing the subset of objects within the image.

getTableRADecKeys(tab)

Returns the column names in the table in which RA, dec coords are stored, after trying a couple of variations.

Parameters

tab (astropy.table.Table) – The table to search.

Returns

Name of RA column, name of dec. column

makeLongName(RADeg, decDeg, prefix='ACT-CL')

Makes a long format object name string from the given object coordinates, following IAU convention.

Parameters
  • RADeg (float) – Right ascension of the object in J2000 decimal degrees.

  • decDeg (float) – Declination of the object in J2000 decimal degrees.

  • prefix (str, optional) – Prefix for the object name.

Returns

Object name string in the format prefix JHHMMSS.s+/-DDMMSS.

makeName(RADeg, decDeg, prefix='ACT-CL')

Makes an object name string from the given object coordinates, following the IAU convention.

Parameters
  • RADeg (float) – Right ascension of the object in J2000 decimal degrees.

  • decDeg (float) – Declination of the object in J2000 decimal degrees.

  • prefix (str, optional) – Prefix for the object name.

Returns

Object name string in the format prefix JHHMM.m+/-DDMM.

makeOptimalCatalog(catalogDict, constraintsList=[])

Identifies common objects between every catalog in the input dictionary of catalogs, and creates a master catalog with one entry per object, keeping only the details of the highest signal-to-noise detection.

Parameters
  • catalogDict (dict) – Dictionary where each key points to a catalog of objects.

  • constraintsList (list, optional) – A list of constraints (for the format, see selectFromCatalog()).

Returns

None - an optimalCatalog key is added to catalogDict in place.

removeCrossMatched(refCatalog, matchCatalog, radiusArcmin=2.5)

Cross matches matchCatalog onto refCatalog for objects found within some angular radius (specified in arcmin), and returns refCatalog with the matching entries removed.

Parameters
  • refCatalog (astropy.table.Table) – The reference catalog.

  • matchCatalog (astropy.table.Table) – The catalog to match onto the reference catalog.

  • radiusArcmin (float, optional) – Cross-match radius in arcmin.

Returns

Cross-matched reference catalog (astropy.table.Table) with matches to matchCatalog removed.

removeDuplicates(tab)

Given an astropy Table object, remove duplicate objects - keeping the highest SNR detection for each duplicate. This routine is used to clean up output of MPI runs (where we have overlapping tiles). This method could be applied elsewhere, if we re-did how we handled catalogs everywhere in nemo.

Returns table with duplicates removed, number of duplicates found, names of duplicates

selectFromCatalog(catalog, constraintsList)

Return a table of objects matching the given constraints from the catalog.

Parameters
  • catalog (astropy.table.Table) – The catalog from which objects will be selected.

  • constraintsList (list) – A list of constraints, where each item is a string of the form “key < value”, “key > value”, etc.. Note that the spaces between the key, operator (e.g. ‘<’), and value are essential.

Returns

An astropy Table object.

tabToCatalogList(tab)

Converts an astropy Table into a list of dictionaries.

Returns catalog list

writeCatalog(catalog, outFileName, constraintsList=[])

Writes the catalog (astropy Table) to disk. The output format depends upon the extension given by outFileName - any format supported by the astropy Table object can be used.

NOTE: For .csv format, tab is used as the delimiter. So, to read these using astropy you will need to use something like: tab=atpy.Table().read(“test.csv”, delimiter = ‘ ‘).

constraintsList works as in the selectFromCatalog function.

12.2. completeness

This module contains tools for calculating the completeness of cluster samples.

class SelFn(selFnDir, SNRCut, configFileName=None, footprintLabel=None, zStep=0.01, tileNames=None, enableDrawSample=False, downsampleRMS=True, applyMFDebiasCorrection=True, setUpAreaMask=False, enableCompletenessCalc=True)

An object that describes a survey selection function, and provides routines for calculating survey completeness for a given set of scaling relation and cosmological parameters.

SNRCut

Completeness will be computed relative to this signal-to-noise selection cut (labelled as fixed_SNR in the catalogs produced by nemo).

Type

float

footprintLabel

Use this to specify a footprint, if any are defined in the Nemo config used to produce the selFn dir (e.g., ‘DES’, ‘HSC’, ‘KiDS’ etc.). The default of None uses the whole survey footprint.

Type

str

applyMFDebiasCorrection

Set to False to disable the Eddington bias correction of mass estimates. Probably only useful for debugging.

Type

bool

selFnDir

Path to a selFn/ directory, as produced by the nemo and nemoSelFn commands. This directory contains information such as the survey noise maps, area masks, and filter mismatch function (Q) fits for clusters.

Type

str

zStep

Use this to set the binning in redshift for completeness calculations.

Type

float

tileNames

The list of tiles used by the SelFn object (default of None uses all tiles).

Type

list

WCSDict

A dictionary indexed by tileName, containing WCS objects that describe the mapping between pixel coords and (RA, dec) coords in each tile.

Type

dict

areaMaskDict

A dictionary containing the survey area masks, indexed by tileName. Values > 0 in these masks define the cluster or source search area.

Type

dict

scalingRelationDict

A dictionary of scaling relation parameters (see example Nemo config files for the format).

Type

dict

tckQFitDict

A dictionary of interpolation spline knots indexed by tileName, that can be used to estimate Q, the filter mismatch function (see nemo.signals.loadQ()).

Type

dict

RMSDict

A dictionary of RMS tables, indexed by tileName. Each RMSTable contains noise level by area, as returned by getRMSTab.

Type

dict

totalAreaDeg2

The total area, as measured from the survey mask, for the given set of tiles and footprint.

Type

float

fRelDict

A dictionary of weights used for relativistic corrections, indexed by tileName.

Type

dict

mockSurvey (

obj: MockSurvey): A Nemo MockSurvey object, used for halo mass function calculations and generating mock catalogs.

Initialise an object that contains a survey selection function.

This class uses the output in the selFn/ directory (produced by the nemo and nemoSelFn commands) to re-calculate the survey completeness for a given signal-to-noise cut on a (log10 M500c, z) grid for a given set of cosmological and scaling relation parameters.

Parameters
  • selFnDir (str) – Path to a selFn/ directory, as produced by the nemo and nemoSelFn commands. This directory contains information such as the survey noise maps, area masks, and filter mismatch function (Q) fits for clusters.

  • SNRCut (float) – Completeness will be computed relative to this signal-to-noise selection cut (labelled as fixed_SNR in the catalogs produced by nemo).

  • configFileName (str, optional) – Path to a Nemo configuration file. If not given, this will be read from selFnDir/config.yml, which is the config file for the nemo run that produced the selFn directory.

  • footprintLabel (str, optional) – Use this to specify a footprint, if any are defined in the Nemo config used to produce the selFn dir (e.g., ‘DES’, ‘HSC’, ‘KiDS’ etc.). The default of None uses the whole survey footprint.

  • zStep (float, optional) – Use this to set the binning in redshift for completeness calculations.

  • tileNames (list, optional) – If given, restrict the SelFn object to use only these tiles.

  • enableDrawSample (bool, optional) – This only needs to be set to True for generating mock catalogs.

  • downsampleRMS (float, optional) – Downsample the resolution of the RMS tables by this factor. The RMS tables are generated from the noise maps, and are just a listing of noise level versus survey area. Downsampling speeds up completeness calculations considerably.

  • applyMFDebiasCorrection (bool, optional) – Set to False to disable the Eddington bias correction of mass estimates. Probably only useful for debugging.

  • setupAreaMask (bool, optional) – If True, read in the area masks so that quick position checks can be done (e.g., by checkCoordsAreInMask).

  • enableCompletenessCalc (bool, optional) – If True, set up the machinery needed to do completeness calculations.

checkCoordsInAreaMask(RADeg, decDeg)

Checks if the given RA, dec coords are in valid regions of the map.

Parameters
  • RADeg (float or numpy array) – RA in decimal degrees

  • decDeg (float or numpy array) – dec in decimal degrees

Returns

True if in the area mask mask, False if not

overrideRMS(RMS, obsFreqGHz=148.0)

Override the RMS noise of the SelFn object - replacing it with constant value RMS (given in uK/arcmin^2). This is translated into a y0~ noise level at the given observing frequency.

After doing this call update() to get an updated (M, z) completeness grid.

projectCatalogToMz(tab)

Project a catalog (as astropy Table) into the (log10 M500, z) grid. Takes into account the uncertainties on y0, redshift - but if redshift error is non-zero, is a lot slower.

Returns (log10 M500, z) grid

update(H0, Om0, Ob0, sigma_8, scalingRelationDict=None)

Re-calculates survey-average selection function given new set of cosmological / scaling relation parameters.

This updates self.mockSurvey at the same time - i.e., this is the only thing that needs to be updated.

Resulting (M, z) completeness grid stored as self.compMz

To apply the selection function and get the expected number of clusters in the survey do e.g.:

selFn.mockSurvey.calcNumClustersExpected(selFn = selFn.compMz)

(yes, this is a bit circular)

calcCompleteness(RMSTab, SNRCut, tileName, mockSurvey, scalingRelationDict, tckQFitDict, plotFileName=None, z=None, method='fast', numDraws=2000000, numIterations=100)

Calculate completeness as a function of (log10 M500c, z) on the mockSurvey grid at the given SNRCut. Intrinsic scatter in the scaling relation is taken into account.

Parameters
  • RMSTab (astropy Table) – Table containing noise level by area, as returned by getRMSTab.

  • SNRCut (float) – Completeness will be calculated for objects relative to this cut in fixed_SNR.

  • tileName (str) – Name of the map tile.

  • mockSurvey (MockSurvey) – A Nemo MockSurvey object used for halo mass function calculations.

  • scalingRelationDict (dict) – A dictionary of scaling relation parameters (see example Nemo config files for the format).

  • tckQFitDict (dict) – A dictionary, with keys corresponding to tile names, containing the spline knots used to obtain interpolated values from the filter mismatch function Q (see nemo.signals.loadQ()).

  • plotFileName (str, optional) – If given, write a plot showing 90% completness limit.

  • z (array, optional) – Redshift grid on which the completeness calculation will be performed. Alternatively, a single redshift can be specified as a float instead.

  • method (str, optional) – Two methods for doing the calculation are available: “fast” (applies the measurement errors and scatter to ‘true’ y0 values on a grid) and “montecarlo” (uses samples drawn from a mock catalog, generated on the fly, to estimate the completeness). Both methods should give consistent results.

  • numDraws (int, optional) – Used by the “montecarlo” method - sets the number of draws from the halo mass function on each iteration.

  • numIterations (int, optional) – Used by the “montecarlo” method - sets the number of iterations, i.e., the number of mock catalogs from which the completeness is estimated.

Returns

2d array of (log10 M500c, z) completeness

calcMassLimit(completenessFraction, compMz, mockSurvey, zBinEdges=[])

Given a completeness (log10M, z) grid as made by calcCompleteness, return the mass limit (units of 10^14 MSun) as a function of redshift. By defualt, the same binning used in the given mockSurvey object - this can be overridden by giving zBinEdges.

calcTileWeightedAverageNoise(tileName, photFilterLabel, selFnDir, diagnosticsDir=None, footprintLabel=None)

Returns weighted average noise value in the tile.

completenessByFootprint(selFnCollection, mockSurvey, diagnosticsDir, additionalLabel='')

Write out average (M, z) grid for all tiles (tileNames) given in selFnCollection (a dictionary with keys corresponding to footprints: ‘full’ is the entire survey), weighted by fraction of total survey area within the footprint. We also produce a bunch of other stats and plots to do with completeness versus redshift.

Output is written to files named e.g. diagnosticsDir/MzCompleteness_label.npz, where label is the footprint (key in selFnCollection); ‘full’ is the default (survey-wide average).

additionalLabel is optional and will be added to the output filenames (use for e.g., tagging with the calcCompleteness method used).

cumulativeAreaMassLimitPlot(z, diagnosticsDir, selFnDir, tileNames)

Make a cumulative plot of 90%-completeness mass limit versus survey area, at the given redshift.

downsampleRMSTab(RMSTab, stepSize=1.0000000000000001e-07)

Downsamples the RMSTab in terms of noise resolution, binning by stepSize.

Returns RMSTab

getRMSTab(tileName, photFilterLabel, selFnDir, diagnosticsDir=None, footprintLabel=None)

Makes a table containing map area in tile pointed to by tileName against RMS values (so this compresses the information in the RMS maps). The first time this is run takes ~200 sec (for a 1000 sq deg tile), but the result is cached.

If making the RMSTab, diagnosticsDir must be given. If the RMSTab exists, it is not needed.

Can optionally take extra masks for specifying e.g. HSC footprint. Here, we assume these have already been made by makeIntersectionMask, and we can load them from the cache, identifying them through footprintLabel

Returns RMSTab

getTileTotalAreaDeg2(tileName, selFnDir, masksList=[], footprintLabel=None)

Returns total area of the tile pointed at by tileName (taking into account survey mask and point source masking).

A list of other survey masks (e.g., from optical surveys like KiDS, HSC) can be given in masksList. These should be file names for .fits images where 1 defines valid survey area, 0 otherwise. If given, this routine will return the area of intersection between the extra masks and the SZ survey.

loadAreaMask(tileName, selFnDir)

Loads the survey area mask (i.e., after edge-trimming and point source masking, produced by nemo) for the given tile.

Returns map array, wcs

loadIntersectionMask(tileName, selFnDir, footprint)

Loads the intersection mask for the given tile and footprint.

Returns map array, wcs

loadMassLimitMap(tileName, diagnosticsDir, z)

Loads the mass limit map for the given tile at the given z.

Returns map array, wcs

loadRMSMap(tileName, selFnDir, photFilter)

Loads the RMS map for the given tile.

Returns map array, wcs

makeFullSurveyMassLimitMapPlot(z, config)

Reprojects tile mass limit maps onto the full map pixelisation, then makes a plot and saves a .fits image.

makeIntersectionMask(tileName, selFnDir, label, masksList=[])

Creates intersection mask between mask files given in masksList.

Calculating the intersection is slow (~30 sec per tile per extra mask for KiDS), so intersection masks are cached; label is used in output file names (e.g., diagnosticsDir/intersect_label#tileName.fits)

Can optionally be called without extraMasksList, IF the intersection mask has already been created and cached.

NOTE: We assume masks have dec aligned with y direction for speed.

Returns intersectionMask as array (1 = valid area, 0 = otherwise)

makeMassLimitMap(SNRCut, z, tileName, photFilterLabel, mockSurvey, scalingRelationDict, tckQFitDict, diagnosticsDir, selFnDir)

Makes a map of 90% mass completeness (for now, this fraction is fixed).

NOTE: The map here is “downsampled” (i.e., binned) in terms of noise resolution - okay for display purposes, may not be for analysis using the map (if anyone should decide to). Much quicker and saves on memory issues (the y0 noise estimates have uncertainties themselves anyway).

makeMassLimitVRedshiftPlot(massLimit_90Complete, zRange, outFileName, title=None)

Write a plot of 90%-completeness mass limit versus z, adding a spline interpolation.

makeMzCompletenessPlot(compMz, log10M, z, title, outFileName)

Makes a (M, z) plot. Here, compMz is a 2d array, and log10M and z are arrays corresponding to the axes.

tidyUp(config)

Tidy up the selFn directory - making MEF files etc. and deleting individual tile files.

We also copy the configuration file into the selFn dir to make things easier for distribution/end users.

12.3. filters

This module defines several filter classes, together with a function that uses them to filter maps.

There are two main classes of filter: MatchedFilter and RealSpaceMatchedFilter.

New base classes can be derived from the overall base class MapFilter. There are also base classes corresponding to filters with different signal templates (BeamFilter, ArnaudModelFilter). The actual filters that can be used are derived from these, e.g.:

class ArnaudModelFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Base class for filters using the GNFW profile as described in Arnaud et al. (2010).

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

makeSignalTemplateMap(beamFileName, amplitude=None)

Makes a model signal template map.

class ArnaudModelMatchedFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

class ArnaudModelRealSpaceMatchedFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

class BeamFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Base class for filters using beam profile files in Matthew + Kavi’s format.

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

makeSignalTemplateMap(beamFileName, amplitude=None)

Makes a beam model signal template map.

class BeamMatchedFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

class BeamRealSpaceMatchedFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

class MapFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Generic map filter base class. Defines common interface.

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

buildAndApply()

Builds and applies the filter to the unfiltered map(s). Returns a dictionary, containing keys ‘data’, ‘wcs’, ‘weights’, ‘obsFreqGHz’. If this routine converts to yc in place, the latter is set to ‘yc’.

loadFRelWeights()

Reads frequency weights used for relativistic corrections from filter header.

makeForegroundsPower()

Returns a Power2D object with foregrounds power from the CMB (this could be extended to add other foregrounds if necessary).

makeNoiseMap(mapData)

Estimate the noise map using local RMS measurements on a grid, over the whole filtered map.

makeRadiansMap()

Makes a map of distance in radians from centre, for the map being filtered.

makeRealSpaceFilterProfile()

Makes a 1d real-space profile of the filter. This is normalised to 1 at the frequency of the first map given in the unfiltered maps list.

Returns profile, arcminRange

saveRealSpaceFilterProfile()

Saves a real-space profile of the filter out as a .png plot and 2d .fits image.

class MatchedFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Multi-frequency matched filter…

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

applyFilter(mapDataToFilter)

Apply the filter to the given map data (must be a cube - with each plane corresponding to a frequency). If the map data is not complex, it will be Fourier transformed. If the map data is not the same shape as the filter, the filter will be interpolated to match.

An optional additional high-pass filter can be applied if ‘bckSub’ and ‘bckSubScaleArcmin’ are given in self.params.

Returns

Filtered map (2d numpy array)

buildAndApply()

Builds and applies the filter to the unfiltered map(s). Returns a dictionary, containing keys ‘data’, ‘wcs’, ‘weights’, ‘obsFreqGHz’. If this routine converts to yc in place, the latter is set to ‘yc’.

loadFilter()

Loads in a previously saved filter.

reshapeFilter(shape)

Use linear interpolation to transform the filter to the given shape.

Returns

Reshaped filter (2d numpy array)

class RealSpaceMatchedFilter(label, unfilteredMapsDictList, paramsDict, tileName='PRIMARY', writeFilter=False, forceRebuild=False, diagnosticsDir=None, selFnDir=None)

Makes a matched-filter kernel using the noise properties of a specified region of the map (e.g. the deepest part) in Fourier space, which is converted into real space and truncated such that the kernel is small enough in footprint to be applied by directly convolving with the map in real space in a short amount of time.

First though we high-pass filter the maps on a similar scale using a Mexican hat (difference of Gaussian smoothed images), to get rid of large scale noise from CMB and atmosphere.

Tne S/N map is constructed, measuring locally on a scale 3 x that of the filter kernel, with 3-sigma clipping applied. This allows us to take into account how the noise level varies across the map.

A rank filter is used to zap noisy regions at the edge of the map, where the RMS values are not accurate.

Optionally, a ‘surveyMask’ can be applied (e.g., for point source masking).

A map of the final area searched for clusters called ‘areaMask.fits’ is written in the selFn/ folder.

Initialises a MapFilter. unfilteredMapsDictList describes the input maps, paramsDict describes the filter options. The convention is that single frequency only filters (e.g. GaussianWienerFilter) only operate on the first map in the unfilteredMapDictList.

label is used to store the filter to save calculating it again. Filters are stored under filters/ dir. If a filter already exists under this dir with the given label, it is loaded rather than recalculated. To force recalculating the filter, set forceRebuild == True.

applyFilter(mapDataToFilter, calcFRelWeights=False)

Apply the kernel to the given map data (must be a cube - with each plane corresponding to a frequency).

NOTE: calcFRelWeights = True should ONLY be set if this routine is being applied to an ideal signal map (when calculating signalNorm). This operation is being done in here to save complications / time from doing unnecessary convolution operations elsewhere.

Returns

Filtered map (2d numpy array)

buildAndApply()

Builds and applies the filter to the unfiltered map(s). Returns a dictionary, containing keys ‘data’, ‘wcs’, ‘weights’, ‘obsFreqGHz’. If this routine converts to yc in place, the latter is set to ‘yc’.

buildKernel(RADecSection, RADeg='centre', decDeg='centre')

Builds the real space kernel itself.

RADeg, decDeg are used for figuring out pixel scales for background subtraction

Returns kern2d, signalNorm, bckSubScaleArcmin

loadFilter()

Loads a previously cached kernel.

Returns kern2d, signalNorm, bckSubScaleArcmin

filterMaps(unfilteredMapsDictList, filterParams, tileName, filteredMapsDir='.', diagnosticsDir='.', selFnDir='.', verbose=True, undoPixelWindow=True, useCachedMaps=True)

Build and applies filters to the unfiltered maps(s). The output is a filtered map in yc or uK (this can be set with outputUnits in the config file). All filter operations are done in the filter objects, even if multifrequency (a change from previous behaviour).

Filtered maps are written to rootOutDir/filteredMaps Filters, if stored, are written to rootOutDir/filters

Returns a dictionary containing the filtered map and S/N map etc.

12.4. gnfw

Perform line-of-sight integral of GNFW model.

Matthew Hasselfield’s code, slightly updated to include P0, c500 (makes life a bit easier for some applications of this).

func(x, params)

The GNFW radial profile - now with added P0, c500.

integrated(b, params={'P0': 8.403, 'alpha': 1.051, 'beta': 5.4905, 'c500': 1.177, 'gamma': 0.3081, 'npts': 100, 'tol': 1e-07})

Return the line of sight integral of the GNFW profile at impact parameter b. Since the GNFW is a smoothly varying power-law from 0+ to +infty, we make a change of variables to

u = ln x

and perform the integral by summing over equally spaced logarithmic bins in x.

xfunc(x, b, params)

The log-scaled integrand for GNFW cylindrical integration along line-of-sight coordinate x, with impact parameter b.

12.5. maps

This module contains tools for manipulating maps (e.g., conversion of units etc.).

addAutoTileDefinitions(parDict, DS9RegionFileName=None, cacheFileName=None)

Runs autotiler to add automatic tile definitions into the parameters dictionary in-place.

Parameters
  • parDict (dict) – Dictionary containing the contents of the Nemo config file.

  • DS9RegionFileName (str, optional) – Path to DS9 regions file to be written.

  • cacheFileName (str, optional) – Path to output a cached .yml file which will can be read instead on repeated runs (for speed).

addWhiteNoise(mapData, noisePerPix)

Adds Gaussian distributed white noise to mapData.

applyPointSourceMask(maskFileName, mapData, mapWCS, mask=0.0, radiusArcmin=2.8)

Given file name pointing to a point source mask (as made by maskOutSources), apply it to given mapData.

autotiler(surveyMask, wcs, targetTileWidth, targetTileHeight)

Given a survey mask (where values > 0 indicate valid area, and 0 indicates area to be ignored), figure out an optimal tiling strategy to accommodate tiles of the given dimensions. The survey mask need not be contiguous (e.g., AdvACT and SO maps, using the default pixelization, can be segmented into three or more different regions).

Parameters
  • surveyMask (numpy.ndarray) – Survey mask image (2d array). Values > 0 will be taken to define valid area.

  • wcs (astWCS.WCS) – WCS associated with survey mask image.

  • targetTileWidth (float) – Desired tile width, in degrees (RA direction for CAR).

  • targetTileHeight (float) – Desired tile height, in degrees (dec direction for CAR).

Returns

Dictionary list defining tiles in same format as config file.

Note

While this routine will try to match the target file sizes, it may not match exactly. Also, makeTileDir will expand tiles by a user-specified amount such that they overlap.

checkMask(fileName)

Checks whether a mask contains negative values (invalid) and throws an exception if this is the case.

Parameters

fileName (str) – Name of the .fits format mask file to check

convertToDeltaT(mapData, obsFrequencyGHz=148)

Converts mapData (in yc) to deltaT (micro Kelvin) at given frequency.

convertToY(mapData, obsFrequencyGHz=148)

Converts mapData (in delta T) at given frequency to y.

convolveMapWithBeam(data, wcs, beam, maxDistDegrees=1.0)

Convolves map defined by data, wcs with the beam.

Parameters
  • data (numpy.ndarray) – Map to convolve, as 2d array.

  • wcs (astWCS.WCS) – WCS corresponding to data (i.e., the map).

  • beam (BeamProfile or str) – Either a BeamProfile object, or a string that gives the path to a text file that describes the beam profile.

  • maxDistDegrees (float) – Sets the size of the convolution kernel, for optimization purposes.

Returns

Beam-convolved map (numpy array).

Note

The pixel scale used to define the convolution kernel is evaluated at the central map pixel. So, this routine should only be used with either pixelisations where the scale is constant or on relatively small tiles.

estimateContamination(contamSimDict, imageDict, SNRKeys, label, diagnosticsDir)

Performs the actual contamination estimate, makes output under diagnosticsDir.

Use label to set a prefix for output (plots / .fits tables), e.g., label = “skySim”

estimateContaminationFromInvertedMaps(config, imageDict)

Run the whole filtering set up again, on inverted maps.

Writes a DS9. reg file, which contains only the highest SNR contaminants (since these are most likely to be associated with artefacts in the map - e.g., point source masking).

Writes a plot and a .fits table to the diagnostics dir.

Runs over both SNR and fixed_SNR values.

Returns a dictionary containing the results

estimateContaminationFromSkySim(config, imageDict)

Estimate contamination by running on source-free sky simulations (CMB plus noise that we generate here on the fly).

This uses the same kernels that were constructed and used on the real maps. The whole filtering and object detection pipeline is run on the simulated maps repeatedly. The number of sky sims used (set by numSkySims in the .yml config file) should be fairly large (~100) for the results to be robust (results on individual sims can vary by a lot).

Parameters
  • config (startUp.NemoConfig) – Nemo configuration object.

  • imageDict (dict) – A dictionary containing the output filtered maps and catalogs from running on the real data (i.e., the output of pipelines.filterMapsAndMakeCatalogs). This will not be modified, but is used for estimating the contamination rate by comparison to the source-free sims.

Returns

A dictionary where each key points to an astropy Table object containing the average contamination estimate corresponding to SNR (maximal estimate) and fixed_SNR (for the chosen reference filter scale).

getPixelAreaArcmin2Map(mapData, wcs)

Returns a map of pixel area in arcmin2

loadTile(pathToTileImages, tileName, returnWCS=False)

Given a path to a directory full of tiles, or a .fits file, return the map array and (optionally) the WCS.

Parameters
  • pathToTileImages (str) – Path to either a .fits file, or a directory full of .fits files named by tileName.

  • tileName (str) – The name of the tile to load.

Returns

Map data (and optionally wcs)

makeModelImage(shape, wcs, catalog, beamFileName, obsFreqGHz=None, GNFWParams='default', cosmoModel=None, applyPixelWindow=True, override=None)

Make a map with the given dimensions (shape) and WCS, containing model clusters or point sources, with properties as listed in the catalog. This can be used to either inject or subtract sources from real maps.

Parameters
  • shape (tuple) – The dimensions of the output map (height, width) that will contain the model sources.

  • wcs (astWCS.WCS) – A WCS object that defines the coordinate system of the map.

  • catalog (astropy.table.Table) – An astropy Table object containing the catalog. This must include columns named ‘RADeg’, ‘decDeg’ that give object coordinates. For point sources, the amplitude in uK must be given in a column named ‘deltaT_c’. For clusters, either ‘M500’ (in units of 10^14 MSun), ‘z’, and ‘fixed_y_c’ must be given (as in a mock catalog), OR the catalog must contain a ‘template’ column, with templates named like, e.g., Arnaud_M1e14_z0p2 (for a z = 0.2, M500 = 1e14 MSun cluster; see the example .yml config files included with nemo).

  • beamFileName – Path to a text file that describes the beam.

  • obsFreqGHz (float, optional) – Used only by cluster catalogs - if given, the returned map will be converted into delta T uK, assuming the given frequency. Otherwise, a y0 map is returned.

  • GNFWParams (str or dict, optional) – Used only by cluster catalogs. If ‘default’, the Arnaud et al. (2010) Universal Pressure Profile is assumed. Otherwise, a dictionary that specifies the profile parameters can be given here (see gnfw.py).

  • override (dict, optional) – Used only by cluster catalogs. If a dictionary containing keys {‘M500’, ‘redshift’} is given, all objects in the model image are forced to have the corresponding angular size. Used by positionRecoveryTest().

  • applyPixelWindow (bool, optional) – If True, apply the pixel window function to the map.

Returns

Map containing injected sources.

makeTileDir(parDict)

Update this later. Revised version. Instead of making a MEF, makes a directory for each map and puts individual tile images there. We’ll only need to edit preprocessMapDict to handle the difference. Why are we doing this? Just in case reading from the same file is gumming up MPI runs on hippo when using lots of nodes.

Makes a tileDir (directory containing .fits images, one per tile) file, if the needed parameters are given in parDict. Adjusts unfilteredMapsDictList accordingly and returns it.

If the options for making a tileDir image aren’t given in parDict, then we pass through a standard single extension file (or rather the path to it, as originally given)

NOTE: If the map given in unfilteredMaps is 3d (enki gives I, Q, U as a datacube), then this will extract only the I (temperature) part and save that in the tileDir file. This will need changing if hunting for polarized sources…

Returns unfilteredMapsDictList [input for filterMaps], list of extension names

NOTE: Under MPI, this should only be called by the rank = 0 process

maskOutSources(mapData, wcs, catalog, radiusArcmin=7.0, mask=0.0, growMaskedArea=1.0)

Given a mapData array and a catalog of source positions, replace the values at the object positions in the map within radiusArcmin with replacement values. If mask == ‘whiteNoise’, this will be white noise with mean and sigma set by the pixel values in an annulus of 1 < r < 2 * radiusArcmin.

growMaskedArea sets factor larger than radiusArcmin to set masked area to in returned mask. This can avoid any weird artefacts making it into source lists.

Returns a dictionary with keys ‘data’ (mapData with mask applied), ‘mask’ (0-1 mask of areas masked).

noiseBiasAnalysis(sourceInjTable, plotFileName, clipPercentile=99.7, sourceInjectionModel=None)

Estimate the noise bias from the ratio of input to recovered flux as a function of signal-to-noise.

Parameters
  • posRecTable (astropy.table.Table) – Table containing recovered position offsets versus fixed_SNR for various cluster/source models (produced by sourceInjectionTest).

  • plotFileName (str) – Path where the plot file will be written.

  • clipPercentile (float, optional) – Clips offset values outside of this percentile of the whole position offsets distribution, to remove a small number of outliers (spurious next-neighbour cross matches) that otherwise bias the contours high for large (99%+) percentile cuts in individual fixed_SNR bins.

  • sourceInjectionModel (str, optional) – If given, restrict analysis to only objects matching this.

Notes

For clusters, bear in mind this only makes sense if any mismatch between the inserted cluster’s shape and the signal assumed by the filter is taken into account. This is done using the Q-function in sourceInjectionTest.

plotContamination(contaminTabDict, diagnosticsDir)

Makes contamination rate plots, output stored under diagnosticsDir

While we’re at it, we write out a text file containing interpolated values for e.g., 5%, 10% contamination levels

positionRecoveryAnalysis(posRecTable, plotFileName, percentilesToPlot=[50, 90, 95, 99], plotRawData=True, rawDataAlpha=0.02, pickleFileName=None, clipPercentile=99.7)

Estimate and plot position recovery accuracy as function of fixed filter scale S/N (fixed_SNR), using the contents of posRecTable (see positionRecoveryTest).

Parameters
  • posRecTable (astropy.table.Table) – Table containing recovered position offsets versus fixed_SNR for various cluster/source models (produced by sourceInjectionTest).

  • plotFileName (str) – Path where the plot file will be written.

  • percentilesToPlot (list, optional) – List of percentiles to plot (some interpolation will be done).

  • plotRawData (bool, optional) – Plot the raw (fixed_SNR, positional offset) data in the background.

  • pickleFileName (string, optional) – Saves the percentile contours data as a pickle file if not None. This is saved a dictionary with top-level keys named according to percentilesToPlot.

  • clipPercentile (float, optional) – Clips offset values outside of this percentile of the whole offsets distribution, to remove a small number of outliers (spurious next-neighbour cross matches) that otherwise bias the contours high for large (99%+) percentile cuts in individual fixed_SNR bins.

preprocessMapDict(mapDict, tileName='PRIMARY', diagnosticsDir=None)

Applies a number of pre-processing steps to the map described by mapDict, prior to filtering.

The first step is to load the map itself and the associated weights. Some other operations that may be applied are controlled by keys added to mapDict. Some of these may be specified in the .yml configuration file, while others are applied by particular filter objects or by routines that generate simulated data. The following keys are understood:

surveyMask (str)

Path to a mask (.fits image; 1 = valid, 0 = masked) that defines the valid object search area.

pointSourceMask (str)

Path to a mask (.fits image; 1 = valid, 0 = masked) that contains holes at the locations of point sources, defining regions that are excluded from the object search area.

RADecSection (list)

Defines a region to extract from the map. Use the format [RAMin, RAMax, decMin, decMax] (units: decimal degrees).

CMBSimSeed (int)

If present, replace the map with a source-free simulated CMB realisation, generated using the given seed number. Used by estimateContaminationFromSkySim().

applyBeamConvolution (:obj: str)

If True, the map is convolved with the beam given in the beamFileName key. This should only be needed when using preliminary y-maps made by tILe-C.

Parameters
  • mapDict (dict) – A dictionary with the same keys as given in the unfilteredMaps list in the .yml configuration file (i.e., it must contain at least the keys ‘mapFileName’, ‘units’, and ‘weightsFileName’, and may contain some of the optional keys listed above).

  • tileName (str) – Name of the map tile (extension name) to operate on.

  • diagnosticsDir (str) – Path to a directory where miscellaneous diagnostic data are written.

Returns

A dictionary with keys that point to the map itself (‘data’), weights (‘weights’), masks (‘surveyMask’, ‘pointSourceMask’), and WCS object (‘wcs’).

saveFITS(outputFileName, mapData, wcs, compressed=False)

Writes a map (2d image array) to a new .fits file.

Parameters
  • outputFileName (str) – Filename of output FITS image.

  • mapData (np.ndarray) – Map data array

  • wcs (astWCS.WCS) – Map WCS object

  • compressed (bool) – If True, writes a compressed image

saveTilesDS9RegionsFile(parDict, DS9RegionFileName)

Writes a DS9 .reg file containing the locations of tiles defined in parDict.

Parameters
  • parDict (dict) – Dictionary containing the contents of the Nemo config file.

  • DS9RegionFileName (str) – Path to DS9 regions file to be written.

shrinkWCS(origShape, origWCS, scaleFactor)

Given an astWCS object and corresponding image shape, scale the WCS by scaleFactor. Used for making downsampled quicklook images (using stitchMaps).

Parameters
  • origShape (tuple) – Shape of the original image.

  • origWCS (astWCS.WCS object) – WCS for the original image.

  • scaleFactor (float) – The factor by which to scale the image WCS.

Returns

shape (tuple), WCS (astWCS.WCS object)

simCMBMap(shape, wcs, noiseLevel=0.0, beamFileName=None, seed=None)

Generate a simulated CMB map, optionally convolved with the beam and with (white) noise added.

Parameters
  • shape (tuple) – A tuple describing the map (numpy array) shape in pixels (height, width).

  • wcs (astWCS.WCS) – An astWCS object.

  • noiseLevel (numpy.ndarray or float) – If a single number, this is taken as sigma (in map units, usually uK) for generating white noise that is added across the whole map. Alternatively, an array with the same dimensions as shape may be used, specifying sigma (in map units) per corresponding pixel. Noise will only be added where non-zero values appear in noiseLevel.

  • beamFileName (str) – The file name of the text file that describes the beam with which the map will be convolved. If None, no beam convolution is applied.

  • seed (int) – The seed used for the random CMB realisation.

Returns

A map (numpy.ndarray)

smoothMap(data, wcs, RADeg='centre', decDeg='centre', smoothScaleDeg=0.08333333333333333)

Smoothes map with Gaussian of given scale.

If RADeg, decDeg = ‘centre’, then the pixel scales used to set the kernel shape will be set from that at the centre of the WCS. Otherwise, they will be taken at the given coords.

Note that wcs is only used to figure out the pixel scales here.

sourceInjectionTest(config, writeRankTable=True)

Insert sources with known positions and properties into the map, apply the filter, and record their offset with respect to the true location as a function of S/N (for the fixed reference scale only). If the inserted sources are clusters, the Q function will be applied to the output fluxes, to account for any mismatch between the reference filter scale and the inserted clusters.

Writes output to the diagnostics/ directory.

Parameters
  • config (nemo.startUp.NemoConfig) – Nemo configuration object.

  • writeRankTable (bool, optional) – If True, saves a table as output for this MPI rank under the diagnostics/ directory. Useful for MPI debugging only.

Returns

An astropy Table containing recovered position offsets and fluxes versus fixed_SNR for inserted sources.

stitchTiles(filePattern, outFileName, outWCS, outShape, fluxRescale=1.0)

Fast routine for stitching map tiles back together. Since this uses interpolation, you probably don’t want to do analysis on the output - this is just for checking / making plots etc.. This routine sums images as it pastes them into the larger map grid. So, if the zeroed (overlap) borders are not handled, correctly, this will be obvious in the output.

NOTE: This assumes RA in x direction, dec in y direction (and CAR projection).

NOTE: This routine only writes output if there are multiple files that match filePattern (to save needless duplicating maps if nemo was not run in tileDir mode).

Output map will be multiplied by fluxRescale (this is necessary if downsampling in resolution).

Takes 10 sec for AdvACT S16-sized downsampled by a factor of 4 in resolution.

subtractBackground(data, wcs, RADeg='centre', decDeg='centre', smoothScaleDeg=0.5)

Smoothes map with Gaussian of given scale and subtracts it, to get rid of large scale power.

If RADeg, decDeg = ‘centre’, then the pixel scales used to set the kernel shape will be set from that at the centre of the WCS. Otherwise, they will be taken at the given coords.

Note that wcs is only used to figure out the pixel scales here.

12.6. MockSurvey

This module defines the MockSurvey class, used for obtaining de-biased cluster mass estimates and selection function calculations.

12.7. nemoCython

convolve()

Convolves image arr with kernel kern. Hopefully faster than ndimage.convolve.

Returns 2d array

makeDegreesDistanceMap()

Fills (in place) the 2d array degreesMap with distance in degrees from the given position, out to some user-specified maximum distance.

Parameters
  • degreesMap (np.ndarray) – Map (2d array) that will be filled with angular distance from the given coordinates. Probably you should feed in an array set to some extreme initial value (e.g., 1e6 everywhere) to make it easy to filter for pixels near the object coords afterwards.

  • wcs (astWCS.WCS) – WCS corresponding to degreesMap.

  • RADeg (float) – RA in decimal degrees of position of interest (e.g., object location).

  • decDeg (float) – Declination in decimal degrees of position of interest (e.g., object location).

  • maxDistDegrees – The maximum radius out to which distance will be calculated.

Returns

A map (2d array) of distance in degrees from the given position, (min x, max x) pixel coords corresponding to maxDistDegrees box, (min y, max y) pixel coords corresponding to maxDistDegrees box

Note

This routine measures the pixel scale local to the given position, then assumes that it does not change. So, this routine may only be accurate close to the given position, depending upon the WCS projection used.

makeLocalSNMap()

Makes local S/N map - measuring noise in an annulus using the standard deviation. Measures background using mean in annulus. Since there is no sorting involved, this is much quicker than using rank filters. But it may not give quite the same results.

makeXYDegreesDistanceMaps()

Returns an array of distance along x, y axes in degrees from given position. maxDistDegrees sets the limit around the given RADeg, decDeg position to which the distance is calculated.

rankFilter()

Rank is relative, with 0 lowest, 1 max.

standardDeviationFilter()

Standard deviation filter using pixels within rPix.

12.8. photometry

This module contains source finding and photometry routines.

addFreqWeightsToCatalog(imageDict, photFilter, diagnosticsDir)

Add relative weighting by frequency for each object in the optimal catalog, extracted from the data cube saved under diagnosticsDir (this is made by makeSZMap in RealSpaceMatchedFilter). This is needed for multi-frequency cluster finding / analysis - for, e.g., weighting relativistic corrections to y0~ (which was estimated from the inverse variance weighted average of y0~ from each frequency map).

NOTE: this is only applied for the reference filter (pointed to by photFilter)

deltaTToJySr(temp, obsFreqGHz)

Convert delta T (uK) to Jy/sr at the given frequency in GHz

findObjects(filteredMapDict, threshold=3.0, minObjPix=3, rejectBorder=10, findCenterOfMass=True, removeRings=True, invertMap=False, objIdent='ACT-CL', longNames=False, verbose=True, useInterpolator=True, measureShapes=False, DS9RegionsPath=None)

Finds objects in the filtered maps pointed to by the imageDict. Threshold is in units of sigma (as we’re using S/N images to detect objects). Catalogs get added to the imageDict.

Throw out objects within rejectBorder pixels of edge of frame if rejectBorder != None

Set invertMap == True to do a test of estimating fraction of spurious sources

Set measureShapes == True to fit ellipses in a similar style to how SExtractor does it (may be useful for automating finding of extended sources).

This now find rings around extremely bright sources and adds them to the survey mask if removeRings is True

getObjectPositions(mapData, threshold, findCenterOfMass=True)

Creates a segmentation map and find objects above the given threshold.

Parameters
  • mapData (numpy.ndarray) – The 2d map to segment.

  • threshold (float) – The threshold above which objects will be selected.

  • findCenterOfMass – If True, return the object center weighted according to the values in mapData. If False, return the pixel that holds the maximum value.

Returns

Array of object ID numbers. objPositions (list): List of corresponding (y, x) positions. objNumPix (numpy.ndarray): Array listing number of pixels per object. segmentationMap (numpy.ndarray): The segmentation map (2d array).

Return type

objIDs (numpy.ndarray)

getPixelsDistanceMap(objDict, data)

Returns an array of same dimensions as data but containing radial distance to the object specified in the objDict in units of pixels

getRadialDistanceMap(objDict, data, wcs)

Returns an array of same dimensions as data but containing radial distance to the object specified in the objDict in units degrees on the sky

getSNRValues(catalog, SNMap, wcs, useInterpolator=True, invertMap=False, prefix='')

Measures SNR values in given map at catalog positions.

Set invertMap == True, to do a test for estimating the spurious source fraction

Use prefix to add prefix before SNR in catalog (e.g., fixed_SNR)

Use template to select a particular filtered map (i.e., label of mapFilter given in .par file, e.g., a filter used to measure properties for all objects at fixed scale).

makeAnnulus(innerScalePix, outerScalePix)

Makes the footprint of an annulus for feeding into the rank filter

Returns annulus, numValues

makeForcedPhotometryCatalog(filteredMapDict, inputCatalogFileName, useInterpolator=True, DS9RegionsPath=None)

Make a catalog on which we can do forced photometry at the locations of objects in it.

measureFluxes(catalog, filteredMapDict, diagnosticsDir, photFilteredMapDict=None, useInterpolator=True)

Add flux measurements to each catalog pointed to in the imageDict. Measured in ‘outputUnits’ specified in the filter definition in the .par file (and written to the image header as ‘BUNIT’).

Maps should have been normalised before this stage to make the value of the filtered map at the object position equal to the flux that is wanted (yc, uK or Jy/beam). This includes taking out the effect of the beam.

Use ‘photFilter’ to choose which filtered map in which to make flux measurements at fixed scale (fixed_delta_T_c, fixed_y_c etc.). Set to None if you don’t want these.

12.9. pipelines

This module defines pipelines - sets of tasks in nemo that we sometimes want to do on different inputs (e.g., real data or simulated data).

filterMapsAndMakeCatalogs(config, rootOutDir=None, copyFilters=False, measureFluxes=True, invertMap=False, verbose=True, useCachedMaps=True)

Runs the map filtering and catalog construction steps according to the given configuration.

Parameters
  • ( (config) – obj: ‘startup.NemoConfig’): Nemo configuration object.

  • rootOutDir (str) – If None, use the default given by config. Otherwise, use this to override where the output filtered maps and catalogs are written.

  • copyFilters (bool, optional) – If True, and rootOutDir is given (not None), then filters will be copied from the default output location (from a pre-existing nemo run) to the appropriate directory under rootOutDir. This is used by, e.g., contamination tests based on sky sims, where the same kernels as used on the real data are applied to simulated maps. If rootOutDir = None, setting copyKernels = True has no effect.

  • measureFluxes (bool, optional) – If True, measure fluxes. If False, just extract S/N values for detected objects.

  • invertMap (bool, optional) – If True, multiply all maps by -1; needed by :meth:maps.estimateContaminationFromInvertedMaps).

Returns

Optimal catalog (keeps the highest S/N detection when filtering at multiple scales).

Note

See bin/nemo for how this pipeline is applied to real data, and maps.estimateContaminationFromSkySim for how this is applied to source-free sims that are generated on the fly.

makeMockClusterCatalog(config, numMocksToMake=1, combineMocks=False, writeCatalogs=True, writeInfo=True, verbose=True)

Generate a mock cluster catalog using the given nemo config.

Returns

List of catalogs (each is an astropy Table object)

makeSelFnCollection(config, mockSurvey)

Makes a collection of selection function dictionaries (one per footprint specified in selFnFootprints in the config file, plus the full survey mask), that contain information on noise levels, area covered, and completeness.

Returns a dictionary (keys: ‘full’ - corresponding to whole survey, plus other keys named by footprint).

12.10. plotSettings

This module contains global plot settings.

For any routine that makes a plot using matplotlib, call plotSettings.update_rcParams() first.

update_rcParams(dict={})

Based on Cristobal’s preferred settings. Updates matplotlib rcParams in place.

12.11. signals

This module contains routines for modeling cluster and source signals.

calcDz(zIn)

Calculate linear growth factor, normalised to D(z) = 1.0 at z = 0.

calcFRel(z, M500, cosmoModel, obsFreqGHz=148.0)

Calculates relativistic correction to SZ effect at specified frequency, given z, M500 in MSun.

This assumes the Arnaud et al. (2005) M-T relation, and applies formulae of Itoh et al. (1998)

As for H13, we return fRel = 1 + delta_SZE (see also Marriage et al. 2011)

calcM500Fromy0(y0, y0Err, z, zErr, tckQFit, mockSurvey, tenToA0=4.95e-05, B0=0.08, Mpivot=300000000000000.0, sigma_int=0.2, applyMFDebiasCorrection=True, applyRelativisticCorrection=True, calcErrors=True, fRelWeightsDict={148.0: 1.0})

Returns M500 +/- errors in units of 10^14 MSun, calculated assuming a y0 - M relation (default values assume UPP scaling relation from Arnaud et al. 2010), taking into account the steepness of the mass function. The approach followed is described in H13, Section 3.2.

Here, mockSurvey is a MockSurvey object. We’re using this to handle the halo mass function calculations (in turn using the Colossus module). Supplying mockSurvey is no longer optional (and handles setting the cosmology anyway when initialised or updated).

tckQFit is a set of spline knots, as returned by fitQ.

If applyMFDebiasCorrection == True, apply correction that accounts for steepness of mass function.

If applyRelativisticCorrection == True, apply relativistic correction (weighted by frequency using the contents of fRelWeightsDict).

If calcErrors == False, error bars are not calculated, they are just set to zero.

fRelWeightsDict is used to account for the relativistic correction when y0~ has been constructed from multi-frequency maps. Weights should sum to 1.0; keys are observed frequency in GHz.

Returns dictionary with keys M500, M500_errPlus, M500_errMinus

calcPM500(y0, y0Err, z, zErr, tckQFit, mockSurvey, tenToA0=4.95e-05, B0=0.08, Mpivot=300000000000000.0, sigma_int=0.2, applyMFDebiasCorrection=True, applyRelativisticCorrection=True, fRelWeightsDict={148.0: 1.0}, return2D=False)

Calculates P(M500) assuming a y0 - M relation (default values assume UPP scaling relation from Arnaud et al. 2010), taking into account the steepness of the mass function. The approach followed is described in H13, Section 3.2. The binning for P(M500) is set according to the given mockSurvey, as are the assumed cosmological parameters.

This routine is used by calcM500Fromy0.

If return2D == True, returns a grid of same dimensions / binning as mockSurvey.z, mockSurvey.log10M, normalised such that the sum of the values is 1.

calcQ(theta500Arcmin, tck)

Returns Q, given theta500Arcmin, and a set of spline fit knots for (theta, Q).

calcR500Mpc(z, M500, cosmoModel)

Given z, M500 (in MSun), returns R500 in Mpc, with respect to critical density.

calcTheta500Arcmin(z, M500, cosmoModel)

Given z, M500 (in MSun), returns angular size equivalent to R500, with respect to critical density.

calcWeightedFRel(z, M500, fRelWeightsDict)

Return fRel for the given (z, M500), weighted by frequency according to fRelWeightsDict

convertM200mToM500c(M200m, z)

Returns M500c (MSun), R500c (Mpc) for given M200m and redshift. Uses the Bhattacharya et al. c-M relation: http://adsabs.harvard.edu/abs/2013ApJ…766…32B

See also Hu & Kravtsov: http://iopscience.iop.org/article/10.1086/345846/pdf

convertM500cToM200m(M500c, z)

Returns M200m given M500c

criticalDensity(z)

Returns the critical density at the given z.

fSZ(obsFrequencyGHz)

Returns the frequency dependence of the (non-relativistic) Sunyaev-Zel’dovich effect.

fitQ(config)

Calculates the filter mismatch function Q on a grid of scale sizes (given in terms of theta500 in arcmin), for each tile in the map. The results are initially cached (with a separate .fits table for each tile) under the selFn directory, before being combined into a single file at the end of a nemo run. Use GNFWParams (in config.parDict) if you want to specify a different cluster profile shape.

Parameters

config (startUp.NemoConfig) – A NemoConfig object.

Note

See loadQ for how to load in the output of this routine.

getFRelWeights(config)

Returns a dictionary of frequency weights used in relativistic corrections for each tile. This is cached in the selFn/ dir after the first time this routine is called.

getM500FromP(P, log10M, calcErrors=True)

Returns M500 as the maximum likelihood value from given P(log10M) distribution, together with 1-sigma error bars (M500, -M500Err, +M500 err).

gz(zIn, zMax=1000, dz=0.1)

Calculates linear growth factor at redshift z. Use Dz if you want normalised to D(z) = 1.0 at z = 0.

See http://www.astronomy.ohio-state.edu/~dhw/A873/notes8.pdf for some notes on this.

loadFRelWeights(fRelWeightsFileName)

Returns a dictionary of frequency weights used in relativistic corrections for each tile (stored in a .fits table, made by getFRelWeights).

loadQ(source, tileNames=None)

Load the filter mismatch function Q (see Hasselfield et al. 2013) as a dictionary of spline fits.

Parameters
  • source (nemo.startUp.NemoConfig or str) – Either the path to a .fits table (containing Q fits for all tiles - this is normally selFn/QFit.fits), or a nemo.startUp.NemoConfig object (from which the path and tiles to use will be inferred).

  • tileNames (optional, list) – A list of tiles for which the Q function spline fit coefficients will be extracted. If source is a nemo.startUp.NemoConfig object, this should be set to None.

Returns

A dictionary (with tilNames as keys), containing spline knots for the Q function for each tile. Q values can then be obtained by using these with scipy.interpolate.splev().

makeArnaudModelProfile(z, M500, GNFWParams='default', cosmoModel=None)

Given z, M500 (in MSun), returns dictionary containing Arnaud model profile (well, knots from spline fit, ‘tckP’ - assumes you want to interpolate onto an array with units of degrees) and parameters (particularly ‘y0’, ‘theta500Arcmin’).

Use GNFWParams to specify a different shape. If GNFWParams = ‘default’, then the default parameters as listed in gnfw.py are used, i.e.,

GNFWParams = {‘gamma’: 0.3081, ‘alpha’: 1.0510, ‘beta’: 5.4905, ‘tol’: 1e-7, ‘npts’: 100}

Otherwise, give a dictionary that specifies the wanted values. This would usually be specified as GNFWParams in the filter params in the nemo .par file (see the example .par files).

If cosmoModel is None, use default (Om0, Ol0, H0) = (0.3, 0.7, 70 km/s/Mpc) cosmology.

Used by ArnaudModelFilter

makeArnaudModelSignalMap(z, M500, degreesMap, wcs, beam, GNFWParams='default', amplitude=None, maxSizeDeg=15.0, convolveWithBeam=True)

Makes a 2d signal only map containing an Arnaud model cluster.

Parameters
  • z (float) – Redshift; used for setting angular size.

  • M500 (float) – Mass within R500, defined with respect to critical density; units are solar masses.

  • degreesMap (numpy.ndarray) – A 2d array containing radial distance measured in degrees from the centre of the model to be inserted. The output map will have the same dimensions and pixel scale (see nemoCython.makeDegreesDistanceMap).

  • GNFWParams (dict, optional) – Used to specify a different profile shape to the default (which follows Arnaud et al. 2010). If GNFWParams = ‘default’, then the default parameters as listed in gnfw.py are used, i.e., GNFWParams = {‘gamma’: 0.3081, ‘alpha’: 1.0510, ‘beta’: 5.4905, ‘tol’: 1e-7, ‘npts’: 100}. Otherwise, give a dictionary that specifies the wanted values. This would usually be specified using the GNFWParams key in the .yml config used when running nemo (see the examples/ directory).

  • amplitude (float, optional) – Amplitude of the cluster, i.e., the central decrement (in map units, e.g., uK), or the central Comptonization parameter (dimensionless), before beam convolution. Not needed for generating filter kernels.

  • maxSizeDeg (float, optional) – Use to limit the region over which the beam convolution is done, for optimization purposes.

  • convolveWithBeam (bool, optional) – If False, no beam convolution is done (it can be quicker to apply beam convolution over a whole source-injected map rather than per object).

Returns

signalMap (np.ndarray).

Note

The pixel window function is not applied here; use pixell.enmap.apply_window to do that (see nemo.filters.filterMaps).

makeBeamModelSignalMap(degreesMap, wcs, beam, amplitude=None)

Makes a 2d signal only map containing the given beam.

Parameters
  • degreesMap (np.ndarray) – Map of angular distance from the object position.

  • wcs (astWCS.WCS) – WCS corresponding to degreesMap.

  • beam (BeamProfile or str) – Either a BeamProfile object, or a string that gives the path to a text file that describes the beam profile.

  • ( (amplitude) – obj: float, optional): Specifies the amplitude of the input signal (in map units, e.g., uK), before beam convolution. This is only needed if this routine is being used to inject sources into maps. It is not needed for making filter kernels.

Returns

signalMap (np.ndarray)

Note

The pixel window function is not applied here; use pixell.enmap.apply_window to do that (see nemo.filters.filterMaps).

makeCombinedQTable(config)

Writes dictionary of tables (containing individual tile Q fits) as a single .fits table.

Returns combined Q astropy table object

meanDensity(z)

Returns the mean density at the given z.

y0FromLogM500(log10M500, z, tckQFit, tenToA0=4.95e-05, B0=0.08, Mpivot=300000000000000.0, sigma_int=0.2, cosmoModel=None, fRelWeightsDict={148.0: 1.0})

Predict y0~ given logM500 (in MSun) and redshift. Default scaling relation parameters are A10 (as in H13).

Use cosmoModel (astropy.cosmology object) to change/specify cosmological parameters.

fRelWeightsDict is used to account for the relativistic correction when y0~ has been constructed from multi-frequency maps. Weights should sum to 1.0; keys are observed frequency in GHz.

Returns y0~, theta500Arcmin, Q

12.12. startUp

This module contains basic set-up stuff (making directories, parsing config etc.) used by all the scripts in bin/ (nemo, nemoMass, nemoSelFn etc.).

class NemoConfig(configFileName, makeOutputDirs=True, setUpMaps=True, selFnDir=None, MPIEnabled=False, divideTilesByProcesses=True, verbose=True, strictMPIExceptions=True)

An object that keeps track of nemo’s configuration, maps, and output directories etc..

parDict

Dictionary containing the contents of the config file.

Type

dict

rootOutDir

Path to the directory where all output will be written.

Type

str

filteredMapsDir

Name of the directory where filtered maps will be written.

Type

str

diagnosticsDir

Path to the directory where miscellaneous diagnostic data (e.g., filter kernel plots) will be written.

Type

str

unfilteredMapsDictList

List of dictionaries corresponding to maps needed.

Type

list

tileNames

List of map tiles (extension names) to operate on.

Type

list

MPIEnabled

If True, use MPI to divide tileNames list among processes.

Type

bool

comm

Used by MPI.

Type

MPI.COMM_WORLD

rank

Used by MPI.

Type

int

size

Used by MPI.

Type

int

Creates an object that keeps track of nemo’s configuration, maps, output directories etc..

Parameters
  • configFileName (str) – Path to a nemo .yml configuration file.

  • makeOutputDirs (bool, optional) – If True, create output directories (where maps, catalogs are stored).

  • setUpMaps (bool, optional) – If True, set-up data structures for handling maps (inc. breaking into tiles if wanted).

  • selFnDir (str, optional) – Path to the selFn directory (use to override the default location).

  • MPIEnabled (bool, optional) – If True, use MPI to divide the map into tiles, distributed among processes. This requires tileDefinitions and tileNoiseRegions to be given in the .yml config file.

  • divideTilesByProcesses (bool, optional) – If True, divides up the list of tiles optimally among the available MPI processes.

  • strictMPIExceptions (bool) – If True, MPI will abort if an Exception is encountered (the downside is that you will not get the full traceback, but at least you will not waste CPU cycles). If False, MPI probably will not abort if an Exception is encountered, but you will get the full traceback (the downside is your MPI job may never complete). These options are a compromise due to how mpi4py handles MPI errors (the default handling for mpi4py corresponds to strictMPIExceptions = False).

  • verbose (bool) – If True, print some info to the terminal while we set-up the config file.

restoreConfig()

Restores the parameters dictionary (self.parDict) to the original state specified in the config .yml file.

parseConfigFile(parDictFileName)

Parse a nemo .yml config file.

Parameters

parDictFileName (str) – Path to a nemo .yml configuration file.

Returns

A dictionary of parameters.