# This script is used reconstruct the script of DBSCAN. # To make is more usefule to use. and can be imported septeratrely import numpy as np from astropy import units as u from spectral_cube import SpectralCube from astropy.table import Table, vstack # need to download myPYTHON.py from myPYTHON import * import os from progressbar import * import gc import sys from skimage.morphology import erosion, dilation from scipy.ndimage import label, generate_binary_structure, binary_erosion, binary_dilation from sklearn.cluster import DBSCAN import requests import configparser from astropy.table import Column doFITS = myFITS() """ Version 7.0 adding control mask """ #version2.0 #Add a colnames to distinguish clouds due to bad channels # in development # print doFITS.weighted_avg_and_std # Version 6.0 # record LBV range for each cloud class MWISPDBSCAN( object ): __author__ = "MWISP" __version__ = 7.0 algorithmName="DBSCAN" #processPath = None #"./" # by default save the product in the current path backupPath = "./DBSCANresults/" processPath = backupPath rawCOFITS = None rmsFITS = None extraTag="" rawLabelFITSName = None rawCatName = None cleanCatName = None cleanLabelFITSName = None cleanTBWithDis = None #used only by matchDisTB, to save the distance table with distance assigned rawCOFITSConf = "rawCOFITS" rmsFITSConf = "rmsFITS" dbscanPathConf = "dbscanPath" lvPathConf = "lvPath" maskPathConf = "maskPath" extraTagConf = "extraTag" rmsData = None # must has the same shape with Ny,and Nz of the data file averageRMS = None coreFITS = None maskFirst= False ################# DBSCAN parameters #### should be fixed cutoff_sigma = 2 minPts = 4 connectivity = 1 ###################### seting of catalog selection saveCore = False minVox = 16 minChannel = 3 hasBeam = 1 minPeakSigma = 5 ### ##### parameters for cloud int cloudIntPath = None cloudSpecPath =None doFore=False #produce Foreround localVindex = None # either 0 or max(Nx), depends on Galacticl longitude doPureInt=False #produce pureInt #### running ID runID= None #only for cloud IT doClean=True #### not always saveIntFITSname = None savePureIntFITSname = None saveMaskFITSname = None saveForegroundFITSname = None cropExtend=2 #in degree regionName = None intFolderName = "cloudInt" specFolderName = "specAvg" areaThresh=0.01 #deg square, configFile = "/home/qzyan/WORK/dataMWISP/project/distanceMC/BEEPII/dbscanRecords.ini" ### by default, local saving #this file keep record DBSCASN that have been processed, you can alter this file #configFile = "/home/qzyan/WORK/dataMWISP/project/distanceMC/BEEPII/regions.ini" ### by default, local saving, if the save name #is config = None writeNoise= False pixelSize=None velResolution=None saveSpec = False controlMaskFITS=None #and additional mask file saveCloudID=None saveCloudIndex=None #### in order to speed up cloudCubeSaveFolder = "cloudSubCubes" def __init__(self,regionName=None,configureFile=None ): # self.getModelTB() print("Running DBSCAN with version ", self.__version__) if configureFile is not None and os.path.isfile(configureFile): print("Configure file found!") self.configFile = configureFile if regionName is not None: self.getDBSCANcase(regionName ) pass #def initDBSCAN(self,rawCOFITS,averageRMS=None,rmsFITS=None,processPath =None): def setCloudIntParameters(self,cloudIntPath=None,doFore=False,doPureInt=False,localVindex=False): """ :param cloudIntPath: :param doFore: :param doPureInt: :return: """ self.cloudIntPath = cloudIntPath self.doFore = doFore self.doPureInt = doPureInt self.localVindex = localVindex def setControlMaskFITS(self,controlMaskFITS): self.controlMaskFITS=controlMaskFITS def setDBSCANParameters(self, cutoff_sigma=2, minPts=4, connectivity=1): self.cutoff_sigma = cutoff_sigma self.minPts = minPts self.connectivity = connectivity def setCatalogSelectionCriteria(self, minVox=None, minChannel=None, hasBeam=None, minPeakSigma=None): if minVox is not None: self.minVox = minVox if minChannel is not None: self.minChannel = minChannel if hasBeam is not None: self.hasBeam = hasBeam if minPeakSigma is not None: self.minPeakSigma = minPeakSigma def sumEdgeByCon1(self, extendMask): # 7 in total # return extendMask[1:-1, 1:-1, 1:-1] + extendMask[0:-2, 1:-1, 1:-1] + extendMask[2:, 1:-1, 1:-1] + extendMask[1:-1, 0: -2, 1:-1]+ extendMask[1:-1, 2:, 1:-1]+ extendMask[1:-1, 1:-1, 0:-2] + extendMask[1:-1, 1:-1, 2:] raw = extendMask[1:-1, 1:-1, 1:-1] leftShiftZ = extendMask[0:-2, 1:-1, 1:-1] rightShiftZ = extendMask[2:, 1:-1, 1:-1] leftShiftY = extendMask[1:-1, 0: -2, 1:-1] rightShiftY = extendMask[1:-1, 2:, 1:-1] leftShiftX = extendMask[1:-1, 1:-1, 0:-2] rightShiftX = extendMask[1:-1, 1:-1, 2:] sumAll = raw + leftShiftZ + rightShiftZ + leftShiftY + rightShiftY + leftShiftX + rightShiftX return sumAll def sumEdgeByCon2(self, extendMask): # 27 in total sumAll = extendMask[1:-1, 1:-1, 1:-1] * 0 Nz, Ny, Nx = sumAll.shape for i in [-1, 0, 1]: for j in [-1, 0, 1]: for k in [-1, 0, 1]: if np.sqrt(abs(i) + abs(j) + abs(k)) > 1.5: continue sumAll = sumAll + extendMask[1 + i:Nz + 1 + i, j + 1:Ny + 1 + j, k + 1: Nx + 1 + k] return sumAll def sumEdgeByCon3(self, extendMask): # 27 in total raw = extendMask[1:-1, 1:-1, 1:-1] Nz, Ny, Nx = raw.shape sumAll = raw * 0 for i in [-1, 0, 1]: for j in [-1, 0, 1]: for k in [-1, 0, 1]: sumAll = sumAll + extendMask[1 + i:Nz + 1 + i, j + 1:Ny + 1 + j, k + 1: Nx + 1 + k] return sumAll ############ def setOutputPath(self, targetPath): self.processPath = targetPath doFITS.checkPath(self.processPath ) self.getOutputNames()### redirect the outpyt path def writeConfig(self,configFile=None): """ #### :return: """ #### first read then add if configFile is not None and os.path.isfile( configFile ): self.configFile=configFile self.config = configparser.ConfigParser() self.config.read( self.configFile ) self.config[self.regionName] = {} self.config[self.regionName] [self.rawCOFITSConf ] = self.rawCOFITS self.config[self.regionName] [self.rmsFITSConf ] = self.rmsFITS self.config[self.regionName] [self.dbscanPathConf ] = self.processPath self.config[self.regionName][self.extraTagConf ] = self.extraTag #### bug fixed#qzyan 2023/10/31 print("Writting config of ",self.regionName, " to " , self.configFile) with open(self.configFile, 'w') as cf: self.config.write(cf) def set_cloud(self,cloudID): self.getCloudIntName(cloudID) def set_region(self,regionName): """ :param regionName: :return: """ if regionName is None: return self.regionName = regionName self.config = configparser.ConfigParser() self.config.read( self.configFile ) if self.regionName not in self.config.sections(): print("No configuration of {} found...please add the configure...{}".format(self.regionName,self.configFile )) else: self.rmsFITS = self.config[self.regionName][self.rmsFITSConf] self.rawCOFITS = self.config[self.regionName][self.rawCOFITSConf] self.extraTag = self.config[self.regionName][self.extraTagConf ] if not os.path.isfile(self.rawCOFITS): print("The CO raw file does not exist for region {}, stop reading configure files.".format( self.regionName) ) return if self.config[self.regionName][self.dbscanPathConf] != "" : self.setOutputPath( self.config[self.regionName][self.dbscanPathConf] ) else: defaultPath = os.path.join(self.defaultDBSCANpath,self.regionName) self.setOutputPath( defaultPath ) #using default path self.initDBSCAN() self.getCloudIntPath() #using defaults def getOutputNames(self): """ create all the outpute names of files, even they have not been created :return: """ print("Initializing {}....".format(self.regionName)) if self.rawCOFITS is None: print("CO FITS is needed") return if not os.path.isfile(self.rawCOFITS): print("The raw COFITS does not exist!") return if self.averageRMS is None and self.rmsFITS is None: print("Please set the average RMS or produce the rmsFITS!") return if self.rmsFITS is not None and not os.path.isfile(self.rmsFITS): print("The rms FITS does not exist!") return self.rawCOFITS=os.path.abspath(self.rawCOFITS) if self.rmsFITS is not None: self.rmsFITS=os.path.abspath(self.rmsFITS) self.processPath=os.path.abspath(self.processPath) doFITS.checkPath(self.processPath) #### baseName = os.path.basename(self.rawCOFITS) saveTag = self.getDBSCANTag() preFix = os.path.join(self.processPath, baseName[0:-5] + saveTag) if self.regionName is None: self.regionName = baseName[0:-5] + saveTag # used to save configures self.rawLabelFITSName = preFix + ".fits" self.cleanLabelFITSName = preFix + "_Clean.fits" self.rawCatName = preFix + ".fit" self.cleanCatName = preFix + "_Clean.fit" self.noiseFITSname = preFix + "_noise.fits" self.getInfo() self.getCloudIntPath() def getInfo(self): """ basis information from header :return: """ dataCO,headCO= doFITS.readFITS(self.rawCOFITS) #### at this stage rawCOFITS is self.pixelSize=headCO["CDELT2"]*60#in armin self.velResolution = headCO["CDELT3"]/1000. #in km/s del dataCO gc.collect() ##### # if self.rmsFITS is None: #no rms fits is provided def getDBSCANTag(self): """ return :param cutoff: #the sigma to which the overall datacube will be cut :param minpts: :param conType: #thre are :return: """ saveName = "dbscanS{}P{}Con{}".format(self.cutoff_sigma, self.minPts, self.connectivity)+self.extraTag return saveName def getUniformRMSData(self, averageRMS, Ny, Nx): return averageRMS + np.zeros((Ny, Nx)) def getCoreArray(self, extendMask): """ :param maskArray: :param contype: :return: """ if self.connectivity == 1: coreArray = self.sumEdgeByCon1(extendMask) if self.connectivity == 2: coreArray = self.sumEdgeByCon2(extendMask) if self.connectivity == 3: coreArray = self.sumEdgeByCon3(extendMask) return coreArray def computeDBSCAN(self, splitN=1): """ There are two steps of computing DBSCAN, compute the DBSCAN label, 2 extract catlog, select post catalog and produce clean fits, would be doe seperatrely :param COdata: :param COHead: :param min_sigma: cutoff :param min_pix: :param connectivity: :param region: :param getMask: :param savePath: :param mimicDendro: :param rmsFITS: :param inputRMS: :return: """ # pass if self.processPath is None: if not os.path.isdir(self.backupPath): os.mkdir(self.backupPath ) self.processPath=self.backupPath self.getOutputNames() if self.rawCOFITS is None: print("rawCOFITS need to be provided, stoping...") return COdata, COHead = doFITS.readFITS(self.rawCOFITS) if self.controlMaskFITS is not None: print("Masking data cube with extra mask file ",self.controlMaskFITS) forceMask,headForeMask = doFITS.readFITS( self.controlMaskFITS ) COdata[forceMask==0]=0 del forceMask gc.collect() if len(COdata.shape) == 4: COdata = COdata[0] Nz, Ny, Nx = COdata.shape if self.averageRMS is None and self.rmsFITS is None: print( "You need either use averaged rms, or use rmsFITS, which can be accurate to each spectrum. Please check you input") if self.rmsFITS == None: # use average rms, te rmsFITS has a higher priority rmsData = self.getUniformRMSData(self.averageRMS, Ny, Nx) else: print("Using rms fits...", self.rmsFITS) rmsData, rmsHead = myFITS.readFITS(self.rmsFITS) self.rmsData = rmsData # record the rms Data rmsCOData = COdata / rmsData del COdata gc.collect() goodValues = rmsCOData >= self.cutoff_sigma if self.maskFirst: extraMask = self.getIsolatedMask(splitN=splitN) goodValues = extraMask # np.logical_and(goodValues,) del rmsCOData gc.collect() s = generate_binary_structure(3, self.connectivity) overlapping = 10 indexEdges = np.linspace(0, Nx - 1, splitN + 1, dtype=np.int32) coreArray = np.zeros_like(goodValues, dtype=np.int32) for i in range(len(indexEdges) - 1): #### print("doing", i) indexRange = [indexEdges[i] - overlapping, indexEdges[i + 1] + overlapping + 1] leftCut = overlapping rightCut = indexEdges[i + 1] - indexEdges[i] + 1 + overlapping if i == 0: indexRange[0] = 0 leftCut = 0 rightCut = indexEdges[i + 1] - indexEdges[i] + 1 if i == len(indexEdges) - 1: indexRange[1] = indexEdges[i + 1] + 1 rightCut = indexEdges[i + 1] - indexEdges[i] + 1 + overlapping subGood = goodValues[:, :, indexRange[0]: indexRange[1]] subNz, subNy, subNx = subGood.shape subExtend = np.zeros([subNz + 2, subNy + 2, subNx + 2], dtype=np.int32) subExtend[1:-1, 1:-1, 1:-1] = subGood subCoreArray = self.getCoreArray(subExtend) coreArray[:, :, indexEdges[i]: indexEdges[i + 1] + 1] = subCoreArray[:, :, leftCut: rightCut] del subCoreArray del subExtend gc.collect() # coreArray = self.getCoreArray(extendMask) coreArray = coreArray >= self.minPts gc.collect() coreArray[ ~goodValues] = False # nan could be, #remove falsely, there is a possibility that, a bad value may have lots of pixels around and clould be # coreArray = coreArray + 0 labeled_core, num_features = label(coreArray, structure=s,output=np.int32 ) # first label core, then expand, otherwise, the expanding would wrongly connected print(num_features, "features found!") del coreArray gc.collect() if self.saveCore: """ """ saveCoreName = self.getLabelCoreFITSName() fits.writeto(saveCoreName, labeled_core, header=COHead, overwrite=True) self.coreFITS = saveCoreName selectExpand = np.logical_and(labeled_core == 0, goodValues) del goodValues gc.collect() # expand labeled_core # coreLabelCopy=labeled_core.copy() expandTry = dilation(labeled_core, s) # first try to expand, then only keep those region that are not occupied previously # it is possible that a molecular cloud may have less pixN than 8, because of the closeness of two # only expanded, one time labeled_core[selectExpand] = expandTry[selectExpand] fits.writeto(self.rawLabelFITSName, labeled_core, header=COHead, overwrite=True) del labeled_core gc.collect() def getExtend(self,dataCube): """ :param dataCube: :return: """ subNz, subNy, subNx = dataCube.shape subExtend = np.zeros([subNz + 2, subNy + 2, subNx + 2], dtype=np.int32) subExtend[1:-1, 1:-1, 1:-1] = dataCube return subExtend def getLabelCoreFITSName(self): """ used to get the DBSCAN label fits name, when you do not want to rerun the DBSCAN process :return: """ baseName = os.path.basename(self.rawCOFITS) saveTag = self.getDBSCANTag() saveLabelFITSNameCore = os.path.join(self.processPath, baseName[0:-5] + saveTag + "_core.fits") return saveLabelFITSNameCore def getIndices(self, Z0, Y0, X0, values1D, choseID,pixN=None ): if pixN is not None: #np.searchsorted(values1D,choseID) #np.where(values1D == choseID) endIndex = pixN cX0 = X0[0:endIndex ] # X0[cloudIndices] cY0 = Y0[0:endIndex ] # Y0[cloudIndices] cZ0 = Z0[0:endIndex ] # Z0[cloudIndices] return tuple([cZ0, cY0, cX0]) , values1D[endIndex+1:], Z0[endIndex:],Y0[endIndex:],X0[endIndex:] else: cloudIndices = np.where(values1D == choseID) cX0 = X0[cloudIndices] cY0 = Y0[cloudIndices] cZ0 = Z0[cloudIndices] return tuple([cZ0, cY0, cX0]) def rmsmap(self, outPUT=None, overwrite=True): """ 3d rms :param outPUT: :param overwrite: :return: """ if self.rawCOFITS is None: print("rawCOFITS need to be provided, stopping...") return COdata, COHead = doFITS.readFITS(self.rawCOFITS) if len(COdata.shape) == 4: COdata = COdata[0] if outPUT is None: writeName = "rmsmap.fits" else: writeName = outPUT fileExist = os.path.isfile(writeName) if overwrite and fileExist: os.remove(writeName) Nz, Ny, Nx = COdata.shape rmsData = np.zeros_like(COdata, dtype=np.float32) for i in range(Nz): channelI = COdata[i] negativeValues = channelI[channelI < 0] sigma = np.std(negativeValues) / np.sqrt(1 - 2. / np.pi) print(sigma) rmsData[i, :, :] = sigma fits.writeto(writeName, rmsData, header=COHead) return fits.open(writeName)[0] def getEmptyCat(self, getEmptyRow=False): """ :return: """ newTB = Table(names=("_idx" ,), dtype=('i8', ) ) ### newTB.add_column( Column( name="area_exact", dtype='f8',unit= u.arcmin*u.arcmin ) ) ##### newTB.add_column( Column( name="v_cen", dtype='f8',unit= u.km/u.s ) ) #### newTB.add_column( Column( name="v_rms", dtype='f8',unit= u.km/u.s ) ) #### newTB.add_column(Column( name="x_cen", dtype='f8',unit= u.deg ) ) #### newTB.add_column( Column( name="y_cen", dtype='f8',unit= u.deg ) ) #### newTB.add_column( Column( name="l_rms", dtype='f8',unit= u.deg ) ) #### newTB.add_column( Column( name="b_rms", dtype='f8',unit= u.deg ) ) newTB.add_column( Column( name="pixN", dtype='i8' ) ) newTB.add_column( Column( name="area_accurate", dtype='f8',unit= u.arcmin*u.arcmin ) ) #### newTB.add_column( Column( name="peak", dtype='f8',unit= u.Kelvin ) ) #### newTB.add_column( Column( name="peakB", dtype='i8' ) ) #### newTB.add_column( Column(name="peakL", dtype='i8' ) ) #### newTB.add_column( Column( name="peakV", dtype='i8' ) ) ####LBV of peak newTB.add_column( Column( name="peakBworld", dtype='f8' ) ) #### newTB.add_column( Column(name="peakLworld", dtype='f8' ) ) #### newTB.add_column( Column( name="peakVworld", dtype='f8' ) ) #### newTB.add_column( Column( name="sum", dtype='f8',unit= u.Kelvin ) ) #### newTB.add_column( Column( name="flux", dtype='f8',unit= u.Kelvin* u.km/u.s *u.arcmin*u.arcmin ) ) #### newTB.add_column( Column( name="lineWidth", dtype='f8',unit= u.km/u.s ) ) #### newTB.add_column(Column( name="meanSNR", dtype='f8')) newTB.add_column(Column( name="meanT", dtype='f8')) newTB.add_column( Column( name="allChannel", dtype='i8' ) ) ######### startLindexCol="startLindex" startBindexCol="startBindex" startVindexCol="startVindex" newTB.add_column(Column( name= startLindexCol, dtype='f8')) newTB.add_column(Column( name= startBindexCol, dtype='f8')) newTB.add_column(Column( name= startVindexCol, dtype='f8')) endLindexCol="endLindex" endBindexCol="endBindex" endVindexCol="endVindex" newTB.add_column(Column( name= endLindexCol, dtype='f8')) newTB.add_column(Column( name= endBindexCol, dtype='f8')) newTB.add_column(Column( name= endVindexCol, dtype='f8')) startLCol="startL" startBCol="startB" startVCol="startV" newTB.add_column(Column( name= startLCol, dtype='f8')) newTB.add_column(Column( name= startBCol, dtype='f8')) newTB.add_column(Column( name= startVCol, dtype='f8')) endLCol="endL" endBCol="endB" endVCol="endV" newTB.add_column(Column( name= endLCol, dtype='f8')) newTB.add_column(Column( name= endBCol, dtype='f8')) newTB.add_column(Column( name= endVCol, dtype='f8')) #### maximum and minimum Z of cloud, exact Index newTB.add_column( Column( name="minZ", dtype='i8' ) ) newTB.add_column( Column( name="maxZ", dtype='i8' ) ) #### newTB.add_column( Column( name="has22", dtype='i4' ) ) #### newTB.add_column( Column( name="touchLB", dtype='i4' ) ) newTB.add_column( Column( name="touchL", dtype='i4' ) ) newTB.add_column( Column( name="touchB", dtype='i4' ) ) newTB.add_column( Column( name="touchV", dtype='i4' ) ) newTB.add_column( Column( name="complete", dtype='i4' ) ) if getEmptyRow: newTB.add_row() return newTB[0] return newTB def getDBSCANCatlog(self, doClean=True): """ Extract catalog from label fits, the minPix and rms is only used for saving #need a TB Model, if TB model does not exist, better downlooad from online This catalog has already been cleaned four criteria minPix" if you want to change the catlog selection criteria, please see setCatalogSelectionCariteria :param labelArray: :param head: :return: """ #### Add a column of line width, which will be used in statistics #### self.getOutputNames() if self.rawLabelFITSName is None or not os.path.isfile(self.rawLabelFITSName): print("The DBSCAN label FITS is needed ") return dataCluster, headCluster = myFITS.readFITS(self.rawLabelFITSName) # calculate the area of one pixel sizeL = headCluster["CDELT1"] * 60 sizeB = headCluster["CDELT2"] * 60 sizeV = headCluster["CDELT3"] /1000. #km/s #### sizeB = abs(sizeB) # in arcmin sizeL = abs(sizeL) # in arcmin pixelArea = sizeB * sizeL Nz, Ny, Nx = dataCluster.shape if self.rawCOFITS is None: print("A COFITS is needed ") return dataCO, headCO = myFITS.readFITS(self.rawCOFITS) # if len(dataCO.shape) == 4: dataCO = dataCO[0] if self.rmsFITS is not None: print("Using rms fits...", self.rmsFITS) rmsData, rmsHead = myFITS.readFITS(self.rmsFITS) self.rmsData = rmsData if len(rmsData.shape) ==3 and rmsData.shape[0]==1 : #### the first cooridiane is wrong self.rmsData=rmsData[0] else: if self.averageRMS is not None and self.rawCOFITS is not None: rmsData = self.getUniformRMSData(self.averageRMS, Ny, Nx) self.rmsData = rmsData if self.rmsData is None: print("A rmsData is needed to select peakSgima") return minV = np.nanmin(dataCluster[0]) # noise mask value wcsCloud = WCS(headCluster, naxis=3) #### wcsCloud.wcs.bounds_check(False, False) # for molecular clouds toward the anti Galacic center clusterIndex1D = np.where(dataCluster > minV) clusterValue1D = dataCluster[clusterIndex1D] if 1: #### sort index sortIndices = np.argsort(clusterValue1D) iZ, iY, iX = clusterIndex1D iZ_S = iZ[sortIndices] iY_S = iY[sortIndices] iX_S = iX[sortIndices] clusterValue1D = clusterValue1D[sortIndices] clusterIndex1D = (iZ_S, iY_S, iX_S) del dataCluster gc.collect() dataZero = np.zeros_like(dataCO) # #### #### zeroProjection = np.zeros((Ny, Nx)) # one zero channel, used to get the projection area and zeroProjectionExtend = np.zeros((Ny + 1, Nx + 1)) idCol = "_idx" # count all clusters # ids,count=np.unique(dataCluster,return_counts=True ) ids, count = np.unique(clusterValue1D, return_counts=True) #print( "{} DBSCAN clusters found".format( len(ids) ) ) # print "Total number of Good Trunks? ",len(GoodIDs) GoodIDs = ids GoodCount = count print( "{} DBSCAN clusters found".format( len(ids) ) ) # print "Total number of Good Trunks? ",len(GoodIDs) # dataCO,headCO=doFITS.readFITS( CO12FITS ) # zeroP # remove any cloud with voxels less than minPix selectVox = GoodCount >= self.minVox #criteria 1 checkN = len(GoodIDs) if doClean: GoodCount = GoodCount[selectVox] GoodIDs = GoodIDs[selectVox] checkN = len(GoodIDs) print(checkN, "cluster has more than {} voxels.".format(self.minVox)) allCounts = np.repeat(count,count) if len(allCounts) !=len(clusterValue1D): print( len(allCounts),len(clusterValue1D),"should be equal???Terminating ..... " ) aaaaaaaaa startLindexCol="startLindex" startBindexCol="startBindex" startVindexCol="startVindex" endLindexCol="endLindex" endBindexCol="endBindex" endVindexCol="endVindex" startLCol="startL" startBCol="startB" startVCol="startV" endLCol="endL" endBCol="endB" endVCol="endV" if 1: #### remove useless indecies selectGoodIndices = allCounts>= self.minVox iZ, iY, iX = clusterIndex1D iZ_S = iZ[selectGoodIndices] iY_S = iY[selectGoodIndices] iX_S = iX[selectGoodIndices] clusterValue1D = clusterValue1D[selectGoodIndices] clusterIndex1D = (iZ_S, iY_S, iX_S) ### remove those clouds with peak less that self.minPeakSigma ### getCOvalues coValues1D = dataCO[clusterIndex1D] rmsValues1D = self.rmsData[(iY_S,iX_S)] SNRValues1D = coValues1D/rmsValues1D emptyTB = self.getEmptyCat( ) emptyTB.add_row() if 1: #sort indices clusterValue1DSearch = clusterValue1D.copy() clusterIndex1DSearch = (clusterIndex1D[0].copy(),clusterIndex1D[1].copy(),clusterIndex1D[2].copy() ) Z0Search, Y0Search, X0Search = clusterIndex1DSearch SNRValues1DSearch = SNRValues1D.copy() while 1 : copyTB = emptyTB.copy() emptyTB= vstack((copyTB,emptyTB)) if len(emptyTB)> checkN: break widgets = ['Calculating cloud parameters: ', Percentage(), ' ', Bar(marker='0', left='[', right=']'), ' ', ETA(), ' ', FileTransferSpeed()] # see docs for other options newTB = self.getEmptyCat(getEmptyRow=False) catTB = newTB pbar = ProgressBar(widgets=widgets, maxval=len(GoodIDs)) pbar.start() idRecord=[] indicesRecord=[] for i in range(len(GoodIDs)): pbar.update(i) # i would be the newID newID = GoodIDs[i] pixN = GoodCount[i] newRow =emptyTB[i] newRow[idCol] = newID #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! becare full cloudIndex ,clusterValue1DSearch , Z0Search,Y0Search,X0Search= self.getIndices(Z0Search, Y0Search, X0Search, clusterValue1DSearch, newID, pixN= pixN ) cloudMaxSNR = np.max(SNRValues1DSearch[0:pixN]) if cloudMaxSNR4: # if 4 in sum22: newRow["has22"] = 1 else: # Serious bug newRow["has22"] = 0 if self.hasBeam and newRow["has22"] == 0 and doClean: # #criteria 4 zeroProjection[projectIndex] = 0 zeroProjectionExtend[0:-1, 0:-1] = zeroProjection continue #### tohere? # newRow["has22"] = 0 # calculate the accurate, this is only used for data at high Galactic latitude, and useless form the MWISP survey # to increase the speed, this part is commentted # indexB2D,indexL2D=np.where(zeroProjection==1 ) # _,BS2D, LS2D = wcsCloud.wcs_pix2world(indexL2D,indexB2D, 0, 0) # area_accurate=np.sum( np.cos( np.deg2rad(BS2D) ) )*0.25 # newRow["area_accurate"]= area_accurate sumCO = np.sum(coValues, dtype=np.float64) # float32 is not enough for large molecular clouds Vcen, Vrms = doFITS.weighted_avg_and_std(cloudV, coValues) Bcen, Brms = doFITS.weighted_avg_and_std(cloudB, coValues) Lcen, Lrms = doFITS.weighted_avg_and_std(cloudL, coValues) # calculate the exact area # LBcore = zip(cloudB, cloudL) # pixelsN= {}.fromkeys(LBcore).keys() #len( set(LBcore) ) # area_exact=len(pixelsN)*0.25 #arc mins square area_exact = np.sum(zeroProjection) * pixelArea newRow["area_exact"] = area_exact # dataClusterNew[cloudIndex] =newID # save values newRow["pixN"] = pixN newRow["peak"] = peak newRow["peakV"] = peakV newRow["peakB"] = peakB newRow["peakL"] = peakL newRow["meanSNR"] = meanSNR # newRow["peak2"]= peak2 newRow["flux"] = sumCO * pixelArea * sizeV #### newRow["sum"] = sumCO newRow["x_cen"], newRow["y_cen"], newRow["v_cen"] = wcsCloud.wcs_pix2world(Lcen, Bcen, Vcen, 0) newRow["v_cen"] = newRow["v_cen"] / 1000. newRow["peakLworld"], newRow["peakBworld"], newRow["peakVworld"] = wcsCloud.wcs_pix2world(peakL, peakB, peakV, 0) newRow["peakVworld"] = newRow["peakVworld"] / 1000. dv = headCluster["CDELT3"] / 1000. # km/s dl = abs(headCluster["CDELT1"]) # deg newRow["v_rms"] = Vrms * dv newRow["l_rms"] = Lrms * dl newRow["b_rms"] = Brms * dl # _, Nchan=np.unique( cloudV, return_counts=True) # newRow["Nchannel"] = np.max(P3Value)# if there is a three consecutive spectra in the cloud newRow["allChannel"] = len(diffVs) ##################################################below are equivalenet legnth dataZero[cloudIndex] = dataCO[cloudIndex] # cropThe cloudRange minY = np.min(cloudB) maxY = np.max(cloudB) ########### minX = np.min(cloudL) maxX = np.max(cloudL) ########### minZ = np.min(cloudV) maxZ = np.max(cloudV) newRow["minZ"] = minZ newRow["maxZ"] = maxZ #### fanny save ######### #### #### #### rawCrop = dataCO[:, minY:maxY + 1, minX:maxX + 1] zeroprojMask = zeroProjection[minY:maxY + 1, minX:maxX + 1] rawCrop = rawCrop*zeroprojMask averageSpectraRaw = np.nansum(rawCrop, axis=(1, 2))/projectionPixN #### find the width of linewidth based on real spectra? #### should we save this spectra? then investigate them if self.saveSpec: saveName= self.getSpecSaveName(newID) np.save(saveName,averageSpectraRaw) cleanCrop = dataZero[:, minY:maxY + 1, minX:maxX + 1] averageSpectraClean = np.nansum(cleanCrop, axis=(1, 2)) # count the number spectra totalSpectral = projectionPixN / 1. meanSpectral = averageSpectraClean / 1. / totalSpectral # if 0: # savefileName = saveSpectralPath + "{}_{}Spectral".format(regionName, testID) # np.save(savefileName, [vAxis, meanSpectral]) spectraPeak = np.max(meanSpectral) area = dv * np.sum(meanSpectral) eqLineWidth = area / spectraPeak newRow["lineWidth"] = eqLineWidth #### record to the sake of clean ##### record indices here idRecord.append(newID) indicesRecord.append(cloudIndex) #catTB.add_row(newRow) #### clear middle fits dataZero[cloudIndex] = 0 zeroProjection[projectIndex] = 0 zeroProjectionExtend[0:-1, 0:-1] = zeroProjection pbar.finish() # save the clouds emptyTB = emptyTB[emptyTB["area_exact"]>0] # fits.writeto(self.regionName+"NewCloud.fits", dataClusterNew,header=headCluster,overwrite=True ) if doClean: print("{} clouds found by DBSCAN!!!!".format(len(emptyTB))) emptyTB.write(self.cleanCatName, overwrite=True) #### what do we know self.saveCloudID = idRecord self.saveCloudIndex = indicesRecord # else: emptyTB.write(self.rawCatName, overwrite=True) def getSpecSaveName(self,cloudID): saveName = "SpecRaw{}.npy".format(cloudID) saveName = os.path.join(self.cloudSpecPath, saveName) #### return saveName def cleanTB(self): """ select clouds that satifies the selection criteria :return: """ if self.rawCatName is None: print("no catalog file provided") return if self.rmsData is None: print("You need to provide the rms FITS data to select peak accordingly") rawTB = Table.read(self.rawCatName) # first voxel filterTB = rawTB[rawTB["pixN"] >= self.minVox] filterTB = filterTB[filterTB["has22"] == self.hasBeam] # select by minCkannel filterTB = filterTB[filterTB["allChannel"] >= self.minChannel] # second peak peakB = list(map(int, filterTB["peakB"])) peakL = list(map(int, filterTB["peakL"])) peakCoordinate = tuple([peakB, peakL]) rmsArray = self.rmsData[peakCoordinate] filterTB = filterTB[filterTB["peak"] >= self.minPeakSigma * rmsArray] filterTB.write(self.cleanCatName, overwrite=True) def getCloudIntName(self,cloudID): """ :return: """ self.runID = cloudID self.saveIntFITSname = os.path.join(self.cloudIntPath, "Cloud{}_int.fits".format(self.runID)) self.savePureIntFITSname = os.path.join(self.cloudIntPath , "Cloud{}_PureInt.fits".format(self.runID)) self.saveMaskFITSname = os.path.join(self.cloudIntPath , "Cloud{}_mask.fits".format( self.runID )) self.saveForegroundFITSname = os.path.join(self.cloudIntPath , "Cloud{}_fore.fits".format( self.runID )) def getCloudIntPath(self ): #print(self.processPath,"??????????????????????????") #doFITS.checkPath(self.processPath) if self.cloudIntPath is None: self.cloudIntPath = os.path.join(self.processPath, self.intFolderName ) doFITS.checkPath(self.cloudIntPath) ########## self.cloudSpecPath = os.path.join(self.processPath, self.specFolderName ) doFITS.checkPath(self.cloudSpecPath) #### else: pass # the output path has been assigned directely def produceCloudIntFITS(self, inputTB=None, outputPath=None, pureInt=False, foreground=False, complete=True, areaControl=True, localVindx=None, onID=None, doCrop=True ,cropExtend=None): """ :param inputTB: in case, a sub set of the DBSCAN clean TB is intested, other would be ignore :param outputPath: outPath of int and mask FITS :param pureInt: produce a int map, with noise voxels masked :param foreground: local V is need to produce a foreground emission map based on the velocity, used for distance examination :param localV: #in km/s, used to produce foreground fits :param onID: # in case only one cloud is interested :param areaThresh: The minimum area for distance examination :param cropExtend: crop The int FITS, incase the wholse area is large :return: """ if self.cleanCatName is None or not os.path.isfile( self.cleanCatName ) : print("Please get the DBSCAN data ready 1") return if self.rawCOFITS is None or not os.path.isfile( self.rawCOFITS ) : print("Please get the DBSCAN data ready 2") return if self.cleanLabelFITSName is None or not os.path.isfile( self.cleanLabelFITSName ) : print("Please get the DBSCAN data ready 3") return if self.processPath is None: print("Please get the DBSCAN data ready 4") return if cropExtend is not None: self.cropExtend=cropExtend if outputPath is None : self.getCloudIntPath() else: doFITS.checkPath(outputPath) self.cloudIntPath = outputPath #self.getCloudIntPath() print("Saving cloud int FITS to ",self.cloudIntPath) if localVindx is not None and localVindx!="auto": self.localVindex = localVindx ##### TB=Table.read(self.cleanCatName) if inputTB is not None : TB = inputTB if complete: #### usually this is true TB = TB[TB["complete"] > 0] # only complete ones TB.sort("area_exact") TB.reverse() if areaControl:# too many or not, it depends ### remove small molecular clouds, other wise there are too many TB = TB[TB["area_exact"] >= self.areaThresh * 3600] ## Remove small molecular clouds print("After removing smaller molecular clouds, {} left".format(len(TB) ) ) else: print("All complete molecular clouds, {} left".format(len(TB) ) ) if onID is not None: #if only cares about one molecular cloud TB=TB[ TB["_idx"] == onID ] dataCluster, headCluster = myFITS.readFITS(self.cleanLabelFITSName) Nz,Ny,Nx= dataCluster.shape WCSco = WCS(headCluster) #### if localVindx == "auto": # calculate local vindex automatically _,_,localV=WCSco.wcs_world2pix(180,0,0,0) self.localVindex= int( np.round(localV) ) self.localVindex=max([0,self.localVindex ]) self.localVindex=min([Nz-1,self.localVindex ]) #### set localVindex, at the index of 0 km/s clusterIndex1D = np.where(dataCluster > 0) clusterValue1D = dataCluster[clusterIndex1D] Z0, Y0, X0 = clusterIndex1D del dataCluster # release the memory gc.collect() dataCO, headCO = myFITS.readFITS( self.rawCOFITS ) coSpec0, vaxis0 = doFITS.getSpectraByIndex(dataCO, headCO, 0, 0) #### used to produce precise foreground emission dataCOMask = np.zeros_like(dataCO, dtype=np.float32) dataCOMask[ clusterIndex1D ] = dataCO[clusterIndex1D] #################### velsolution = vaxis0[1] - vaxis0[0] ################ projection0 = np.zeros_like(dataCO[0] ) widgets = ['Integrating cloud fits:', Percentage(), ' ', Bar(marker='0', left='[', right=']'), ' ', ETA(), ' ', FileTransferSpeed()] # see docs for other options pbar = ProgressBar(widgets=widgets, maxval=len(TB)) pbar.start() indexRun = 0 zeroCube = np.zeros_like(dataCO, dtype=np.float32) for eachDBRow in TB: indexRun = indexRun + 1 pbar.update(indexRun) cloudID = eachDBRow["_idx"] self.getCloudIntName(cloudID) #### cloudIndex = self.getIndices(Z0, Y0, X0, clusterValue1D, cloudID) iz, iy, ix = cloudIndex projectIndices = tuple([iy, ix]) ############ determine v range centerV = eachDBRow["v_cen"] dv = eachDBRow["v_rms"] vRange = [centerV - dv * 3, centerV + dv * 3] startV0, endVindex = self.getVindexRange(vaxis0, vRange) cropCOCube = dataCO[startV0:endVindex + 1] ###int fits sumCO = np.sum(cropCOCube, axis=0, dtype=float) * velsolution fits.writeto(self.saveIntFITSname , sumCO, header=headCluster, overwrite=True) ##mask fits projection0[projectIndices] = 1 fits.writeto(self.saveMaskFITSname, projection0, header=headCluster, overwrite=True) ##foreground fits if pureInt: # usually do no use this #startV0 = min(iz) #endV0 = max(iz) zeroCube[cloudIndex] = dataCO[cloudIndex] cropCOCubePure = zeroCube[startV0:endVindex + 1] sumCOPure = np.sum(cropCOCubePure, axis=0, dtype=float) * velsolution #print(self.savePureIntFITSname,"What is this?") fits.writeto(self.savePureIntFITSname, sumCOPure, header=headCluster, overwrite=True) #print("Why this does not work") if foreground: # the way of generating foreground fits is different, the following code is used to produce foreground fits for Q2 #### using the close one if abs(endVindex-self.localVindex) < abs(startV0-self.localVindex): starForeI = np.min( [endVindex, self.localVindex ] ) endForeI = np.max( [endVindex, self.localVindex ] ) #print(starForeI,endForeI,"Fore range",startV0,endVindex,"cloud range") foreCOCube = dataCOMask[starForeI : endForeI ] else: starForeI = np.min( [startV0, self.localVindex ] ) endForeI = np.max( [startV0, self.localVindex ] ) #print(starForeI,endForeI ,"Fore range",startV0,endVindex,"cloud range") foreCOCube = dataCOMask[starForeI : endForeI ] sumForeground = np.sum(foreCOCube, axis=0, dtype=float) * velsolution fits.writeto(self.saveForegroundFITSname , sumForeground, header=headCluster, overwrite=True) if doCrop: #crop the fits cropLrange, cropBrange = doFITS.getLBrangeWithMask(self.saveMaskFITSname,extendByDeg= self.cropExtend ) doFITS.cropFITS2D(self.saveIntFITSname,outFITS=self.saveIntFITSname,Lrange=cropLrange,Brange=cropBrange, overWrite=True) doFITS.cropFITS2D(self.saveMaskFITSname,outFITS=self.saveMaskFITSname,Lrange=cropLrange,Brange=cropBrange, overWrite=True) if os.path.isfile(self.savePureIntFITSname): doFITS.cropFITS2D(self.savePureIntFITSname,outFITS=self.savePureIntFITSname,Lrange=cropLrange,Brange=cropBrange, overWrite=True) if os.path.isfile(self.saveForegroundFITSname): doFITS.cropFITS2D(self.saveForegroundFITSname,outFITS=self.saveForegroundFITSname,Lrange=cropLrange,Brange=cropBrange, overWrite=True) projection0[projectIndices] = 0 zeroCube[cloudIndex] = 0 # zeroCluster[cloudIndex] = cloudID # remove this cluster and do not record this cluster pbar.finish() def getVindexRange(self, vaxis, vRange): """ :param vaxis: :param vRange: :return: """ minValue = min(vRange) maxValue = max(vRange) indexV0 = doFITS.find_nearestIndex(vaxis, minValue) indexV1 = doFITS.find_nearestIndex(vaxis, maxValue) return [indexV0, indexV1] def produceCleanFITS(self): """ remove noise cluster according to the label fits and catFITS name :return: """ if self.rawLabelFITSName is None: print("You need to provde the label FITS ") return if self.cleanCatName is not None and os.path.isfile(self.cleanCatName): targetTBFile = self.cleanCatName # table has alread been cleaned if not os.path.isfile(self.cleanCatName) and not os.path.isfile(self.rawCatName): print("You need to provde a raw or a cleaned fits first to produce a clean fits") return if os.path.isfile(self.rawCatName) and not os.path.isfile(self.cleanCatName): # catalog has been calculated but not cleaned, need to clean table first, then do clean self.cleanTB() targetTBFile = self.cleanCatName # table has alread been cleaned TB = Table.read(targetTBFile) dataCluster, headCluster = doFITS.readFITS(self.rawLabelFITSName) noiseLabel = np.min(dataCluster[0]) clusterIndex1D = np.where(dataCluster > noiseLabel) clusterValue1D = dataCluster[clusterIndex1D] Z0, Y0, X0 = clusterIndex1D # if self.saveCloudID is not None and self.saveCloudIndex is not None: print("Using saved index to speed up!") widgets = ['Cleaning label fits:', Percentage(), ' ', Bar(marker='0', left='[', right=']'), ' ', ETA(), ' ', FileTransferSpeed()] # see docs for other options pbar = ProgressBar(widgets=widgets, maxval=len(TB)) pbar.start() indexRun = 0 returnCluster = np.zeros_like(dataCluster) #+ noiseLabel if self.saveCloudID is not None and self.saveCloudIndex is not None: totalN = len(TB) for i in range(totalN): cloudIndex=self.saveCloudIndex[i] cloudID = self.saveCloudID[i ] returnCluster[cloudIndex] = cloudID # remove this cluster and do not record this cluster else: for eachDBRow in TB: indexRun = indexRun + 1 pbar.update(indexRun) cloudID = eachDBRow["_idx"] cloudIndex = self.getIndices(Z0, Y0, X0, clusterValue1D, cloudID) returnCluster[cloudIndex] = cloudID # remove this cluster and do not record this cluster pbar.finish() fits.writeto(self.cleanLabelFITSName, returnCluster, header=headCluster, overwrite=True) dataCO,headCO=doFITS.readFITS(self.rawCOFITS) if len( dataCO.shape )==4: dataCO = dataCO[0] dataCO[ np.where(returnCluster>0) ]=np.nan if self.writeNoise: fits.writeto(self.noiseFITSname, dataCO, header=headCluster, overwrite=True) def produceMask(self, COFITS, labelFITS, outFITS=None): """ produce a masked COFITS, with labelsFITS :param COFITS: :param labelFITS: :return: """ dataCO, headCO = doFITS.readFITS(COFITS) dataLabel, headLabel = doFITS.readFITS(labelFITS) noiseLabel = np.min(dataLabel[0]) if outFITS is None: outFITS = "mask_" + COFITS # maskCO dataCO[dataLabel == noiseLabel] = np.nan fits.writeto(outFITS, dataCO, header=headCO, overwrite=True) def getLBVrangeOfClouds(self): """ A function of get the LBV range of clouds """ cloudTB=Table.read(self.cleanCatName) catNameWithLBVrange = doFITS.addSuffix(self.cleanCatName,"_WithLBVRange") dataCluster, headCluster = myFITS.readFITS(self.cleanLabelFITSName) minV = np.nanmin(dataCluster[0]) wcsCloud = WCS(headCluster, naxis=3) ### usd to convert index to lbv clusterIndex1D = np.where(dataCluster > minV) clusterValue1D = dataCluster[clusterIndex1D] Z0, Y0, X0 = clusterIndex1D #### delete dataCLuster to recollect the memory del dataCluster gc.collect() #### startLindexCol="startLindex" startBindexCol="startBindex" startVindexCol="startVindex" cloudTB[startLindexCol]= cloudTB["peakL"]*0 cloudTB[startBindexCol]= cloudTB["peakL"]*0 cloudTB[startVindexCol]= cloudTB["peakL"]*0 endLindexCol="endLindex" endBindexCol="endBindex" endVindexCol="endVindex" cloudTB[endLindexCol] = cloudTB["peakL"] * 0 cloudTB[endBindexCol] = cloudTB["peakL"] * 0 cloudTB[endVindexCol] = cloudTB["peakL"] * 0 startLCol="startL" startBCol="startB" startVCol="startV" cloudTB[startLCol] = cloudTB["x_cen"] * 0 cloudTB[startBCol] = cloudTB["y_cen"] * 0 cloudTB[startVCol] = cloudTB["v_cen"] * 0 endLCol="endL" endBCol="endB" endVCol="endV" cloudTB[endLCol] = cloudTB["x_cen"] * 0 cloudTB[endBCol] = cloudTB["y_cen"] * 0 cloudTB[endVCol] = cloudTB["v_cen"] * 0 #twelve colnames # # print cloudTB.colnames for eachC in cloudTB: cloudID = eachC["_idx"] cloudIndex = self.getIndices(Z0, Y0, X0, clusterValue1D, cloudID) cloudZ0, cloudY0, cloudX0 = cloudIndex eachC[ startVindexCol ] = min(cloudZ0) eachC[ endVindexCol ] = max(cloudZ0) eachC[ startBindexCol ] = min(cloudY0) eachC[ endBindexCol ] = max(cloudY0) eachC[ startLindexCol ] = min(cloudX0) eachC[ endLindexCol ] = max(cloudX0) l0,b0,v0=wcsCloud.wcs_pix2world(eachC[ endLindexCol ],eachC[ startBindexCol ] ,eachC[ startVindexCol ] , 0 ) eachC[startLCol] = l0 eachC[startBCol] = b0 eachC[startVCol] = v0/1000. l1,b1,v1=wcsCloud.wcs_pix2world(eachC[ startLindexCol ],eachC[ endBindexCol ] ,eachC[ endVindexCol ] , 0 ) eachC[endLCol] = l1 eachC[endBCol] = b1 eachC[endVCol] = v1/1000. cloudTB.write(catNameWithLBVrange,overwrite=True) def getCloudCubeFITS(self, cloudID, savePath=None): """ """ if savePath is None: savePath = os.path.join(self.processPath, self.cloudCubeSaveFolder ) saveName = "cloud{}cube_{}.fits".format(cloudID ,self.extraTag) return os.path.join(savePath, saveName ) def produceIndividualClouds(self, savePath=None, noiseMask=0, inputTB=None,cropNoise=False ): """ Revisement: By qzyan for adding a noise cube for each individual clouds #output all data cubes for each cloud :return: """ ################# # savePath="" self.getOutputNames() if savePath is None: savePath = os.path.join(self.processPath, self.cloudCubeSaveFolder ) if os.path.isdir(savePath): pass else: os.makedirs(savePath) ##### TB=Table.read(self.cleanCatName) if inputTB is not None : cloudTB = inputTB else: TB.sort("area_exact") TB.reverse() cloudTB = TB[TB["area_exact"] >= self.areaThresh * 3600] ## Remove small molecular clouds print("Removing clouds with area less than {} deg".format(self.areaThresh)) print("Producing Individual clouds for {} clouds".format(len(cloudTB))) dataCluster, headCluster = myFITS.readFITS(self.cleanLabelFITSName) dataCO, headCO = myFITS.readFITS(self.rawCOFITS) # print cloudTB minV = np.nanmin(dataCluster[0]) wcsCloud = WCS(headCluster, naxis=3) clusterIndex1D = np.where(dataCluster > minV) clusterValue1D = dataCluster[clusterIndex1D] Z0, Y0, X0 = clusterIndex1D #### delete dataCLuster to recollect the memory del dataCluster gc.collect() #### fitsZero = np.zeros_like(dataCO, dtype=np.float32) fitsZero = fitsZero + noiseMask if cropNoise: noiseMaskZero = np.zeros_like(dataCO, dtype=np.int) noiseData= dataCO.copy() noiseData[clusterIndex1D] = np.nan #all singla regions are masked #yxProj=np.zeros_like(dataCO[0,:,:],dtype=np.int64 ) #zxProj=np.zeros_like(dataCO[:,0,:],dtype=np.int64 ) #zyProj=np.zeros_like(dataCO[:,:,0],dtype=np.int64 ) # print cloudTB.colnames for eachC in cloudTB: cloudID = eachC["_idx"] saveName = self.getCloudCubeFITS(cloudID,savePath=savePath) #"cloud{}cube.fits".format(cloudID) saveNameNoise = "cloud{}cube_noise.fits".format(cloudID) cloudIndex = self.getIndices(Z0, Y0, X0, clusterValue1D, cloudID) fitsZero[cloudIndex] = dataCO[cloudIndex] cloudZ0, cloudY0, cloudX0 = cloudIndex minZ = np.min(cloudZ0) maxZ = np.max(cloudZ0) minY = np.min(cloudY0) maxY = np.max(cloudY0) minX = np.min(cloudX0) maxX = np.max(cloudX0) cropWCS = wcsCloud[minZ:maxZ + 1, minY:maxY + 1, minX:maxX + 1] cropData = fitsZero[minZ:maxZ + 1, minY:maxY + 1, minX:maxX + 1] fits.writeto(saveName , cropData, header=cropWCS.to_header(), overwrite=True) fitsZero[:] = 0 continue fitsZero[cloudIndex] = 1 # rest to maskfits #### deal Z axis xyIndex= (Y0, X0) yxProj[xyIndex] = 1 spectraIndicesYX = np.where(yxProj==1) for i in range( len(spectraIndicesYX[0] ) ): spectraYindex = spectraIndicesYX[0][i] spectraXindex = spectraIndicesYX[1][i] spectraMask = fitsZero[:, spectraYindex, spectraXindex ] expandMaskV = self.expandMask( spectraMask ) if cropNoise: noiseMaskZero[ : , spectraYindex, spectraXindex ] = noiseMaskZero[ : , spectraYindex, spectraXindex ] + expandMaskV #### deal Y axis zxIndex= (Z0, X0) zxProj[zxIndex] = 1 spectraIndicesZX = np.where( zxProj==1 ) for i in range( len(spectraIndicesZX[0] )): spectraZindex = spectraIndicesZX[0][i] spectraXindex = spectraIndicesZX[1][i] spectraMask = fitsZero[spectraZindex, :, spectraXindex ] expandMaskV = self.expandMask( spectraMask ) if cropNoise: noiseMaskZero[ spectraZindex, :, spectraXindex ] = noiseMaskZero[ spectraZindex, :, spectraXindex ] + expandMaskV #### deal X axis zyIndex= (Z0, Y0) zyProj[zyIndex] = 1 spectraIndicesZY = np.where( zyProj==1 ) for i in range( len( spectraIndicesZY[0] )): spectraZindex = spectraIndicesZY[0][i] spectraYindex = spectraIndicesZY[1][i] spectraMask = fitsZero[spectraZindex, spectraYindex, : ] expandMaskV = self.expandMask( spectraMask ) if cropNoise: noiseMaskZero[ spectraZindex, spectraYindex, : ] = noiseMaskZero[ spectraZindex, spectraYindex, : ] + expandMaskV if cropNoise: ###crop the expanded mask regions noiseMaskZero[cloudIndex] =0 noiseIndices = np.where(noiseMaskZero>0) noiseZ0, noiseY0, noiseX0 = noiseIndices minZnoise = np.min(noiseZ0) maxZnoise = np.max(noiseZ0) minYnoise = np.min(noiseY0) maxYnoise = np.max(noiseY0) minXnoise = np.min(noiseX0) maxXnoise = np.max(noiseX0) cropWCSnoise = wcsCloud[minZnoise:maxZnoise + 1, minYnoise:maxYnoise + 1, minXnoise:maxXnoise + 1] cropDataNoise = noiseData[minZnoise:maxZnoise + 1, minYnoise:maxYnoise + 1, minXnoise:maxXnoise + 1] cropDataNoise= cropDataNoise.copy() cropNoiseMask = noiseMaskZero[minZnoise:maxZnoise + 1, minYnoise:maxYnoise + 1, minXnoise:maxXnoise + 1] cropDataNoise[cropNoiseMask==0] = np.nan fits.writeto(os.path.join(savePath,saveNameNoise), cropDataNoise, header=cropWCSnoise.to_header(), overwrite=True) del cropDataNoise gc.collect( ) ### initialize noiseMaskZero[:] = 0 fitsZero[:] = 0 yxProj[:] = 0 zxProj[:] = 0 zyProj[:] = 0 print("Cloud fits writing done!") def expandMask(self, spectraMask ): """ extend the range of mask, the mask is actually the mask of single, expand this maks to close regions :param range: :param size: :return: """ signalMask = np.ma.masked_array(spectraMask,mask=spectraMask) #masked are regions length = len(spectraMask) signalRegions = np.ma.clump_masked( signalMask ) for eachRegion in signalRegions: startI = eachRegion.start endI= eachRegion.stop sizeRegion = endI-startI ##### back mask backEnd= min([endI+sizeRegion,length-1] ) signalMask[endI:backEnd]=np.ma.masked #### fore mask foreStar=max([0, startI-sizeRegion] ) signalMask[foreStar:startI]=np.ma.masked return signalMask.mask*1 def getCloudIDByRow(self, eachC): if "_idx" in eachC.colnames: return eachC["_idx"] else: cloudName = eachC["sourceName"] id = cloudName.split("Cloud")[-1] return int(id) def removeKM34(self): """ Remove bad channels around 34 km/s :return: """ #### #### particularly for cloud in the third Galactic quadrant #### #### #### # this is only a rough ###alias def pipeLine(self, redoDBSCAN=True, onlyCat=False, writeconfig=False, splitN=1,getIntFITS=False,getCubeFITS=False): """ :param rawCOFITS: :param redoDBSCAN: :param outPath: :param rmsFITS: :param rmsMean: :param splitN: to overcome the memory problem :param getIntFITS: :param getCubeFITS: :return: """ self.initDBSCAN( ) if writeconfig: self.writeConfig() #do we need to write down the case? if redoDBSCAN: print('Redo DBSCAN') self.computeDBSCAN(splitN= splitN ) self.getDBSCANCatlog(doClean=self.doClean) self.produceCleanFITS() if onlyCat: self.getDBSCANCatlog(doClean=self.doClean) self.produceCleanFITS() if not redoDBSCAN and not onlyCat: print('Recalling old DBSCAN') self.getOutputNames( ) if getIntFITS: self.produceCloudIntFITS( ) if getCubeFITS: self.produceIndividualClouds() def getLrangeList(self, lRange, splitN, overLap): """ get a list of lRanges :return: """ lRangeEdges = np.linspace(min(lRange), max(lRange), splitN + 1) N = len(lRangeEdges) lList = [] for i in range(N - 1): if i == 0: lList.append([lRangeEdges[i], lRangeEdges[i + 1] + overLap / 2.]) continue if i == N - 2: lList.append([lRangeEdges[i] - overLap / 2., lRangeEdges[i + 1]]) continue lList.append([lRangeEdges[i] - overLap / 2., lRangeEdges[i + 1] + overLap / 2.]) return lList def mosaicLabel(self, rawCOFITS, rmsFITS, processPath, splitN=3, overlapL=1): """ #split the rawCOFITS into many smaller regions, do DBSCAN speratelay and then merge then togher, need to star from left to right, because the index in the right is used :param rawCOFITS: :param rmsFITS: :return: """ lRange, Brange, Vrange = doFITS.getLBVRange(rawCOFITS) lList = self.getLrangeList(lRange, splitN, overlapL) lList = lList[::-1] # from left to rights #### ### #### # crop the fits along the velocity cropListRaw = [] rmsList = [] #### baseNameRaw = os.path.basename(rawCOFITS) baseNameRMS = os.path.basename(rmsFITS) for i in range(splitN): cropListRaw.append(os.path.join(processPath, "split_{}_".format(i) + baseNameRaw)) rmsList.append(os.path.join(processPath, "split_{}_".format(i) + baseNameRMS)) for i in range(splitN): lRnage = lList[i] doFITS.cropFITS(rawCOFITS, Lrange=lRnage, outFITS=cropListRaw[i], overWrite=True) doFITS.cropFITS2D(rmsFITS, Lrange=lRnage, outFITS=rmsList[i], overWrite=True) #### # do dbscan for each labelFITSList = [] coreFITSList = [] for i in range(splitN): self.rawCOFITS = cropListRaw[i] self.rmsFITS = rmsList[i] self.processPath = processPath self.saveCore = True labelFITSi = self.computeDBSCAN() labelFITSList.append(labelFITSi) coreFITSList.append(self.coreFITS) leftFITSRaw = labelFITSList[0] leftFITSCore = coreFITSList[0] for i in range(splitN - 1): # merge if i == splitN - 1 - 1: mergei = os.path.join(processPath, "merge_All_rawlabel.fits") mergeCore = os.path.join(processPath, "merge_All_corelabel.fits") else: mergei = os.path.join(processPath, "merge_{}_rawlabel.fits".format(i)) mergeCore = os.path.join(processPath, "merge_{}_corelabel.fits".format(i)) self.getMergeData(leftFITSRaw, labelFITSList[i + 1], saveName=mergei) leftFITSRaw = mergei self.getMergeData(leftFITSCore, coreFITSList[i + 1], saveName=mergeCore) leftFITSCore = mergeCore # self.mergeTwoLabelFITS( leftFITS , labelFITSList[i+1], outFITSName = mergei ) #### s = generate_binary_structure(3, self.connectivity) labelFITSData, headlabel = doFITS.readFITS(mergei) dataCore, headCore = doFITS.readFITS(mergeCore) dataCore = dataCore > 0 gc.collect() labeled_core, num_features = label(dataCore, structure=s) gc.collect() labeled_core = dilation(labeled_core, s) gc.collect() labeled_core[labelFITSData == 0] = 0 gc.collect() fits.writeto("mergeTestFinale.fits", labeled_core, header=headlabel, overwrite=True) def mergeTwoLabelFITS(self, fits1, fits2, outFITSName=None): """ :return: mergeTwo labe FITS, by default, the overlapping area is along the Galactic plane # the left fits is the main fits """ print(fits1, fits2) #### getLBVRange l1, b1, v1 = doFITS.getLBVRange(fits1) l2, b2, v2 = doFITS.getLBVRange(fits2) if np.max(l1) < np.max(l2): fits1, fits2 = fits2, fits1 l1, b1, v1 = doFITS.getLBVRange(fits1) l2, b2, v2 = doFITS.getLBVRange(fits2) ######## breakL = np.mean([min(l1), max(l2)]) dL = abs(min(l1) - max(l2)) mergeLRange = [breakL + dL / 4, breakL - dL / 4] # check if the data in the range are equal fits1Crop = "./tmp/merge1Crop.fits" fits2Crop = "./tmp/merge2Crop.fits" doFITS.cropFITS(fits1, Lrange=mergeLRange, outFITS=fits1Crop, overWrite=True) doFITS.cropFITS(fits2, Lrange=mergeLRange, outFITS=fits2Crop, overWrite=True) data1, head1 = doFITS.readFITS(fits1Crop) data2, head2 = doFITS.readFITS(fits2Crop) maskData1 = data1 > 0 maskData1 = maskData1 + 0 maskData2 = data2 > 0 maskData2 = maskData2 + 0 if not np.sum(maskData2 - maskData1) == 0: print("DBSCAN results between common area are not equal, exist") ## start to relabel the fits Part1 = "./tmp/part1Merge.fits" Part2 = "./tmp/part2Merge.fits" doFITS.cropFITS(fits1, Lrange=[max(l1), breakL - dL / 4], outFITS=Part1, overWrite=True) doFITS.cropFITS(fits2, Lrange=[min(l2), breakL + dL / 4], outFITS=Part2, overWrite=True) dataPart1, headPart1 = doFITS.readFITS(Part1) dataPart2, headPart2 = doFITS.readFITS(Part2) ### # data1,data2, should be the same #### maxPart1 = np.max(dataPart1) dataPart2 = dataPart2 + maxPart1 dataPart2[dataPart2 == maxPart1] = 0 data2 = data2 + maxPart1 data2[data2 == maxPart1] = 0 #### # #### fruther modify part2 indicesOverlap1 = np.where(data1 > 0) overlappingData1 = data1[indicesOverlap1] overLappingCloud1 = np.unique(overlappingData1) overZ1, overY1, overX1 = indicesOverlap1 indicesOverlap2 = np.where(data2 > 0) overlappingData2 = data2[indicesOverlap2] overLappingCloud2 = np.unique(overlappingData2) overZ2, overY2, overX2 = indicesOverlap2 indicesPart2 = np.where(dataPart2 > 0) part2Z, part2Y, part2X = indicesPart2 part2Values = dataPart2[indicesPart2] #### clould be very large.... indicesPart1 = np.where(dataPart1 > 0) part1Z, part1Y, part1X = indicesPart1 part1Values = dataPart1[indicesPart1] #### first, check multiplicity for eachMC2 in overLappingCloud2: # index in mergeData overlapPart2Indices = self.getIndices(overZ2, overY2, overX2, overlappingData2, eachMC2) overlap1Values = data1[overlapPart2Indices] if len(overlap1Values) >= 2: # two values print(overlap1Values, " corresponds ", eachMC2, "in data 2") # need to rewrite the ########## for eachMC in overLappingCloud1: # index in mergeData part1Indices = self.getIndices(overZ1, overY1, overX1, overlappingData1, eachMC) MC12 = data2[part1Indices] # to check if the values are the same correspondMCs = np.unique(MC12) existingMCIndex = correspondMCs # np.unique(indexPart2) if 0 in existingMCIndex: print("Wrong, uncompatible fits") for cloudParat2 in existingMCIndex: indicesDataPart2 = self.getIndices(part2Z, part2Y, part2X, part2Values, cloudParat2) dataPart2[indicesDataPart2] = eachMC fits.writeto(Part2, dataPart2, header=headPart2, overwrite=True) ################ #### self.getMergeData(Part1, Part2, breakL, saveName=outFITSName) ###### # self.labelFITSName= labelFITS # self.getCatFromLabelArray(doClean= True ) def getMergeData(self, fits1, fits2, breakL=None, saveName=None): #### """ :param fits1: :param fits2: :param breakL: :return: """ l1, b1, v1 = doFITS.getLBVRange(fits1) l2, b2, v2 = doFITS.getLBVRange(fits2) if np.max(l1) < np.max(l2): fits1, fits2 = fits2, fits1 l1, b1, v1 = doFITS.getLBVRange(fits1) l2, b2, v2 = doFITS.getLBVRange(fits2) ######## breakL = np.mean([min(l1), max(l2)]) Part1 = "./tmp/part1ForMerging.fits" Part2 = "./tmp/part2ForMerging.fits" doFITS.cropFITS(fits1, Lrange=[max(l1), breakL], outFITS=Part1, overWrite=True) doFITS.cropFITS(fits2, Lrange=[min(l2), breakL], outFITS=Part2, overWrite=True) # one column is overlapping, removing it dataPart1, head1 = doFITS.readFITS(Part1) dataPart2, head2 = doFITS.readFITS(Part2) maxPart1 = np.max(dataPart1) dataPart2 = dataPart2 + maxPart1 dataPart2[dataPart2 == maxPart1] = 0 Nz1, Ny1, Nx1 = dataPart1.shape wcs1 = WCS(head1, naxis=2) wcs2 = WCS(head2, naxis=2) l1Last, _ = wcs1.wcs_pix2world(Nx1 - 1, Ny1 / 2, 0) startPart2Index, _ = wcs2.wcs_world2pix(l1Last, np.mean(b2), 0) if startPart2Index < 0: print("---wrong...existing...") return startPart2Index = int(round(startPart2Index)) # cutposition of dataPart2 dataPart2 = dataPart2[:, :, startPart2Index + 1:] newData = np.concatenate((dataPart1, dataPart2), axis=2) # newData[:,:,Nx1:]=0 # return newData if saveName is None: saveName = "mergedFITS.fits" # # s = generate_binary_structure(3, self.connectivity ) # newData[newData>0] =1 # labeled_core, num_features = label(newData, structure=s) # first label core, then expand, otherwise, the expanding would wrongly connected fits.writeto(saveName, newData, header=head1, overwrite=True) return saveName def setConfigFile(self,configFile): """ :param configFile: :return: """ print("checking ", configFile ) if configFile is not None and os.path.isfile( configFile ): self.configFile=configFile return #else pathConfigure,fileConfigure=os.path.split(configFile) if os.path.isdir(pathConfigure): ####create the files with open(configFile,'w') as file: pass self.configFile = configFile #else,test def getDBSCANcase(self,regionName,configFile=None ): if configFile is not None: self.setConfigFile(configFile) self.set_region(regionName) #### def ZZZ(self): pass ### alisas initDBSCAN = getOutputNames pipeline=pipeLine """ Use example from mwispDBSCANV5dev import MWISPDBSCAN doDBSCAN = MWISPDBSCAN() doDBSCAN.setConfigFile(None) doDBSCAN.extraTag="" doDBSCAN.rawCOFITS="" doDBSCAN.rmsFITS="" doDBSCAN.pipeLine(writeconfig=True ) """ ### recall DBSCAN results """ Use example from mwispDBSCANV5dev import MWISPDBSCAN doDBSCAN = MWISPDBSCAN() doDBSCAN.setConfigFile(None) doDBSCAN.set_region() """