diff --git a/bin/zCluster b/bin/zCluster index bcf7a7d..284b119 100644 --- a/bin/zCluster +++ b/bin/zCluster @@ -305,25 +305,25 @@ def runOnCatalog(catalog, retriever, retrieverOptions, photoRedshiftEngine, outD gobj['r-z']=gobj['r']-gobj['z'] wantedKeys=['id', 'RADeg', 'decDeg', 'u', 'uErr', 'g', 'gErr', 'r', 'rErr', 'i', - 'iErr', 'z', 'zErr', 'Ks', 'KsErr', 'zPhot', 'odds', 'r-i', 'r-z'] - if os.path.exists(objOutDir+os.path.sep+"galaxyCatalog_%s.fits" % (obj['name'].replace(" ", "_"))) == True: - os.remove(objOutDir+os.path.sep+"galaxyCatalog_%s.fits" % (obj['name'].replace(" ", "_"))) + 'iErr', 'z', 'zErr', 'Ks', 'KsErr', 'w1', 'w1Err', 'w2', 'w2Err', 'r-i', 'r-z', 'zPhot', 'odds'] tab=atpy.Table() for key in wantedKeys: arr=[] + allSentinelVals=True for gobj in galaxyCatalog: try: arr.append(gobj[key]) + allSentinelVals=False except: arr.append(99) - tab.add_column(atpy.Column(np.array(arr), key)) + if allSentinelVals == False: + tab[key]=arr # NOTE: to cut down on disk space this takes, include only galaxies within some radius # Chosen one is just over 1.5 Mpc at z = 0.1 - tab.add_column(atpy.Column(astCoords.calcAngSepDeg(obj['RADeg'], obj['decDeg'], - np.array(tab['RADeg']), np.array(tab['decDeg'])), "rDeg")) + tab['rDeg']=astCoords.calcAngSepDeg(obj['RADeg'], obj['decDeg'], tab['RADeg'].data, tab['decDeg'].data) #tab=tab[np.where(tab['rDeg'] < 14./60.0)] tab.table_name="zCluster" - tab.write(objOutDir+os.path.sep+"galaxyCatalog_%s.fits" % (obj['name'].replace(" ", "_"))) + tab.write(objOutDir+os.path.sep+"galaxyCatalog_%s.fits" % (obj['name'].replace(" ", "_")), overwrite = True) try: catalogs.catalog2DS9(galaxyCatalog, objOutDir+os.path.sep+"galaxyCatalog_%s.reg" % (obj['name'].replace(" ", "_")), @@ -637,6 +637,8 @@ if __name__ == '__main__': if cacheDir is not None: retrieverOptions['altCacheDir']=cacheDir + + retrieverOptions['fetchAndCacheOnly']=args.fetchAndCacheOnly photoRedshiftEngine=PhotoRedshiftEngine.PhotoRedshiftEngine(absMagCut, passbandSet = passbandSet, ZPError = ZPError, ZPOffsets = ZPOffsets, diff --git a/bin/zField b/bin/zField index 817d52b..835441e 100644 --- a/bin/zField +++ b/bin/zField @@ -16,6 +16,7 @@ from zCluster import retrievers from zCluster import PhotoRedshiftEngine from zCluster import clusters from zCluster import catalogs +import pickle import urllib import time @@ -77,6 +78,10 @@ def makeParser(): masses will be estimated at the given redshift. Otherwise, the stellar masses will be\ estimated at the maximum likelihood redshift of each galaxy (although that isn't\ actually implemented yet).") + parser.add_argument("-p", "--pickle-galaxy-catalog", dest="pickleGalaxyCatalog", default = False, + action = "store_true", help = "If given, pickle the galaxy catalog after estimating\ + photo-zs. Use this if you want to e.g. read p(z) for each galaxy from disk for some\ + other analysis.") return parser @@ -162,6 +167,16 @@ if __name__ == '__main__': if args.stellarMassModelDir is not None: photoRedshiftEngine.estimateStellarMasses(galaxyCatalog, args.stellarMassModelDir, z = z) + # Dump pickle file for p(z), optional + if args.pickleGalaxyCatalog is True: + pickleFileName=args.outDir+os.path.sep+"galaxyCatalog.pkl" + with open(pickleFileName, "wb") as pickleFile: + pickler=pickle.Pickler(pickleFile) + pickler.dump(galaxyCatalog) + # with open(pickleFileName, "rb") as pickleFile: + # unpickler=pickle.Unpickler(pickleFile) + # galaxyCatalog=unpickler.load() + # Write catalog and region file wantedKeys=['id', 'RADeg', 'decDeg', 'zPhot', 'odds'] if args.stellarMassModelDir is not None: diff --git a/setup.cfg b/setup.cfg index e90d31f..a834d83 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,21 @@ +[metadata] +name = zCluster +author = Matt Hilton + zCluster Contributors +author_email = matt.hilton@wits.ac.za +url = https://github.com/ACTCollaboration/zCluster +description = A code for measuring galaxy cluster photometric redshifts +long_description = file: README.rst +keywords = astrophysics +license = GPL v3 +classifiers = + Programming Language :: Python :: 3 + Environment :: Console + Intended Audience :: Science/Research + License :: OSI Approved :: GNU General Public License v3 (GPLv3) + Natural Language :: English + Operating System :: POSIX + Topic :: Scientific/Engineering :: Astronomy + [versioneer] VCS = git style = pep440 @@ -6,3 +24,11 @@ versionfile_build = zCluster/_version.py tag_prefix = v parentdir_prefix = zCluster- +[options] +install_requires = + astropy >= 4.0 + numpy >= 1.19 + matplotlib >= 2.0 + astLib >= 0.11.10 + scipy >= 1.0 + requests diff --git a/setup.py b/setup.py index 9f12dff..c60799a 100644 --- a/setup.py +++ b/setup.py @@ -4,35 +4,14 @@ import glob from setuptools import setup from setuptools import Extension -from setuptools.command.install import install -import stat -#from Cython.Distutils import build_ext -import numpy import versioneer -cmdclass=versioneer.get_cmdclass() -#cmdclass['build_ext']=build_ext - setup(name='zCluster', version=versioneer.get_version(), - cmdclass=cmdclass, - url="https://github.com/ACTCollaboration/zCluster", author='Matt Hilton + zCluster contributors', - author_email='matt.hilton@mykolab.com', - classifiers=[], - description='A code for measuring galaxy cluster photometric redshifts.', - long_description="""A code for measuring galaxy cluster photometric redshifts. Runs on both large scale - public survey data (e.g., SDSS) and user-supplied photometric catalogs.""", + author_email='matt.hilton@wits.ac.za', packages=['zCluster'], package_data={'zCluster': ['data/*', 'SED/CWW/*', 'SED/BR07/*', 'SED/EAZY_v1.0/*', 'passbands/*']}, scripts=['bin/zCluster', 'bin/zField', 'bin/zClusterBCG', 'bin/zClusterComparisonPlot'], - #ext_modules=[Extension("zClusterCython", ["zCluster/zClusterCython.pyx"], include_dirs=[numpy.get_include()])], - install_requires=["astropy >= 4.0", - "numpy >= 1.19", - "matplotlib >= 2.0", - "astLib >= 0.11.7", - "scipy >= 1.0", - #"cython", - "requests"] ) diff --git a/zCluster/PhotoRedshiftEngine.py b/zCluster/PhotoRedshiftEngine.py index 5d21ebf..12448df 100644 --- a/zCluster/PhotoRedshiftEngine.py +++ b/zCluster/PhotoRedshiftEngine.py @@ -24,7 +24,12 @@ class PhotoRedshiftEngine: """ def __init__(self, absMagCut, passbandSet = 'SDSS+Ks', zMin = 0.01, zMax = 3.0, zStep = 0.01, - ZPError = 0.0, ZPOffsets = None, templatesDir = None): + ZPError = 0.0, ZPOffsets = None, templatesDir = None, EBMinusVList = [0.0], + emLinesScaleList = [0.0]): + # EBMinusVList = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5], + # emLinesScaleList = [0.0, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0]): + # EBMinusVList = list(np.linspace(0, 0.48, 13))): + # EBMinusVList = list(np.linspace(0, 1.5, 16)),): """Sets up the stuff we would otherwise calculate every time, i.e., the templates. """ @@ -133,7 +138,7 @@ def __init__(self, absMagCut, passbandSet = 'SDSS+Ks', zMin = 0.01, zMax = 3.0, self.templateIndex=unpickler.load() self.modelFlux2=self.modelFlux**2 else: - self.numModels=len(self.SEDFiles) + self.numModels=len(self.SEDFiles)*len(EBMinusVList)*len(emLinesScaleList) i=0 t=0 self.templateIndex=[] @@ -141,17 +146,23 @@ def __init__(self, absMagCut, passbandSet = 'SDSS+Ks', zMin = 0.01, zMax = 3.0, s=astSED.SED() s.loadFromFile(f) t=t+1 - for z in self.zRange: - s.redshift(z) - modelSEDDict=s.getSEDDict(self.passbandsList) - modelSEDDict['E(B-V)']=None - modelSEDDict['ageGyr']=0.0 - modelSEDDict['z']=z - modelSEDDict['fileName']=f - modelSEDDict['modelListIndex']=i - modelSEDDict['SED']=s.copy() - self.modelSEDDictList.append(modelSEDDict) - self.templateIndex.append(t) + for emScaleFactor in emLinesScaleList: + if emScaleFactor > 0: + s.addEmissionLines(emScaleFactor) + for EBMinusV in EBMinusVList: + if EBMinusV > 0: + s.extinctionCalzetti(EBMinusV) # modifies z0flux, so this should be okay + for z in self.zRange: + s.redshift(z) + modelSEDDict=s.getSEDDict(self.passbandsList) + modelSEDDict['E(B-V)']=EBMinusV + modelSEDDict['ageGyr']=0.0 + modelSEDDict['z']=z + modelSEDDict['fileName']=f + modelSEDDict['modelListIndex']=i + # modelSEDDict['SED']=s.copy() # Only add this if we need [e.g. rest frame colours] - uses lots of RAM + self.modelSEDDictList.append(modelSEDDict) + self.templateIndex.append(t) i=i+1 del s self.templateIndex=np.array(self.templateIndex) diff --git a/zCluster/retrievers.py b/zCluster/retrievers.py index 3012701..03210b6 100644 --- a/zCluster/retrievers.py +++ b/zCluster/retrievers.py @@ -19,6 +19,10 @@ import time import zCluster import requests +try: + from dl import queryClient as qc +except: + print("WARNING: Failed to import dl module - retrievers that use NOAO Data Lab will not work.") from astropy.io.votable import parse_single_table #------------------------------------------------------------------------------------------------------------- @@ -111,7 +115,11 @@ def getRetriever(database, maxMagError = 0.2): elif database == 'DELVEDR2': retriever=DELVEDR2Retriever retrieverOptions={'maxMagError': maxMagError} - elif database.find("DECaLS") != -1: + elif database == 'DL_DECaLSDR10': + retriever=DL_DECaLSDR10Retriever + passbandSet='DECaLS' + retrieverOptions={'maxMagError': maxMagError} + elif database.find("DECaLS") != -1 and database != 'DECaLSDR10DL': if database == "DECaLS": raise Exception("Specify either 'DECaLSDR8', 'DECaLSDR9', or 'DECaLSDR10' instead of DECaLS") if database == 'DECaLSDR8': @@ -924,7 +932,7 @@ def SDSSRetriever(RADeg, decDeg, halfBoxSizeDeg = 18.0/60.0, DR = 7, optionsDict # For some reason, SDSS changed their whole web API in ~May 2016 without calling it a new DR #url='http://skyserver.sdss.org/dr12/en/tools/search/x_sql.aspx' #url='http://skyserver.sdss.org/dr12/en/tools/search/x_results.aspx' - url='http://skyserver.sdss.org/dr12/en/tools/search/x_results.aspx?searchtool=SQL&TaskName=Skyserver.Search.SQL&syntax=NoSyntax&ReturnHtml=false&' + url='https://skyserver.sdss.org/dr12/en/tools/search/x_results.aspx?searchtool=SQL&TaskName=Skyserver.Search.SQL&syntax=NoSyntax&ReturnHtml=false&' outFileName=cacheDir+os.path.sep+"SDSSDR12_%.4f_%.4f_%.4f.csv" % (RADeg, decDeg, halfBoxSizeDeg) lineSkip=2 @@ -988,7 +996,7 @@ def SDSSRetriever(RADeg, decDeg, halfBoxSizeDeg = 18.0/60.0, DR = 7, optionsDict try: response=urllib.request.urlopen(url+'%s' % (params)) except: - print("Network down? Waiting 30 sec...") + print("Network down? Waiting 30 sec... - if this persists, probably the query URL has changed.") time.sleep(30) # Some faffing about here because of python2 -> python3 @@ -1261,6 +1269,76 @@ def DECaLSDR10Retriever(RADeg, decDeg, halfBoxSizeDeg = 36.0/60.0, DR = None, op return stuff +#------------------------------------------------------------------------------------------------------------- +def DL_DECaLSDR10Retriever(RADeg, decDeg, halfBoxSizeDeg = 36.0/60.0, DR = None, optionsDict = {}): + """DECaLS DR10 retriever, using NOAO datalab. + + """ + + makeCacheDir() + if 'altCacheDir' in list(optionsDict.keys()): + cacheDir=optionsDict['altCacheDir'] + else: + cacheDir=CACHE_DIR + if os.path.exists(cacheDir) == False: + os.makedirs(cacheDir, exist_ok = True) + + outFileName=cacheDir+os.path.sep+"DL_DECaLSDR10_%.4f_%.4f_%.2f.fits" % (RADeg, decDeg, halfBoxSizeDeg) + if os.path.exists(outFileName) == False: + RAMin, RAMax, decMin, decMax=astCoords.calcRADecSearchBox(RADeg, decDeg, halfBoxSizeDeg) + try: + result=qc.query(sql='select objid, ra, dec, dered_mag_g, dered_mag_r, dered_mag_i, dered_mag_z, dered_mag_w1, dered_mag_w2,\ + snr_g, snr_r, snr_i, snr_z, snr_w1, snr_w2, type, maskbits from ls_dr10.tractor where\ + RA BETWEEN %.6f AND %.6f AND DEC BETWEEN %.6f and %.6f' % (RAMin, RAMax, decMin, decMax), + fmt = 'table') + result.write(outFileName, overwrite = True) + except: + result=None + print("... WARNING: datalab query failed to get %s" % (outFileName)) + else: + if 'fetchAndCacheOnly' in optionsDict.keys() and optionsDict['fetchAndCacheOnly'] == True: + print("... already retrieved: %s ..." % (outFileName)) + return None + print("... reading from cache: %s ..." % (outFileName)) + result=atpy.Table().read(outFileName) + + if result is None: + return None + + # DECaLS redshifts go very wrong when there are stars bright in W1, W2 in the vicinity + # This should fix - we'll also throw out PSF-shaped sources too + tab=result + if len(tab) == 0: + return None + tab=tab[tab['maskbits'] != 2**1] + tab=tab[tab['maskbits'] < 2**11] + tab=tab[np.where(tab['type'] != 'PSF')] + tab=tab[np.where(tab['type'] != 'PSF ')] # Trailing space - probably not an issue on datalab + + # WISE fluxes are available... i-band added in DECaLS DR10 + bands=['g', 'r', 'i', 'z', "w1", "w2"] + catalog=[] + for row in tab: + photDict={} + photDict['id']=row['objid'] + photDict['RADeg']=row['ra'] + photDict['decDeg']=row['dec'] + for b in bands: + if row['snr_%s' % (b)] > 0: + photDict[b]=row['dered_mag_%s' % (b)] + photDict[b+"Err"]=1/row['snr_%s' % (b)] + else: + photDict[b]=99.0 + photDict[b+'Err']=99.0 + if 'maxMagError' in list(optionsDict.keys()): + keep=checkMagErrors(photDict, optionsDict['maxMagError'], bands = bands) + else: + keep=True + if keep == True: + catalog.append(photDict) + + return catalog + #------------------------------------------------------------------------------------------------------------- def DELVEDR2Retriever(RADeg, decDeg, halfBoxSizeDeg = 36.0/60.0, DR = None, optionsDict = {}): """DELVE DR2 retriever, using NOAO datalab.