From 5e3e7533d7cfd3c0b82589e429a249e343c041ce Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Fri, 10 Jul 2020 17:52:56 +0200 Subject: [PATCH 1/6] New smart (binned) option for noise maps This works... doesn't grid, just uses areas of weight maps with similar levels, and estimates noise in filtered maps in those regions. Added example config file. --- examples/pointSources/PS_S18d_f220_202006.sh | 9 + examples/pointSources/PS_S18d_f220_202006.yml | 177 ++++++++++++++++++ nemo/filters.py | 129 +++++++++---- 3 files changed, 278 insertions(+), 37 deletions(-) create mode 100644 examples/pointSources/PS_S18d_f220_202006.sh create mode 100644 examples/pointSources/PS_S18d_f220_202006.yml diff --git a/examples/pointSources/PS_S18d_f220_202006.sh b/examples/pointSources/PS_S18d_f220_202006.sh new file mode 100644 index 0000000..288c181 --- /dev/null +++ b/examples/pointSources/PS_S18d_f220_202006.sh @@ -0,0 +1,9 @@ +#!/bin/sh +#SBATCH --nodes=7 +#SBATCH --ntasks-per-node=20 +#SBATCH --mem=64000 +#SBATCH --time=23:59:00 + +source ~/.bashrc +time mpiexec nemo PS_S18d_f220_202006.yml -M -n + diff --git a/examples/pointSources/PS_S18d_f220_202006.yml b/examples/pointSources/PS_S18d_f220_202006.yml new file mode 100644 index 0000000..c9058d3 --- /dev/null +++ b/examples/pointSources/PS_S18d_f220_202006.yml @@ -0,0 +1,177 @@ +# Nemo config file +# YAML format +# - use null to return None in Python +# - note that YAML is fussy about large numbers: use e.g. 1.0e+14 for M500MSun (not 1e14) + +# Valid units are uK or Jy/sr +# this should be a list of maps at different frequencies +# NOTE: surveyMask is optional +unfilteredMaps: + - {mapFileName: "maps/Jun2020/act_s08_s18_cmb_f220_daynight_map.fits", + weightsFileName: "maps/Jun2020/act_s08_s18_cmb_f220_daynight_ivar.fits", + obsFreqGHz: 220.0, units: 'uK', + beamFileName: "maps/Jun2020/beams/s17_pa4_f220_nohwp_night_beam_profile_jitter.txt"} +# - {mapFileName: "maps/Mar2020/act_s08_s18_cmb_f090_daynight_map.fits", +# weightsFileName: "maps/Mar2020/act_s08_s18_cmb_f090_daynight_ivar.fits", +# obsFreqGHz: 98.3, units: 'uK', +# beamFileName: "Beams/190809/b20190809_s16_pa3_f090_nohwp_night_beam_profile_jitter_cmb.txt"} + +# Masks +surveyMask: "maps/Jun2020/AdvACTSurveyMask_v7_S18.fits" +#surveyMask: 'maps/Sep2019/surveyMask_v7_S18_inc_extsrc.fits' + +# Instead of giving point source mask, can give catalog instead +#maskPointSourcesFromCatalog: +# - "PS_S18_f150_auto/PS_S18_f150_auto_optimalCatalog.fits" +# - "customPSMask_S18/customPSCatalog_S18.fits" + +# Detection/catalog options +# Set useInterpolator; True for sub-pixel flux and SNR measurements +thresholdSigma: 4.0 +minObjPix: 1 +findCenterOfMass: True +useInterpolator: True +rejectBorder: 0 +objIdent: 'ACT-S' +longNames: False + +# Photometry options +#photFilter: 'Arnaud_M2e14_z0p4' + +# Optionally override the GNFW parameters - if not present, Arnaud et al. (2010) parameters are used +# The example below is for the Planck Pressure Profile (PPP) +#GNFWParams: {P0: 6.41, c500: 1.81, gamma: 0.31, alpha: 1.33, beta: 4.13, tol: 1e-7, npts: 100} + +# Mass measurement options - used by nemoMass and nemoSelFn scripts +# Writes out .fits file to nemoOutDir/nemoOutDir_M500.fits +# redshiftCatalog: A .fits table containing name, RADeg, decDeg, redshift, redshiftErr columns +# forcedPhotometry: If True, calc mass based on extracted y0~ in 'photFilter' map at RADeg, decDeg as given in redshiftCatalog +# If False, cross match redshiftCatalog with optimal catalog made by nemo +# Q: If 'H13', use fit to Q from results presented in H13 +# If 'fit', use fit to (theta, Q) done by nemo for 'photFilter' kernel +# tenToA0, B0, Mpivot, sigma_int: Fixed scaling relation options (see H13 or ACTPol paper) +# rescaleFactor, rescaleFactorErr: For MCal masses, as in the ACTPol paper (i.e., just rescales M500 results by 1/rescaleFactor) +#massOptions: {tenToA0: 4.95e-5, +# B0: 0.08, +# Mpivot: 3.0e+14, +# sigma_int: 0.2, +# relativisticCorrection: True, +# rescaleFactor: 0.68, +# rescaleFactorErr: 0.11, +# redshiftCatalog: "AdvACT_redshifts.fits", +# forcedPhotometry: False, +# Q: 'fit'} + +# Selection function options +# NOTE: could eventually add 'completenessFraction' to 'massLimitMaps', which is why that's a dictionary list +# Use selFnFootprints to calculate average completeness in given sky areas - e.g., overlap with optical surveys +#calcSelFn: True +#selFnOptions: {fixedSNRCut: 5.0, +# massLimitMaps: [{z: 0.5}]} + +#selFnFootprints: +# - {label: "HSC", +# maskList: ["HSCCheckAndSelFn/s19a_fdfc_CAR_contarea_ziy-gt-5.fits"]} +# - {label: "KiDS", +# maskList: ["KiDSSelFn/mask_KiDSN.fits", "KiDSSelFn/mask_KiDSS.fits"]} +# - {label: "DES", +# maskList: ["DESY3/AdvACT_y3a2_footprint_griz_1exp_v2.0.fits"]} + +# Filter definitions: +# allFilters is a dictionary of parameters that will be copied into all mapFilters +# (these can be overridden by keys with the same name in mapFilters) +#allFilters: {class: "ArnaudModelMatchedFilter", +# params: {noiseParams: {method: "dataMap", +# noiseGridArcmin: 40.}, +# saveFilteredMaps: False, +# saveRMSMap: False, +# savePlots: False, +# saveDS9Regions: False, +# outputUnits: 'yc', +# edgeTrimArcmin: 0.0} +# } + +# mapFilters is a list of all the different filters to apply +# (keys in mapFilters with the same name as those in allFilters take priority) +#mapFilters: +# - {label: "Arnaud_M1e14_z0p2", +# params: {M500MSun: 1.0e+14, z: 0.2}} +# - {label: "Arnaud_M2e14_z0p2", +# params: {M500MSun: 2.0e+14, z: 0.2}} +# - {label: "Arnaud_M4e14_z0p2", +# params: {M500MSun: 4.0e+14, z: 0.2}} +# - {label: "Arnaud_M8e14_z0p2", +# params: {M500MSun: 8.0e+14, z: 0.2}} +# - {label: "Arnaud_M1e14_z0p4", +# params: {M500MSun: 1.0e+14, z: 0.4}} +# - {label: "Arnaud_M2e14_z0p4", +# params: {M500MSun: 2.0e+14, z: 0.4, +# saveFilteredMaps: True, +# savePlots: True}} +# - {label: "Arnaud_M4e14_z0p4", +# params: {M500MSun: 4.0e+14, z: 0.4}} +# - {label: "Arnaud_M8e14_z0p4", +# params: {M500MSun: 8.0e+14, z: 0.4}} +# - {label: "Arnaud_M1e14_z0p8", +# params: {M500MSun: 1.0e+14, z: 0.8}} +# - {label: "Arnaud_M2e14_z0p8", +# params: {M500MSun: 2.0e+14, z: 0.8}} +# - {label: "Arnaud_M4e14_z0p8", +# params: {M500MSun: 4.0e+14, z: 0.8}} +# - {label: "Arnaud_M8e14_z0p8", +# params: {M500MSun: 8.0e+14, z: 0.8}} +# - {label: "Arnaud_M1e14_z1p2", +# params: {M500MSun: 1.0e+14, z: 1.2}} +# - {label: "Arnaud_M2e14_z1p2", +# params: {M500MSun: 2.0e+14, z: 1.2}} +# - {label: "Arnaud_M4e14_z1p2", +# params: {M500MSun: 4.0e+14, z: 1.2}} +# - {label: "Arnaud_M8e14_z1p2", +# params: {M500MSun: 8.0e+14, z: 1.2}} + +# Set this to True to generate a sky sim (with noise), run all the filters over it, and measure contamination +# Set numSkySims to number required - we need to average over many as results vary a fair bit +estimateContaminationFromSkySim: False +numSkySims: 10 + +# Set this to True to estimate contamination by running cluster finder over inverted maps +# This is sensitive to how well point source masking is done +estimateContaminationFromInvertedMaps: False + +# Run position recovery test +positionRecoveryTest: False +posRecIterations: 1 +posRecSourcesPerTile: 200 +posRecModels: + - {redshift: 0.8, M500: 2.0e+14} + - {redshift: 0.4, M500: 2.0e+14} + - {redshift: 0.2, M500: 2.0e+14} + - {redshift: 0.1, M500: 2.0e+14} + +# tileDir options - cut-up each map into smaller sections +makeTileDir: True +makeQuickLookMaps: True +tileOverlapDeg: 1.0 +tileDefLabel: 'auto' +tileDefinitions: {mask: 'maps/Jun2020/AdvACTSurveyMask_v7_S18.fits', + targetTileWidthDeg: 10.0, + targetTileHeightDeg: 5.0} + +# Filter definitions: +mapFilters: + - {label: "Beam", + class: "BeamMatchedFilter", + params: {noiseParams: {method: "model", + noiseGridArcmin: "smart", + numNoiseBins: 100}, + saveFilteredMaps: True, + outputUnits: 'uK', + edgeTrimArcmin: 0.0}} + +# If this is given, only the named tiles will be processed (useful for testing) +#tileNameList: + #- '1_10_7' # powerful f150 source; do as set - sensitive to point-source mask threshold + #- '1_10_8' # J2327 (next to a source); do as set - sensitive to point-source mask threshold + #- '2_0_7' # powerful f150 source + #- '2_2_8' # powerful f150 source + #- '3_0_1' # powerful f150 source diff --git a/nemo/filters.py b/nemo/filters.py index 6672f2c..1e116ea 100644 --- a/nemo/filters.py +++ b/nemo/filters.py @@ -308,39 +308,31 @@ def makeNoiseMap(self, mapData): """Estimate the noise map using local RMS measurements on a grid, over the whole filtered map. """ - - # Grid method - #print "... making SN map ..." - gridSize=int(round((self.params['noiseParams']['noiseGridArcmin']/60.)/self.wcs.getPixelSizeDeg())) - #gridSize=rIndex*3 - overlapPix=int(gridSize/2) - numXChunks=mapData.shape[1]/gridSize - numYChunks=mapData.shape[0]/gridSize - yChunks=np.linspace(0, mapData.shape[0], int(numYChunks+1), dtype = int) - xChunks=np.linspace(0, mapData.shape[1], int(numXChunks+1), dtype = int) - #SNMap=np.zeros(mapData.shape) - apodMask=np.not_equal(mapData, 0) - # We could make below behaviour default if match photFilter? Would need to see photFilter though... - #if 'saveRMSMap' in self.params['noiseParams'] and self.params['noiseParams']['saveRMSMap'] == True: - RMSMap=np.zeros(mapData.shape) - for i in range(len(yChunks)-1): - for k in range(len(xChunks)-1): - y0=yChunks[i]-overlapPix - y1=yChunks[i+1]+overlapPix - x0=xChunks[k]-overlapPix - x1=xChunks[k+1]+overlapPix - if y0 < 0: - y0=0 - if y1 > mapData.shape[0]: - y1=mapData.shape[0] - if x0 < 0: - x0=0 - if x1 > mapData.shape[1]: - x1=mapData.shape[1] - chunkValues=mapData[y0:y1, x0:x1] - goodAreaMask=np.greater_equal(apodMask[y0:y1, x0:x1], 1.0) - + # 'smart option' - measure noise in areas with similar weights (however weights defined) + if self.params['noiseParams']['noiseGridArcmin'] == "smart": + # average all weight maps + # doesn't matter about spectral weighting, we just want areas with similar characteristics + medWeights=[] + for mapDict in self.unfilteredMapsDictList: + medWeights.append(mapDict['weights']) + medWeights=np.median(np.array(medWeights), axis = 0) + # Binning by weights - should make numBins adjustable + try: + numBins=self.params['noiseParams']['numNoiseBins'] + except: + raise Exception("Need to give numNoiseBins in noiseParams when using noiseGridArcmin = 'smart'") + binEdges=np.linspace(medWeights.min(), medWeights.max(), numBins) + RMSMap=np.zeros(medWeights.shape) + apodMask=np.not_equal(mapData, 0) + for i in range(len(binEdges)-1): + # Find area of similar weight + binMin=binEdges[i] + binMax=binEdges[i+1] + weightMask=np.logical_and(medWeights > binMin, medWeights < binMax) + # Measure noise in that area from filtered map, and fill in noise map + chunkValues=mapData[weightMask] + goodAreaMask=np.greater_equal(apodMask[weightMask], 1.0) if 'RMSEstimator' in self.params['noiseParams'].keys() and self.params['noiseParams']['RMSEstimator'] == 'biweight': if goodAreaMask.sum() >= 10: # Astropy version is faster but gives identical results @@ -353,7 +345,7 @@ def makeNoiseMap(self, mapData): else: # Default: 3-sigma clipped stdev if np.not_equal(chunkValues, 0).sum() != 0: - goodAreaMask=np.greater_equal(apodMask[y0:y1, x0:x1], 1.0) + goodAreaMask=np.greater_equal(apodMask[weightMask], 1.0) chunkMean=np.mean(chunkValues[goodAreaMask]) chunkRMS=np.std(chunkValues[goodAreaMask]) sigmaClip=3.0 @@ -365,9 +357,69 @@ def makeNoiseMap(self, mapData): chunkRMS=np.std(chunkValues[mask]) else: chunkRMS=0. - if chunkRMS > 0: - RMSMap[y0:y1, x0:x1]=chunkRMS + RMSMap[weightMask]=chunkRMS + + else: + # Grid method + #print "... making SN map ..." + gridSize=int(round((self.params['noiseParams']['noiseGridArcmin']/60.)/self.wcs.getPixelSizeDeg())) + #gridSize=rIndex*3 + overlapPix=int(gridSize/2) + numXChunks=mapData.shape[1]/gridSize + numYChunks=mapData.shape[0]/gridSize + yChunks=np.linspace(0, mapData.shape[0], int(numYChunks+1), dtype = int) + xChunks=np.linspace(0, mapData.shape[1], int(numXChunks+1), dtype = int) + #SNMap=np.zeros(mapData.shape) + apodMask=np.not_equal(mapData, 0) + # We could make below behaviour default if match photFilter? Would need to see photFilter though... + #if 'saveRMSMap' in self.params['noiseParams'] and self.params['noiseParams']['saveRMSMap'] == True: + RMSMap=np.zeros(mapData.shape) + for i in range(len(yChunks)-1): + for k in range(len(xChunks)-1): + y0=yChunks[i]-overlapPix + y1=yChunks[i+1]+overlapPix + x0=xChunks[k]-overlapPix + x1=xChunks[k+1]+overlapPix + if y0 < 0: + y0=0 + if y1 > mapData.shape[0]: + y1=mapData.shape[0] + if x0 < 0: + x0=0 + if x1 > mapData.shape[1]: + x1=mapData.shape[1] + chunkValues=mapData[y0:y1, x0:x1] + + goodAreaMask=np.greater_equal(apodMask[y0:y1, x0:x1], 1.0) + + if 'RMSEstimator' in self.params['noiseParams'].keys() and self.params['noiseParams']['RMSEstimator'] == 'biweight': + if goodAreaMask.sum() >= 10: + # Astropy version is faster but gives identical results + chunkRMS=apyStats.biweight_scale(chunkValues[goodAreaMask], c = 9.0, modify_sample_size = True) + #chunkRMS=astStats.biweightScale(chunkValues[goodAreaMask], 6.0) + else: + chunkRMS=0. + elif 'RMSEstimator' in self.params['noiseParams'].keys() and self.params['noiseParams']['RMSEstimator'] == 'percentile': + chunkRMS=np.percentile(abs(chunkValues[goodAreaMask]), 68.3) + else: + # Default: 3-sigma clipped stdev + if np.not_equal(chunkValues, 0).sum() != 0: + goodAreaMask=np.greater_equal(apodMask[y0:y1, x0:x1], 1.0) + chunkMean=np.mean(chunkValues[goodAreaMask]) + chunkRMS=np.std(chunkValues[goodAreaMask]) + sigmaClip=3.0 + for c in range(10): + mask=np.less(abs(chunkValues), abs(chunkMean+sigmaClip*chunkRMS)) + mask=np.logical_and(goodAreaMask, mask) + if mask.sum() > 0: + chunkMean=np.mean(chunkValues[mask]) + chunkRMS=np.std(chunkValues[mask]) + else: + chunkRMS=0. + + if chunkRMS > 0: + RMSMap[y0:y1, x0:x1]=chunkRMS return RMSMap @@ -564,8 +616,11 @@ def buildAndApply(self): # NOTE: Now point source mask is applied above, we fill the holes back in here when finding edges if 'edgeTrimArcmin' in self.params.keys() and self.params['edgeTrimArcmin'] > 0: trimSizePix=int(round((self.params['edgeTrimArcmin']/60.)/self.wcs.getPixelSizeDeg())) - else: - gridSize=int(round((self.params['noiseParams']['noiseGridArcmin']/60.)/self.wcs.getPixelSizeDeg())) + elif 'noiseGridArcmin' in self.params['noiseParams']: + try: + gridSize=int(round((self.params['noiseParams']['noiseGridArcmin']/60.)/self.wcs.getPixelSizeDeg())) + except: + gridSize=int(round((40./60.)/self.wcs.getPixelSizeDeg())) trimSizePix=int(round(gridSize*3.0)) edgeCheck=ndimage.rank_filter(abs(filteredMap+(1-psMask)), 0, size = (trimSizePix, trimSizePix)) edgeCheck=np.array(np.greater(edgeCheck, 0), dtype = float) From 290283462ced4151293b9f60c7d13a50e95ca87b Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Sat, 11 Jul 2020 10:28:03 +0200 Subject: [PATCH 2/6] Map stitching properly implemented --- bin/nemo | 38 ++++++++++++++++++--------------- nemo/completeness.py | 6 +++--- nemo/filters.py | 15 +++++++------ nemo/maps.py | 51 +++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 82 insertions(+), 28 deletions(-) diff --git a/bin/nemo b/bin/nemo index 6d899f4..cc8c4b7 100644 --- a/bin/nemo +++ b/bin/nemo @@ -89,23 +89,6 @@ if __name__ == '__main__': addInfo = addInfo, color = "cyan") else: if config.rank == 0: print("... already made catalog %s ..." % (optimalCatalogFileName)) - - # Stitch together map tiles - these are 'quicklook' images (downsampled by factor 4 in resolution) - # The stitchTiles routine will only write output if there are multiple maps matching the file pattern - if config.rank == 0 and 'makeQuickLookMaps' in config.parDict.keys() and config.parDict['makeQuickLookMaps'] == True: - # RMS (we don't save quicklook to selFnDir - users should use full-fat maps to be sure) - if os.path.exists(config.diagnosticsDir+os.path.sep+"quicklook_RMSMap.fits") == False: - maps.stitchTiles(config.selFnDir+os.path.sep+"RMSMap*.fits", - config.diagnosticsDir+os.path.sep+"quicklook_RMSMap.fits", - config.quicklookWCS, config.quicklookShape, fluxRescale = config.quicklookScale) - else: - print("... already made %s ..." % (config.diagnosticsDir+os.path.sep+"quicklook_RMSMap.fits")) - # S/N maps at reference filter scale - quicklookSNMapPath=config.filteredMapsDir+os.path.sep+"quicklook_%s_SNMap.fits" % (config.parDict['photFilter']) - if config.parDict['photFilter'] is not None and os.path.exists(quicklookSNMapPath) == False: - maps.stitchTiles(config.filteredMapsDir+os.path.sep+"*"+os.path.sep+"%s*SNMap.fits" % (config.parDict['photFilter']), - quicklookSNMapPath, config.quicklookWCS, config.quicklookShape, - fluxRescale = config.quicklookScale) # Q function (filter mismatch) - if needed options have been given # We may as well do this here to save having to run nemoMass separately (though we still can...) @@ -185,6 +168,27 @@ if __name__ == '__main__': # Tidying up etc. if config.rank == 0: + + # Stitch together map tiles - full fat versions (this will only work 'saveFilteredMaps: True') + if 'stitchTiles' in config.parDict.keys() and config.parDict['stitchTiles'] == True: + maps.stitchTiles(config) + + # Stitch together map tiles - 'quicklook' images (downsampled by factor 4 in resolution) + # The stitchTilesQuickLook routine will only write output if there are multiple maps matching the file pattern + if 'makeQuickLookMaps' in config.parDict.keys() and config.parDict['makeQuickLookMaps'] == True: + # RMS (we don't save quicklook to selFnDir - users should use full-fat maps to be sure) + if os.path.exists(config.diagnosticsDir+os.path.sep+"quicklook_RMSMap.fits") == False: + maps.stitchTilesQuickLook(config.selFnDir+os.path.sep+"RMSMap*.fits", + config.diagnosticsDir+os.path.sep+"quicklook_RMSMap.fits", + config.quicklookWCS, config.quicklookShape, fluxRescale = config.quicklookScale) + else: + print("... already made %s ..." % (config.diagnosticsDir+os.path.sep+"quicklook_RMSMap.fits")) + # S/N maps at reference filter scale + quicklookSNMapPath=config.filteredMapsDir+os.path.sep+"quicklook_%s_SNMap.fits" % (config.parDict['photFilter']) + if config.parDict['photFilter'] is not None and os.path.exists(quicklookSNMapPath) == False: + maps.stitchTilesQuickLook(config.filteredMapsDir+os.path.sep+"*"+os.path.sep+"%s*SNMap.fits" % (config.parDict['photFilter']), + quicklookSNMapPath, config.quicklookWCS, config.quicklookShape, + fluxRescale = config.quicklookScale) # Plot tile-averaged position recovery test if sourceInjTable is not None: diff --git a/nemo/completeness.py b/nemo/completeness.py index 530fbb0..631045f 100644 --- a/nemo/completeness.py +++ b/nemo/completeness.py @@ -1119,9 +1119,9 @@ def makeFullSurveyMassLimitMapPlot(z, config): config.quicklookShape, config.quicklookWCS=maps.shrinkWCS(config.origShape, config.origWCS, config.quicklookScale) outFileName=config.diagnosticsDir+os.path.sep+"reproj_massLimitMap_z%s.fits" % (str(z).replace(".", "p")) - maps.stitchTiles(config.diagnosticsDir+os.path.sep+"*"+os.path.sep+"massLimitMap_z%s#*.fits" % (str(z).replace(".", "p")), - outFileName, config.quicklookWCS, config.quicklookShape, - fluxRescale = config.quicklookScale) + maps.stitchTilesQuickLook(config.diagnosticsDir+os.path.sep+"*"+os.path.sep+"massLimitMap_z%s#*.fits" % (str(z).replace(".", "p")), + outFileName, config.quicklookWCS, config.quicklookShape, + fluxRescale = config.quicklookScale) # Make plot if os.path.exists(outFileName) == True: diff --git a/nemo/filters.py b/nemo/filters.py index 1e116ea..e3c4eb4 100644 --- a/nemo/filters.py +++ b/nemo/filters.py @@ -616,14 +616,15 @@ def buildAndApply(self): # NOTE: Now point source mask is applied above, we fill the holes back in here when finding edges if 'edgeTrimArcmin' in self.params.keys() and self.params['edgeTrimArcmin'] > 0: trimSizePix=int(round((self.params['edgeTrimArcmin']/60.)/self.wcs.getPixelSizeDeg())) - elif 'noiseGridArcmin' in self.params['noiseParams']: - try: - gridSize=int(round((self.params['noiseParams']['noiseGridArcmin']/60.)/self.wcs.getPixelSizeDeg())) - except: - gridSize=int(round((40./60.)/self.wcs.getPixelSizeDeg())) + elif 'noiseGridArcmin' in self.params['noiseParams'] and self.params['noiseParams']['noiseGridArcmin'] != "smart": trimSizePix=int(round(gridSize*3.0)) - edgeCheck=ndimage.rank_filter(abs(filteredMap+(1-psMask)), 0, size = (trimSizePix, trimSizePix)) - edgeCheck=np.array(np.greater(edgeCheck, 0), dtype = float) + else: + trimSizePix=0.0 + if trimSizePix > 0: + edgeCheck=ndimage.rank_filter(abs(filteredMap+(1-psMask)), 0, size = (trimSizePix, trimSizePix)) + edgeCheck=np.array(np.greater(edgeCheck, 0), dtype = float) + else: + edgeCheck=np.ones(filteredMap.shape) filteredMap=filteredMap*edgeCheck apodMask=np.not_equal(filteredMap, 0) surveyMask=edgeCheck*surveyMask*psMask diff --git a/nemo/maps.py b/nemo/maps.py index 32be051..5b19aa1 100644 --- a/nemo/maps.py +++ b/nemo/maps.py @@ -522,9 +522,58 @@ def checkMask(fileName): if hdu.data is not None: if np.less(hdu.data, 0).sum() > 0: raise Exception("Mask file '%s' contains negative values - please fix your mask." % (fileName)) + +#------------------------------------------------------------------------------------------------------------- +def stitchTiles(config): + """Stitches together full size filtered maps, SN maps, area maps, and noise maps that have been previously + been saved as tiles. + """ + + # Defining the maps to stitch together and where they will go + stitchDictList=[{'pattern': config.filteredMapsDir+os.path.sep+"$TILENAME"+os.path.sep+"$FILTLABEL#$TILENAME_filteredMap.fits", + 'outFileName': config.filteredMapsDir+os.path.sep+"stitched_$FILTLABEL_filteredMap.fits", + 'compressed': False}, + {'pattern': config.filteredMapsDir+os.path.sep+"$TILENAME"+os.path.sep+"$FILTLABEL#$TILENAME_filteredMap.fits", + 'outFileName': config.filteredMapsDir+os.path.sep+"stitched_$FILTLABEL_SNMap.fits", + 'compressed': False}, + {'pattern': config.selFnDir+os.path.sep+"areaMask#$TILENAME.fits", + 'outFileName': config.selFnDir+os.path.sep+"stitched_areaMask.fits", + 'compressed': True}, + {'pattern': config.selFnDir+os.path.sep+"RMSMap_$FILTLABEL#$TILENAME.fits", + 'outFileName': config.selFnDir+os.path.sep+"stitched_RMSMap_$FILTLABEL.fits", + 'compressed': True}] + + tileCoordsDict=config.tileCoordsDict + for filterDict in config.parDict['mapFilters']: + if filterDict['params']['saveFilteredMaps'] == True: + for stitchDict in stitchDictList: + pattern=stitchDict['pattern'] + outFileName=stitchDict['outFileName'].replace("$FILTLABEL", filterDict['label']) + compressed=stitchDict['compressed'] + d=np.zeros([config.origWCS.header['NAXIS2'], config.origWCS.header['NAXIS1']]) + wcs=config.origWCS + for tileName in tileCoordsDict.keys(): + f=pattern.replace("$TILENAME", tileName).replace("$FILTLABEL", filterDict['label']) + if os.path.exists(f) == True: + # Handle compressed or otherwise + with pyfits.open(f) as img: + tileData=img[0].data + if tileData is None: + for extName in img: + tileData=img[extName].data + if tileData is not None: + break + assert tileData is not None + else: + continue + areaMask, areaWCS=completeness.loadAreaMask(tileName, config.selFnDir) + minX, maxX, minY, maxY=config.tileCoordsDict[tileName]['clippedSection'] + d[minY:maxY, minX:maxX]=d[minY:maxY, minX:maxX]+areaMask*tileData + saveFITS(outFileName, d, wcs, compressed = compressed) + #------------------------------------------------------------------------------------------------------------- -def stitchTiles(filePattern, outFileName, outWCS, outShape, fluxRescale = 1.0): +def stitchTilesQuickLook(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, From 648ad7adbd49012e68b77268f9bf0ee00d5c14fd Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Tue, 14 Jul 2020 09:52:18 +0200 Subject: [PATCH 3/6] Increased tolerance for RMSTab and area mask check To 0.003 deg^2 difference, from 0.001. Must be rounding somewhere. Some tiles with 'smart' method fail at the 0.0013 level of tolerance but this will never make a difference to anything. --- nemo/completeness.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nemo/completeness.py b/nemo/completeness.py index 631045f..6e70001 100644 --- a/nemo/completeness.py +++ b/nemo/completeness.py @@ -652,7 +652,7 @@ def getRMSTab(tileName, photFilterLabel, selFnDir, diagnosticsDir = None, footpr RMSTab.add_column(atpy.Column(tileArea, 'areaDeg2')) RMSTab.add_column(atpy.Column(RMSValues, 'y0RMS')) # Sanity checks - these should be impossible but we have seen (e.g., when messed up masks) - tol=1e-3 + tol=0.003 if abs(RMSTab['areaDeg2'].sum()-areaMapSqDeg.sum()) > tol: raise Exception("Mismatch between area map and area in RMSTab for tile '%s'" % (tileName)) if np.less(RMSTab['areaDeg2'], 0).sum() > 0: From f960f0fbbfadcbb8f81a19ffed76b9ae96a42fc3 Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Mon, 20 Jul 2020 16:37:09 +0200 Subject: [PATCH 4/6] Added missing gridSize --- nemo/filters.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nemo/filters.py b/nemo/filters.py index e3c4eb4..282a064 100644 --- a/nemo/filters.py +++ b/nemo/filters.py @@ -617,6 +617,7 @@ def buildAndApply(self): if 'edgeTrimArcmin' in self.params.keys() and self.params['edgeTrimArcmin'] > 0: trimSizePix=int(round((self.params['edgeTrimArcmin']/60.)/self.wcs.getPixelSizeDeg())) elif 'noiseGridArcmin' in self.params['noiseParams'] and self.params['noiseParams']['noiseGridArcmin'] != "smart": + gridSize=int(round((self.params['noiseParams']['noiseGridArcmin']/60.)/self.wcs.getPixelSizeDeg())) trimSizePix=int(round(gridSize*3.0)) else: trimSizePix=0.0 From 4389fb5c59edd4df2863cb55f1e26a255d4b6704 Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Tue, 21 Jul 2020 17:24:24 +0200 Subject: [PATCH 5/6] Fix for nemoMass forced photometry mode To stop writing tile area maps again when not necessary --- nemo/pipelines.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nemo/pipelines.py b/nemo/pipelines.py index c89f9b8..1afde60 100644 --- a/nemo/pipelines.py +++ b/nemo/pipelines.py @@ -144,9 +144,10 @@ def filterMapsAndMakeCatalogs(config, rootOutDir = None, copyFilters = False, me DS9RegionsPath = DS9RegionsPath) # We write area mask here, because it gets modified by findObjects if removing rings + # NOTE: condition added to stop writing tile maps again when running nemoMass in forced photometry mode maskFileName=config.selFnDir+os.path.sep+"areaMask#%s.fits" % (tileName) surveyMask=np.array(filteredMapDict['surveyMask'], dtype = int) - if os.path.exists(maskFileName) == False: + if os.path.exists(maskFileName) == False and os.path.exists(config.selFnDir+os.path.sep+"areaMask.fits") == False: maps.saveFITS(maskFileName, surveyMask, filteredMapDict['wcs'], compressed = True) if measureFluxes == True: From b3bd5f9b840bee36b8c7a32e1df5245dfd48a232 Mon Sep 17 00:00:00 2001 From: Matt Hilton Date: Wed, 22 Jul 2020 06:37:51 +0200 Subject: [PATCH 6/6] Fix for dumb pattern matching bug in map stitch --- nemo/maps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nemo/maps.py b/nemo/maps.py index 5b19aa1..b022425 100644 --- a/nemo/maps.py +++ b/nemo/maps.py @@ -534,7 +534,7 @@ def stitchTiles(config): stitchDictList=[{'pattern': config.filteredMapsDir+os.path.sep+"$TILENAME"+os.path.sep+"$FILTLABEL#$TILENAME_filteredMap.fits", 'outFileName': config.filteredMapsDir+os.path.sep+"stitched_$FILTLABEL_filteredMap.fits", 'compressed': False}, - {'pattern': config.filteredMapsDir+os.path.sep+"$TILENAME"+os.path.sep+"$FILTLABEL#$TILENAME_filteredMap.fits", + {'pattern': config.filteredMapsDir+os.path.sep+"$TILENAME"+os.path.sep+"$FILTLABEL#$TILENAME_SNMap.fits", 'outFileName': config.filteredMapsDir+os.path.sep+"stitched_$FILTLABEL_SNMap.fits", 'compressed': False}, {'pattern': config.selFnDir+os.path.sep+"areaMask#$TILENAME.fits",