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/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/completeness.py b/nemo/completeness.py index 530fbb0..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: @@ -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 6672f2c..282a064 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,11 +616,16 @@ 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: + 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)) - 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..b022425 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_SNMap.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, 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: