diff --git a/.gitignore b/.gitignore index ef25bf1..5558c1a 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,11 @@ *.pyc *.o *.so +venv +*.txt +*.csv +*.obj +__pycache__ +.ipynb_checkpoints/ +.radiopadre +.radiopadre-session diff --git a/ClassDynSpecMS.py b/ClassDynSpecMS.py index 7797ff7..1ffd72c 100644 --- a/ClassDynSpecMS.py +++ b/ClassDynSpecMS.py @@ -10,7 +10,7 @@ log=logger.getLogger("DynSpecMS") from DDFacet.Array import shared_dict from DDFacet.Other import AsyncProcessPool - + from DDFacet.Other import Multiprocessing from DDFacet.Other import ModColor from DDFacet.Other.progressbar import ProgressBar @@ -21,7 +21,7 @@ import os from killMS.Other import reformat from DDFacet.Other import AsyncProcessPool -from .dynspecms_version import version +from dynspecms_version import version import glob from astropy.io import fits from astropy.wcs import WCS @@ -35,6 +35,7 @@ #from killMS.Data import ClassJonesDomains import DDFacet.Other.ClassJonesDomains import psutil +import ClassGiveCatalog def AngDist(ra0,ra1,dec0,dec1): AC=np.arccos @@ -67,8 +68,8 @@ def __init__(self, ColName="DATA", TChunkHours=0., ModelName="PREDICT_KMS", - UVRange=[1.,1000.], - ColWeights=None, + UVRange=[1.,1000.], + ColWeights=None, SolsName=None, FileCoords="Transient_LOTTS.csv", Radius=3., @@ -82,8 +83,11 @@ def __init__(self, SourceCatOff=None, SourceCatOff_FluxMean=None, SourceCatOff_dFluxMean=None, - options=None): + options=None,SubSet=None): + self.options=options + self.SubSet=SubSet + # CutGainsMinMax if BeamModel=="None": @@ -112,8 +116,8 @@ def __init__(self, if ListMSName is None: print(ModColor.Str("WORKING IN REPLOT MODE"), file=log) self.Mode="Plot" - - + + self.Radius=Radius self.ImageI = ImageI self.ImageV = ImageV @@ -123,14 +127,14 @@ def __init__(self, # identify version in logs print("DynSpecMS version %s starting up" % version(), file=log) self.FileCoords=FileCoords - + if self.Mode=="Spec": self.ListMSName = sorted(ListMSName)#[0:2] self.nMS = len(self.ListMSName) self.OutName = self.ListMSName[0].split("/")[-1].split("_")[0] self.ReadMSInfos() self.InitFromCatalog() - + elif self.Mode=="Plot": self.OutName = self.BaseDirSpecs.split("_")[-1] self.InitFromSpecs() @@ -150,10 +154,10 @@ def InitFromSpecs(self): self.NDirSelected=len(ListTargetFits) NOrig=len(ListTargetFits) - + ListOffFits=glob.glob("%s/OFF/*.fits"%self.BaseDirSpecs) NOff=len(ListOffFits) - + self.PosArray=np.zeros((self.NDirSelected+NOff,),dtype=[('Name','S200'),("ra",np.float64), ("dec",np.float64),('Type','S200'), ('iTessel',int),('iFacet',int)]) @@ -172,7 +176,7 @@ def InitFromSpecs(self): if ra<0.: ra+=2.*np.pi self.PosArray.ra[iDir]=ra dec=self.PosArray.dec[iDir]=float(F[0].header['DEC_RAD']) - + # print File,rad2hmsdms(ra,Type="ra").replace(" ",":"),rad2hmsdms(dec,Type="dec").replace(" ",":") # if self.PosArray.Type[iDir]=="Off": stop for iPol in range(4): @@ -195,187 +199,18 @@ def InitFromSpecs(self): continue self.PosArray.Type[iDir]=PosArrayTarget.Type[iS] self.PosArray.Name[iDir]=PosArrayTarget.Name[iS] - + def InitFromCatalog(self): FileCoords=self.FileCoords - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] # should we use the surveys DB? - DoProperMotionCorr=False - if self.options.UseLoTSSDB: - print("Using the surveys database", file=log) - from surveys_db import SurveysDB - with SurveysDB() as sdb: - sdb.cur.execute('select * from transients') - result=sdb.cur.fetchall() - # convert to a list, then to ndarray, then to recarray - l=[] - for r in result: - l.append((r['id'],r['ra'],r['decl'],r['type'])) - - if FileCoords is not None and FileCoords!="": - print('Adding data from file '+FileCoords, file=log) - additional=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] - if len(additional.shape)==0: additional=additional.reshape((1,)) - if not additional.shape: - # deal with a one-line input file - additional=np.array([additional],dtype=dtype) - for r in additional: - l.append(tuple(r)) - - self.PosArray=np.asarray(l,dtype=dtype) - print("Created an array with %i records" % len(result), file=log) - elif self.options.UseGaiaDB is not None: - from astroquery.gaia import Gaia - rac_deg,decc_deg=self.ra0*180/np.pi, self.dec0*180/np.pi - Radius_deg=self.Radius - Dmax,NMax=self.options.UseGaiaDB.split(",") - Dmax,NMax=float(Dmax),int(NMax) - Parallax_min=1./(Dmax*1e-3) - query=f"""SELECT TOP 10000 gaia_source.designation,gaia_source.source_id,gaia_source.ref_epoch,gaia_source.ra,gaia_source.dec,gaia_source.parallax, - gaia_source.pmra,gaia_source.pmdec,gaia_source.phot_g_mean_mag,gaia_source.bp_rp - FROM gaiadr3.gaia_source - WHERE - CONTAINS( - POINT('ICRS',gaiadr3.gaia_source.ra,gaiadr3.gaia_source.dec), - CIRCLE('ICRS',{rac_deg},{decc_deg},{Radius_deg}) - )=1 AND (gaiadr3.gaia_source.parallax>={Parallax_min})""" - log.print(f"Sending du Gaia server") - print(f"{query}") - job = Gaia.launch_job(query) - result = job.get_results() - - log.print(f"Query has returned {len(result)}") - - - - l=[] - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), - ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64), - ("parallax",np.float64),("GaiaDistance",np.float64), - ("g",np.float64),("G",np.float64),("b_r",np.float64), - ('Type','S200')] - for r in result: - parallax=r['parallax'] - g=r['g'] - D=1./(parallax*1e-3) - G=g+5*np.log10(parallax)-10 - l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'], - r['parallax'],D, - g,G,r['b_r'], - b"Gaia DR3")) - self.PosArray=np.asarray(l,dtype=dtype) - - - CGGS=ClassGiveGaiaSample((rac_deg,decc_deg,Radius_deg),self.PosArray,RefCat="/data/cyril.tasse/Analyse_DataDynSpec_Jan23_TestRM/MergedCat.npz.MergeGaia.npz") - indExo,indGaia=CGGS.buildRandGaiaSample() - stop - # ################################## - F=fits.open(self.options.FitsCatalog) - d=F[1].data - l=[] - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), - ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), - ('Type','S200')] - for r in d: - l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],r['Type'])) - - if FileCoords is not None and FileCoords!="": - print('Adding data from file '+FileCoords, file=log) - dtype0=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] - additional=np.genfromtxt(FileCoords,dtype=dtype0,delimiter=",")[()] - if len(additional.shape)==0: additional=additional.reshape((1,)) - if not additional.shape: - # deal with a one-line input file - additional=np.array([additional],dtype=dtype0) - additional1=np.zeros((additional.shape[0],),dtype=dtype) - additional1["Name"][:]=additional["Name"][:] - additional1["ra"][:]=additional["ra"][:] - additional1["dec"][:]=additional["dec"][:] - additional1["Type"][:]=additional["Type"][:] - - for r in additional1: - l.append(tuple(r)) - - - self.PosArray=np.asarray(l,dtype=dtype) - log.print("Created an array with %i records" % len(l)) - DoProperMotionCorr=True - - - # import pylab - # pylab.clf() - # pylab.subplot(1,2,1) - # pylab.scatter(self.PosArray["ra"],self.PosArray["dec"]) - # pylab.scatter([rac_deg],[decc_deg]) - # pylab.subplot(1,2,2) - # D=1./(self.PosArray["parallax"]*1e-3) - # pylab.hist(D,bins=100) - # pylab.draw() - # pylab.show() - DoProperMotionCorr=True - if self.PosArray.size>NMax: - ind=np.int64(np.random.rand(NMax)*self.PosArray.size) - self.PosArray=self.PosArray[ind] - - log.print("Created an array with %i records" % self.PosArray.size) - - elif self.options.FitsCatalog: - print("Using the fits catalog: %s"%self.options.FitsCatalog, file=log) - F=fits.open(self.options.FitsCatalog) - d=F[1].data - l=[] - dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), - ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), - ('Type','S200')] - for r in d: - l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],r['Type'])) - - if FileCoords is not None and FileCoords!="": - print('Adding data from file '+FileCoords, file=log) - dtype0=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] - additional=np.genfromtxt(FileCoords,dtype=dtype0,delimiter=",")[()] - if len(additional.shape)==0: additional=additional.reshape((1,)) - if not additional.shape: - # deal with a one-line input file - additional=np.array([additional],dtype=dtype0) - additional1=np.zeros((additional.shape[0],),dtype=dtype) - additional1["Name"][:]=additional["Name"][:] - additional1["ra"][:]=additional["ra"][:] - additional1["dec"][:]=additional["dec"][:] - additional1["Type"][:]=additional["Type"][:] - - for r in additional1: - l.append(tuple(r)) - - - self.PosArray=np.asarray(l,dtype=dtype) - log.print("Created an array with %i records" % len(l)) - DoProperMotionCorr=True - else: - - #FileCoords="Transient_LOTTS.csv" - if FileCoords is None: - if not os.path.isfile(FileCoords): - ssExec="wget -q --user=anonymous ftp://ftp.strw.leidenuniv.nl/pub/tasse/%s -O %s"%(FileCoords,FileCoords) - print("Downloading %s"%FileCoords, file=log) - print(" Executing: %s"%ssExec, file=log) - os.system(ssExec) - log.print("Reading cvs file: %s"%FileCoords) - #self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] - self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",") - if len(self.PosArray.shape)==0: self.PosArray=self.PosArray.reshape((1,)) - - self.PosArray=self.PosArray.view(np.recarray) - self.PosArray.ra*=np.pi/180. - self.PosArray.dec*=np.pi/180. - - Radius=self.Radius - NOrig=self.PosArray.Name.shape[0] - Dist=AngDist(self.ra0,self.PosArray.ra,self.dec0,self.PosArray.dec) - ind=np.where(Dist<(Radius*np.pi/180))[0] - self.PosArray=self.PosArray[ind] + CGC=ClassGiveCatalog.ClassGiveCatalog(self.options, + self.ra0,self.dec0, + self.Radius, + self.FileCoords) + self.PosArray=CGC.giveCat(SubSet=self.SubSet) + DoProperMotionCorr=CGC.DoProperMotionCorr if DoProperMotionCorr: Lra,Ldec=[],[] @@ -383,31 +218,31 @@ def InitFromCatalog(self): for iPCat,PCat in enumerate(self.PosArray): #print(iPCat,len(self.PosArray)) ra,dec=PCat["ra"],PCat["dec"] - + ra1,dec1=ProperMotionCorrection(PCat["ra"],PCat["dec"], PCat["pmra"],PCat["pmdec"],PCat["ref_epoch"],PCat["parallax"],self.tmin) if np.isnan(ra1) or np.isnan(dec1): ra1,dec1=ra,dec log.print(str((PCat["ra"],PCat["dec"],PCat["pmra"],PCat["pmdec"],PCat["ref_epoch"],PCat["parallax"],ra1,dec1))) - + Lra.append(ra1-ra) Ldec.append(dec1-dec) self.PosArray["ra"][iPCat]=ra1 self.PosArray["dec"][iPCat]=dec1 - print(np.array(Lra)*3600.*180/np.pi,np.array(Ldec)*3600.*180/np.pi) + # print(np.array(Lra)*3600.*180/np.pi,np.array(Ldec)*3600.*180/np.pi) + - self.NDirSelected=self.PosArray.shape[0] - print("Selected %i target [out of the %i in the original list]"%(self.NDirSelected,NOrig), file=log) + print("Selected %i target [out of the %i in the original list]"%(self.NDirSelected,CGC.NOrig), file=log) if self.NDirSelected==0: print(ModColor.Str(" Have found no sources - returning"), file=log) self.killWorkers() return - + NOff=self.NOff - + if NOff==-1: NOff=self.PosArray.shape[0]*2 if NOff is not None: @@ -421,7 +256,7 @@ def InitFromCatalog(self): theta=np.linspace(0,2*np.pi,10) xC0=np.cos(theta)*Rrad0 yC0=np.sin(theta)*Rrad0 - + def give_iFacet_iTessel(l,m): #Plm=Polygon.Polygon(np.array([[l,m]])) Plm=Polygon.Polygon(np.array([xC0+l,yC0+m]).T) @@ -434,11 +269,11 @@ def give_iFacet_iTessel(l,m): ArrMax=Arr iFacetMax=iFacet return iFacetMax,self.DFacet[iFacetMax]["iSol"] - + self.PosArray=RecArrayOps.AppendField(self.PosArray,'iFacet',int) self.PosArray=RecArrayOps.AppendField(self.PosArray,'iTessel',int) l, m = self.CoordMachine.radec2lm(self.PosArray.ra, self.PosArray.dec) - + for iThis in range(self.PosArray.size): iFacet,iSol=give_iFacet_iTessel(l[iThis],m[iThis]) self.PosArray.iFacet[iThis]=iFacet @@ -451,7 +286,7 @@ def give_iFacet_iTessel(l,m): # pylab.scatter(self.PosArray.ra[ind],self.PosArray.dec[ind])#,c=iTessel) # pylab.draw() # pylab.show() - + self.NDir=self.PosArray.shape[0] print("For a total of %i targets"%(self.NDir), file=log) @@ -466,7 +301,7 @@ def give_iFacet_iTessel(l,m): self.DicoGrids["GridWeight2"] = np.zeros((self.NDir,self.NChan, self.NTimesGrid, 4), np.complex128) - + self.DoJonesCorr_kMS =False self.DicoJones=None @@ -477,7 +312,7 @@ def give_iFacet_iTessel(l,m): if self.BeamModel: self.DoJonesCorr_Beam=True - AsyncProcessPool.APP=None + #AsyncProcessPool.APP=None # AsyncProcessPool.init(ncpu=self.NCPU, # num_io_processes=1, # affinity="disable") @@ -486,14 +321,14 @@ def give_iFacet_iTessel(l,m): num_io_processes=1, verbose=0) self.APP=AsyncProcessPool.APP - + self.APP.registerJobHandlers(self) self.APP.startWorkers() - - - + + + def GiveOffPosArray(self,NOff): print("Making random off catalog with %i directions"%NOff, file=log) CatOff=np.zeros((NOff,),self.PosArray.dtype) @@ -512,13 +347,13 @@ def GiveOffPosArray(self,NOff): ind=np.where( (Fd.Isl_Total_flux>S0) & (Fd.Isl_Total_flux1: stop + #if dtBin_.size>1: stop dtBin = dtBin_.flat[0] - - + + # if not np.allclose(ThisTimes, self.times): # raise ValueError("should have the same times") @@ -746,15 +581,15 @@ def ReadMSInfos(self): else: chSlice=slice(None) RevertChans=False - - + + if tmin is None: tmin=ThisTimes.min() tmax=ThisTimes.max() else: tmin=np.min([tmin,ThisTimes.min()]) tmax=np.max([tmax,ThisTimes.max()]) - + DicoMSInfos[iMS] = {"MSName": MSName, "ChanFreq": tf.getcol("CHAN_FREQ").ravel()[chSlice], # Hz "ChanWidth": np.abs(tf.getcol("CHAN_WIDTH").ravel()), # Hz @@ -767,21 +602,21 @@ def ReadMSInfos(self): "deltaTime": (ThisTimes[-1] - ThisTimes[0])/3600., # h "RevertChans": RevertChans, "Readable": True} - + if DicoMSInfos[iMS]["ChanWidth"][0] != self.ChanWidth: raise ValueError("should have the same chan width") pBAR.render(iMS+1, self.nMS) - + self.NTimesGrid=int(np.ceil((tmax-tmin)/dtBin)) self.timesGrid=tmin+np.arange(self.NTimesGrid)*dtBin self.tmin=tmin self.tmax=tmax - + for iMS in range(self.nMS): if not DicoMSInfos[iMS]["Readable"]: print(ModColor.Str("Problem reading %s"%MSName), file=log) print(ModColor.Str(" %s"%DicoMSInfos[iMS]["Exception"]), file=log) - + t.close() tf.close() @@ -795,12 +630,12 @@ def ReadMSInfos(self): #dt=np.median(dtArr) - + f0, f1 = self.Freq_minmax self.NChan = int((f1 - f0)/self.ChanWidth) + 1 # Fill properties - + self.tStart = Time(tmin/(24*3600.), format='mjd', scale='utc').isot self.tStop = Time(tmax/(24*3600.), format='mjd', scale='utc').isot self.fMin = self.Freq_minmax[0] @@ -825,24 +660,24 @@ def ReadMSInfos(self): LJob.append((iMS,iChunk)) self.LJob=LJob - + def LoadMS(self,iJob): iMS,iChunk=self.LJob[iJob] T0,T1=self.T0s[iChunk],self.T1s[iChunk] - - - if not self.DicoMSInfos[iMS]["Readable"]: + + + if not self.DicoMSInfos[iMS]["Readable"]: print("Skipping [%i/%i]: %s"%(iMS+1, self.nMS, self.ListMSName[iMS]), file=log) return "NotRead" print("Reading [%i/%i] %s:%s"%(iMS+1, self.nMS, self.ListMSName[iMS],self.ColName), file=log) MSName=self.ListMSName[iMS] - + t = table(MSName, ack=False) times = t.getcol("TIME") t0=self.tmin - + ind=np.where((times>=(T0+t0))&(times<(T1+t0)))[0] if ind.size==0: print("No Data in requested interval %f -> %f h time interval"%(T0/3600,T1/3600), file=log) @@ -852,17 +687,17 @@ def LoadMS(self,iJob): if ROW0!=0 or NROW!=t.nrows(): print(" Reading chunk in %.3f -> %.3f h"%(T0/3600,T1/3600), file=log) - + nch = self.DicoMSInfos[iMS]["ChanFreq"].size npol = self.DicoMSInfos[iMS]["npol"] #chSlice=self.DicoMSInfos[iMS]["chSlice"] RevertChans=self.DicoMSInfos[iMS]["RevertChans"] - + data = np.zeros((NROW,nch,npol),np.complex64) t.getcolnp(self.ColName,data,ROW0,NROW) if RevertChans: data=data[:,::-1,:] - + if self.ModelName: print(" Substracting %s from %s"%(self.ModelName,self.ColName), file=log) model=np.zeros((NROW,nch,npol),np.complex64) @@ -872,8 +707,8 @@ def LoadMS(self,iJob): data-=model del(model) - - + + if self.ColWeights: print(" Reading weight column %s"%(self.ColWeights), file=log) sW=t.getcol(self.ColWeights,0,1).shape @@ -884,8 +719,8 @@ def LoadMS(self,iJob): else: weights=np.zeros((NROW,nch),np.float32) t.getcolnp(self.ColWeights,weights,ROW0,NROW) - - + + if RevertChans: weights=weights[:,::-1] else: nrow,nch,_=data.shape @@ -894,13 +729,13 @@ def LoadMS(self,iJob): flag=np.zeros((NROW,nch,npol),np.bool) t.getcolnp("FLAG",flag,ROW0,NROW) if RevertChans: flag=flag[:,::-1] - + # data[:,:,:]=0 # data[:,0:300,0]=1 # flag.fill(0) # weights.fill(1) - + times = t.getcol("TIME",ROW0,NROW) A0, A1 = t.getcol("ANTENNA1",ROW0,NROW), t.getcol("ANTENNA2",ROW0,NROW) @@ -935,7 +770,7 @@ def LoadMS(self,iJob): DicoDATA["w"]=w0 DicoDATA["uniq_times"]=np.unique(DicoDATA["times"]) - + if self.DoJonesCorr_kMS or self.DoJonesCorr_Beam: self.setJones(DicoDATA) @@ -953,7 +788,7 @@ def setJones(self,DicoDATA): SolsName=SolsName.split(",") GD={"Beam":{"Model":self.BeamModel, "PhasedArrayMode":"A", - "At":"facet", + "At":"tessel", "DtBeamMin":5., "NBand":self.BeamNBand, "CenterNorm":1}, @@ -966,10 +801,27 @@ def setJones(self,DicoDATA): } print("Reading Jones matrices solution file:", file=log) + RADEC_ForcedBeamDirs=None + if not self.DoJonesCorr_kMS and self.DoJonesCorr_Beam: + RADEC_ForcedBeamDirs=(self.PosArray.ra,self.PosArray.dec) + + ms=ClassMS.ClassMS(self.DicoMSInfos[iMS]["MSName"],GD=GD,DoReadData=False,) JonesMachine = ClassJones.ClassJones(GD, ms, CacheMode=False) - JonesMachine.InitDDESols(DicoDATA) + JonesMachine.InitDDESols(DicoDATA,RADEC_ForcedBeamDirs=RADEC_ForcedBeamDirs) + + # print("================") + # RAJoneskillMS=DicoDATA["killMS"]["Dirs"]["ra"] + # DECJoneskillMS=DicoDATA["killMS"]["Dirs"]["dec"] + # for ThisRA,ThisDEC in zip(RAJoneskillMS,DECJoneskillMS): + # print(rad2hmsdms(ThisRA,Type="ra").replace(" ",":"),rad2hmsdms(ThisDEC,Type="dec").replace(" ",":")) + # print("================") + # RAJonesBeam=DicoDATA["Beam"]["Dirs"]["ra"] + # DECJonesBeam=DicoDATA["Beam"]["Dirs"]["dec"] + # for ThisRA,ThisDEC in zip(RAJonesBeam,DECJonesBeam): + # print(rad2hmsdms(ThisRA,Type="ra").replace(" ",":"),rad2hmsdms(ThisDEC,Type="dec").replace(" ",":")) + # print("================") CutGainsMinMax=StrToList(self.options.CutGainsMinMax) @@ -980,37 +832,46 @@ def setJones(self,DicoDATA): c0,c1=CutGainsMinMax G[G>c1]=0 G[G0)[0] #f = DicoDATA["flag"][indRow, :, :] #d = DicoDATA["data"][indRow, :, :] @@ -1189,18 +1054,18 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): d = np.array((DicoDATA["data"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol))).copy() f = np.array((DicoDATA["flag"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol))).copy() T.timeit("first") - + # for i in range(10): # d = (DicoDATA["data"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol)).copy() # f = (DicoDATA["flag"].flat[indArr.flat[:]]).reshape((nRowOut,nch,npol)).copy() # T.timeit("first %i"%i) - + nrow,nch,_=d.shape #weights = (DicoDATA["weights"][indRow, :]).reshape((nrow,nch,1)) - + indArr=nch*np.int64(indR)+np.int64(indCh) weights = np.array((DicoDATA["weights"].flat[indArr.flat[:]]).reshape((nRowOut,nch,1))).copy() - + A0s = DicoDATA["A0"][indRow].copy() A1s = DicoDATA["A1"][indRow].copy() u0 = DicoDATA["u"][indRow].reshape((-1,1,1)).copy() @@ -1214,14 +1079,14 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): # kk = np.exp( -2.*np.pi*1j* f/const.c.value *(u0*l + v0*m + w0*(n-1)) ) # Phasing term #print iTime,iDir ChanFreqs=np.array(self.DicoMSInfos[iMS]["ChanFreq"][0]).copy() - + iTimeGrid=np.argmin(np.abs(self.timesGrid-self.DicoMSInfos[iMS]["times"][iTime])) - + dcorr=d.copy() f0, _ = self.Freq_minmax ich0 = int( (ChanFreqs - f0)/self.ChanWidth ) OneMinusF=(1-f).copy() - + W=np.zeros((nRowOut,nch,npol),np.float32) for ipol in range(npol): W[:,:,ipol]=weights[:,:,0] @@ -1230,7 +1095,7 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): # weights=weights*np.ones((1,1,npol)) # W=weights - + kk=np.zeros_like(d) T.timeit("third") for iDir in range(self.NDir): @@ -1240,7 +1105,7 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): l, m = self.radec2lm(ra, dec,ra0,dec0) n = np.sqrt(1. - l**2. - m**2.) - + T.timeit("lmn") kkk = np.exp(-2.*np.pi*1j* chfreq/const.c.value *(u0*l + v0*m + w0*(n-1)) ) # Phasing term T.timeit("kkk") @@ -1258,21 +1123,22 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): # pylab.draw() # pylab.show(False) # pylab.pause(0.1) - - - + + + #DicoMSInfos = self.DicoMSInfos - + #_,nch,_=DicoDATA["data"].shape - + dcorr[:]=d[:] W=Wc.copy() + #W2=Wc.copy() dcorr*=W wdcorr=np.ones(dcorr.shape,np.float64) #kk=kk*np.ones((1,1,npol)) - + T.timeit("corr") - + if self.DoJonesCorr_kMS or self.DoJonesCorr_Beam: T1=ClassTimeIt.ClassTimeIt(" DoJonesCorr") T1.disable() @@ -1282,19 +1148,19 @@ def Stack_SingleTimeAllDir(self,iJob,iTime): tm = DicoJones['tm'] # Time slot for the solution iTJones=np.argmin(np.abs(tm-ThisTime))#self.timesGrid[iTime])) - + #iDJones=np.argmin(AngDist(ra,DicoJones['ra'],dec,DicoJones['dec'])) lJones, mJones = self.CoordMachine.radec2lm(DicoJones['ra'], DicoJones['dec']) iDJones=np.argmin(np.sqrt((l-lJones)**2+(m-mJones)**2)) - + _,nchJones,_,_,_,_=DicoJones['G'].shape T1.timeit("argmin") - - + + for iFJones in range(nchJones): - + nu0,nu1=DicoJones['FreqDomains'][iFJones] fData=self.DicoMSInfos[iMS]["ChanFreq"].ravel() indCh=np.where((fData>=nu0) & (fData1.]=1. + D[D<-1.]=-1. + else: + if D>1.: D=1. + if D<-1.: D=-1. + return AC(D) + +class ClassGiveCatalog(): + def __init__(self,options,ra0,dec0,Radius,FileCoords): + self.options=options + self.ra0=ra0 + self.dec0=dec0 + self.Radius=Radius + self.FileCoords=FileCoords + + + def giveCat(self,SubSet=None): + FileCoords=self.FileCoords + dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] + + self.DoProperMotionCorr=False + if self.options.UseLoTSSDB: + print("Using the surveys database", file=log) + from surveys_db import SurveysDB + with SurveysDB() as sdb: + sdb.cur.execute('select * from transients') + result=sdb.cur.fetchall() + # convert to a list, then to ndarray, then to recarray + l=[] + for r in result: + l.append((r['id'],r['ra'],r['decl'],r['type'])) + + if FileCoords is not None and FileCoords!="": + print('Adding data from file '+FileCoords, file=log) + additional=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] + if len(additional.shape)==0: additional=additional.reshape((1,)) + if not additional.shape: + # deal with a one-line input file + additional=np.array([additional],dtype=dtype) + for r in additional: + l.append(tuple(r)) + + self.PosArray=np.asarray(l,dtype=dtype) + elif self.options.UseGaiaDB is not None: + from astroquery.gaia import Gaia + rac_deg,decc_deg=self.ra0*180/np.pi, self.dec0*180/np.pi + Radius_deg=self.Radius + Dmax,NMax=self.options.UseGaiaDB.split(",") + Dmax,NMax=float(Dmax),int(NMax) + Parallax_min=1./(Dmax*1e-3) + query=f"""SELECT TOP 10000 gaia_source.designation,gaia_source.source_id,gaia_source.ref_epoch,gaia_source.ra,gaia_source.dec,gaia_source.parallax,gaia_source.pmra,gaia_source.pmdec + FROM gaiadr3.gaia_source + WHERE + CONTAINS( + POINT('ICRS',gaiadr3.gaia_source.ra,gaiadr3.gaia_source.dec), + CIRCLE('ICRS',{rac_deg},{decc_deg},{Radius_deg}) + )=1 AND (gaiadr3.gaia_source.parallax>={Parallax_min})""" + log.print(f"Sending du Gaia server") + print(f"{query}") + job = Gaia.launch_job(query) + result = job.get_results() + + log.print(f"Query has returned {len(result)}") + + + + l=[] + dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), + ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), + ('Type','S200')] + for r in result: + l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],b"Gaia DR3")) + self.PosArray=np.asarray(l,dtype=dtype) + + # import pylab + # pylab.clf() + # pylab.subplot(1,2,1) + # pylab.scatter(self.PosArray["ra"],self.PosArray["dec"]) + # pylab.scatter([rac_deg],[decc_deg]) + # pylab.subplot(1,2,2) + # D=1./(self.PosArray["parallax"]*1e-3) + # pylab.hist(D,bins=100) + # pylab.draw() + # pylab.show() + + self.DoProperMotionCorr=True + if self.PosArray.size>NMax: + ind=np.int64(np.random.rand(NMax)*self.PosArray.size) + self.PosArray=self.PosArray[ind] + + + elif self.options.FitsCatalog: + print("Using the fits catalog: %s"%self.options.FitsCatalog, file=log) + F=fits.open(self.options.FitsCatalog) + d=F[1].data + l=[] + dtype=[('Name','S200'),("ra",np.float64),("dec",np.float64), + ("pmra",np.float64),("pmdec",np.float64),("ref_epoch",np.float64),("parallax",np.float64), + ('Type','S200')] + for r in d: + l.append((r['DESIGNATION'],r['ra'],r['dec'],r['pmra'],r['pmdec'],r['ref_epoch'],r['parallax'],r['Type'])) + + if FileCoords is not None and FileCoords!="": + print('Adding data from file '+FileCoords, file=log) + dtype0=[('Name','S200'),("ra",np.float64),("dec",np.float64),('Type','S200')] + additional=np.genfromtxt(FileCoords,dtype=dtype0,delimiter=",")[()] + if len(additional.shape)==0: additional=additional.reshape((1,)) + if not additional.shape: + # deal with a one-line input file + additional=np.array([additional],dtype=dtype0) + additional1=np.zeros((additional.shape[0],),dtype=dtype) + additional1["Name"][:]=additional["Name"][:] + additional1["ra"][:]=additional["ra"][:] + additional1["dec"][:]=additional["dec"][:] + additional1["Type"][:]=additional["Type"][:] + + for r in additional1: + l.append(tuple(r)) + + + self.PosArray=np.asarray(l,dtype=dtype) + self.DoProperMotionCorr=True + else: + + #FileCoords="Transient_LOTTS.csv" + if FileCoords is None: + if not os.path.isfile(FileCoords): + ssExec="wget -q --user=anonymous ftp://ftp.strw.leidenuniv.nl/pub/tasse/%s -O %s"%(FileCoords,FileCoords) + print("Downloading %s"%FileCoords, file=log) + print(" Executing: %s"%ssExec, file=log) + os.system(ssExec) + log.print("Reading cvs file: %s"%FileCoords) + #self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",")[()] + self.PosArray=np.genfromtxt(FileCoords,dtype=dtype,delimiter=",") + if len(self.PosArray.shape)==0: self.PosArray=self.PosArray.reshape((1,)) + + self.PosArray=self.PosArray.view(np.recarray) + self.PosArray.ra*=np.pi/180. + self.PosArray.dec*=np.pi/180. + + Radius=self.Radius + self.NOrig=self.PosArray.Name.shape[0] + Dist=AngDist(self.ra0,self.PosArray.ra,self.dec0,self.PosArray.dec) + ind=np.where(Dist<(Radius*np.pi/180))[0] + self.PosArray=self.PosArray[ind] + + print("Created an array with %i records" % self.PosArray.size, file=log) + + if SubSet is not None: + random.seed(42) + i,n=SubSet + ind=random.sample(range(0, self.PosArray.size), self.PosArray.size) + ii=np.int64(np.linspace(0,self.PosArray.size+1,n+1)) + L=np.sort(ind[ii[i]:ii[i+1]]) + self.PosArray=self.PosArray[L] + print("Selected %i objects for subset %i/%i" % (self.PosArray.size,i+1,n), file=log) + + return self.PosArray diff --git a/ClassSaveResults.py b/ClassSaveResults.py index a6e4670..13c5211 100644 --- a/ClassSaveResults.py +++ b/ClassSaveResults.py @@ -20,7 +20,7 @@ import matplotlib.pyplot as pylab from DDFacet.Other.progressbar import ProgressBar from pyrap.images import image -from .dynspecms_version import version +from dynspecms_version import version import DDFacet.Other.MyPickle def GiveMAD(X): @@ -34,8 +34,8 @@ def __init__(self, DynSpecMS,DIRNAME=None): self.DIRNAME="DynSpecs_%s"%self.DynSpecMS.OutName else: self.DIRNAME="%s_DynSpecs_%s"%(self.DIRNAME,self.DynSpecMS.OutName) - - + + #image = self.DynSpecMS.Image #self.ImageData=np.squeeze(fits.getdata(image, ext=0)) @@ -44,7 +44,7 @@ def __init__(self, DynSpecMS,DIRNAME=None): self.im=self.imI=image(self.DynSpecMS.ImageI) self.ImageIData=self.imI.getdata()[0,0] - + self.ImageV=self.DynSpecMS.ImageV if self.ImageV and os.path.isfile(self.ImageV): self.imV=image(self.DynSpecMS.ImageV) @@ -59,7 +59,7 @@ def __init__(self, DynSpecMS,DIRNAME=None): ("IDTessel",np.int32),("IDFacet",np.int32), ("FluxI",np.float32),("FluxV",np.float32),("sigFluxI",np.float32),("sigFluxV",np.float32)]) self.CatFlux=self.CatFlux.view(np.recarray) - + os.system("rm -rf %s"%self.DIRNAME) os.system("mkdir -p %s/TARGET"%self.DIRNAME) os.system("mkdir -p %s/OFF"%self.DIRNAME) @@ -82,12 +82,18 @@ def WriteFits(self): if self.DynSpecMS.DoJonesCorr_kMS: self.CatFlux.IDFacet[:]=self.DynSpecMS.PosArray.iFacet[:] self.CatFlux.IDTessel[:]=self.DynSpecMS.PosArray.iTessel[:] - - + + pBAR = ProgressBar(Title="Saving fits") + NJobs=self.DynSpecMS.NDir #Selected + iDone=0 + pBAR.render(0, NJobs) + for iDir in range(self.DynSpecMS.NDir): self.WriteFitsThisDir(iDir) self.WriteFitsThisDir(iDir,Weight="Weight") self.WriteFitsThisDir(iDir,Weight="W2") + iDone+=1 + pBAR.render(iDone, NJobs) def SaveCatalog(self): FileName = "%s/%s.npy"%(self.DIRNAME,"Catalog") @@ -96,28 +102,28 @@ def SaveCatalog(self): if self.DynSpecMS.DFacet is not None: DDFacet.Other.MyPickle.Save(self.DynSpecMS.DFacet,"%s/%s.npy"%(self.DIRNAME,"DDF.DicoFacet")) self.radecToReg() - + def GiveSubDir(self,Type,Weight=False): SubDir="OFF" if Type!=b"Off": SubDir="TARGET" if Weight=="Weight" or Weight=="W2": SubDir+="_W" - + return SubDir def radecToReg(self): FName="%s/%s.reg"%(self.DIRNAME,self.DynSpecMS.OutName) ra,dec=self.DynSpecMS.PosArray.ra,self.DynSpecMS.PosArray.dec Type=self.DynSpecMS.PosArray.Type - + log.print(("Writting target reg file: %s"%FName)) f=open(FName,"w") - + f.write("""# Region file format: DS9 version 4.1\n""") f.write("""global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1\n""") f.write("""fk5\n""") - + sRA0=rad2hmsdms(self.DynSpecMS.ra0,Type="ra").replace(" ",":") sDEC0=rad2hmsdms(self.DynSpecMS.dec0,Type="dec").replace(" ",":") f.write("""circle(%s,%s,%f" # color=%s\n"""%(sRA0,sDEC0,self.DynSpecMS.Radius*3600,"green")) @@ -126,7 +132,7 @@ def radecToReg(self): ra0,dec0=ra[iTarget],dec[iTarget] sRA0=rad2hmsdms(ra0,Type="ra").replace(" ",":") sDEC0=rad2hmsdms(dec0,Type="dec").replace(" ",":") - + if Type[iTarget].decode('ASCII')=="Off": color="blue" else: @@ -138,25 +144,25 @@ def radecToReg(self): # f.write("""point(%s,%s # text={[F%i_S%i]} point=cross 5 color=%s\n"""%(sRA0,sDEC0,iFacet,iTessel,color)) f.close() - + def WriteFitsThisDir(self,iDir,Weight=False): """ Store the dynamic spectrum in a FITS file """ ra,dec=self.DynSpecMS.PosArray.ra[iDir],self.DynSpecMS.PosArray.dec[iDir] - + strRA=rad2hmsdms(ra,Type="ra").replace(" ",":") strDEC=rad2hmsdms(dec,Type="dec").replace(" ",":") - + fitsname = "%s/%s/%s_%s_%s.fits"%(self.DIRNAME,self.GiveSubDir(self.DynSpecMS.PosArray.Type[iDir],Weight=Weight),self.DynSpecMS.OutName, strRA, strDEC) self.CatFlux.FileName[iDir]=fitsname if Weight=="Weight": fitsname = "%s/%s/%s_%s_%s.W.fits"%(self.DIRNAME,self.GiveSubDir(self.DynSpecMS.PosArray.Type[iDir],Weight=Weight),self.DynSpecMS.OutName, strRA, strDEC) elif Weight=="W2": fitsname = "%s/%s/%s_%s_%s.W2.fits"%(self.DIRNAME,self.GiveSubDir(self.DynSpecMS.PosArray.Type[iDir],Weight=Weight),self.DynSpecMS.OutName, strRA, strDEC) - - print("#%i %s %s"%(iDir,self.DynSpecMS.PosArray.Type[iDir].decode("ascii"),fitsname),file=log) + + #print("#%i %s %s"%(iDir,self.DynSpecMS.PosArray.Type[iDir].decode("ascii"),fitsname),file=log) # Create the fits file - prihdr = fits.Header() + prihdr = fits.Header() prihdr.set('CTYPE1', 'Time', 'Time') prihdr.set('CRPIX1', 1., 'Reference') prihdr.set('CRVAL1', 0., 'Time at the reference pixel (sec since OBS-STAR)') @@ -192,15 +198,15 @@ def WriteFitsThisDir(self,iDir,Weight=False): if Weight: ThisType+="_%s"%Weight prihdr.set('SRC-TYPE', ThisType, 'Type of the source in the source list') - + prihdr.set('ORIGIN', 'DynSpecMS '+version(),'Created by') if "iFacet" in self.DynSpecMS.PosArray.dtype.fields.keys(): iFacet=self.DynSpecMS.PosArray.iFacet[iDir] iTessel=self.DynSpecMS.PosArray.iTessel[iDir] prihdr.set('FACET', iFacet, 'ID of the facet') prihdr.set('TESSEL', iTessel, 'ID of the Tessel') - - + + if Weight=="Weight": Gn = self.DynSpecMS.DicoGrids["GridWeight"][iDir,:, :, 0:1].real # dir, time, freq, pol elif Weight=="W2": @@ -221,7 +227,7 @@ def PlotSpec(self,Prefix=""): # Pdf file of target positions pdfname = "%s/%s_TARGET%s.pdf"%(self.DIRNAME,self.DynSpecMS.OutName,Prefix) print("Making pdf overview: %s"%pdfname, file=log) - pBAR = ProgressBar(Title="Making pages") + pBAR = ProgressBar(Title="Making pages") NPages=self.DynSpecMS.NDirSelected #Selected iDone=0 pBAR.render(0, NPages) @@ -256,7 +262,7 @@ def PlotSpec(self,Prefix=""): # # Pdf smoothed of target positions # pdfname = "%s/%s_TARGET_Smoothed%s.pdf"%(self.DIRNAME,self.DynSpecMS.OutName,Prefix) # print>>log,"Making pdf overview: %s"%pdfname - # pBAR = ProgressBar(Title="Making pages") + # pBAR = ProgressBar(Title="Making pages") # NPages=self.DynSpecMS.NDirSelected #Selected # iDone=0 # pBAR.render(0, NPages) @@ -270,7 +276,7 @@ def PlotSpec(self,Prefix=""): # iDone+=1 # pBAR.render(iDone, NPages) - + def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): label = ["I", "Q", "U", "V"] @@ -279,7 +285,7 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): pylab.rc('text', usetex=True) font = {'family':'serif', 'serif': ['Times']} pylab.rc('font', **font) - + # Figure properties bigfont = 8 @@ -320,7 +326,7 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): cmap = copy.copy(matplotlib.cm.get_cmap(cmap))#.copy() cmap.set_bad(color='black') - + ax1 = pylab.subplot(4, 1, ipol+1) spec = pylab.pcolormesh(times, freqs, Gn, cmap=cmap,#'bone_r', vmin=mean-3*sig, vmax=mean+10*sig, rasterized=True,shading='auto') @@ -343,11 +349,11 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): sig = GiveMAD(Gn) mean = np.median(Gn) #spec = pylab.pcolormesh(times, freqs, Gn, cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, rasterized=True) - spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6),rasterized=True) + spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6),rasterized=True) axspec.axis('tight') cbar = pylab.colorbar(fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=smallfont) - cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') + cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') pylab.text(times[-1]-0.02*(times[-1]-times[0]), freqs[-1]-0.1*(freqs[-1]-freqs[0]), 'I', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) pylab.text(times[0]+0.02*(times[-1]-times[0]), freqs[0]+0.1*(freqs[-1]-freqs[0]), r"$\sigma =$ %.3f Jy"%sig, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) pylab.xlabel("Time (min since %s)"%(t0.iso), fontsize=bigfont) @@ -360,13 +366,13 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): AG=np.abs(Gn) sig = GiveMAD(Gn) - mean = np.median(Gn) + mean = np.median(Gn) #spec = pylab.pcolormesh(times, freqs, Gn, cmap='bone_r', vmin=0, vmax=mean+10*sig, rasterized=True) - spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) + spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) axspec.axis('tight') cbar = pylab.colorbar(fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=smallfont) - cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') + cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') pylab.text(times[-1]-0.02*(times[-1]-times[0]), freqs[-1]-0.1*(freqs[-1]-freqs[0]), 'L', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) pylab.text(times[0]+0.02*(times[-1]-times[0]), freqs[0]+0.1*(freqs[-1]-freqs[0]), r"$\sigma =$ %.3f Jy"%sig, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) pylab.xlabel("Time (min since %s)"%(t0.iso), fontsize=bigfont) @@ -379,13 +385,13 @@ def PlotSpecSingleDir(self, iDir=0, BoxArcSec=300.): AG=np.abs(Gn) sig = GiveMAD(Gn) - mean = np.median(Gn) + mean = np.median(Gn) #spec = pylab.pcolormesh(times, freqs, Gn, cmap='bone_r', vmin=mean-3*sig, vmax=mean+10*sig, rasterized=True) - spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-5*sig, vmax=mean+5*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) + spec = pylab.imshow(Gn, interpolation="nearest", cmap='bone_r', vmin=mean-5*sig, vmax=mean+5*sig, extent=(times[0],times[-1],self.DynSpecMS.fMin*1.e-6,self.DynSpecMS.fMax*1.e-6), rasterized=True) axspec.axis('tight') cbar = pylab.colorbar(fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=smallfont) - cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') + cbar.set_label(r'Flux density (Jy)', fontsize=8, horizontalalignment='center') pylab.text(times[-1]-0.02*(times[-1]-times[0]), freqs[-1]-0.1*(freqs[-1]-freqs[0]), 'V', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) pylab.text(times[0]+0.02*(times[-1]-times[0]), freqs[0]+0.1*(freqs[-1]-freqs[0]), r"$\sigma =$ %.3f Jy"%sig, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) pylab.xlabel("Time (min since %s)"%(t0.iso), fontsize=bigfont) @@ -451,7 +457,7 @@ def giveBounded(x): x1=giveBounded(xc+box) y0=giveBounded(yc-box) y1=giveBounded(yc+box) - + DataBoxed=self.ImageIData[y0:y1,x0:x1] FluxI=self.ImageIData[yc,xc] sigFluxI=GiveMAD(DataBoxed) @@ -467,7 +473,7 @@ def giveBounded(x): im = pylab.imshow(DataBoxed, interpolation="nearest", cmap='bone_r', aspect="auto", vmin=vMin, vmax=vMax, origin='lower', rasterized=True) #pylab.text((ra_crop[1]-ra_crop[0])/16, (dec_crop[1]-dec_crop[0])/16, r"$\sigma =$ %.3f mJy"%rms, horizontalalignment='left', verticalalignment='center', fontsize=bigfont+2) cbar = pylab.colorbar()#(fraction=0.046*2., pad=0.01*4.) - + ax1.set_xlabel(r'RA (J2000)') raax = ax1.coords[0] raax.set_major_formatter('hh:mm:ss') @@ -483,13 +489,13 @@ def giveBounded(x): cbar.ax.tick_params(labelsize=smallfont) pylab.setp(ax1.get_xticklabels(), rotation='horizontal', fontsize=smallfont) pylab.setp(ax1.get_yticklabels(), rotation='horizontal', fontsize=smallfont) - + ra_cen, dec_cen = wcs.wcs_world2pix(np.degrees(self.DynSpecMS.PosArray.ra[iDir]), np.degrees(self.DynSpecMS.PosArray.dec[iDir]), 1) #pylab.plot(ra_cen, dec_cen, 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target #pylab.plot(newra_cen, newdec_cen, 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target pylab.plot(DataBoxed.shape[1]/2., DataBoxed.shape[0]/2., 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target - - pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'I', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) + + pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'I', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) # ---- Image V ---- @@ -524,7 +530,7 @@ def giveBounded(x): self.CatFlux.Type[iDir]=self.DynSpecMS.PosArray.Type[iDir] self.CatFlux.ra[iDir]=self.DynSpecMS.PosArray.ra[iDir] self.CatFlux.dec[iDir]=self.DynSpecMS.PosArray.dec[iDir] - + newra_cen, newdec_cen = wcs.wcs_pix2world( (x1+x0)/2., (y1+y0)/2., 1) wcs.wcs.crpix = [ DataBoxed.shape[1]/2., DataBoxed.shape[0]/2. ] # update the WCS object wcs.wcs.crval = [ newra_cen, newdec_cen ] @@ -549,7 +555,7 @@ def giveBounded(x): pylab.setp(ax1.get_yticklabels(), rotation='horizontal', fontsize=smallfont) ra_cen, dec_cen = wcs.wcs_world2pix(np.degrees(self.DynSpecMS.PosArray.ra[iDir]), np.degrees(self.DynSpecMS.PosArray.dec[iDir]), 1) pylab.plot(DataBoxed.shape[1]/2., DataBoxed.shape[0]/2., 'o', markerfacecolor='none', markeredgecolor='red', markersize=bigfont) # plot a circle at the target - pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'V', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) + pylab.text(DataBoxed.shape[0]*0.9, DataBoxed.shape[1]*0.9, 'V', horizontalalignment='center', verticalalignment='center', fontsize=bigfont+2) #pylab.subplots_adjust(wspace=0.15, hspace=0.30) diff --git a/dynamic_spectra_analysis.ipynb b/dynamic_spectra_analysis.ipynb new file mode 100644 index 0000000..b4e019f --- /dev/null +++ b/dynamic_spectra_analysis.ipynb @@ -0,0 +1,316 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import fits\n", + "from astropy.wcs import WCS\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from astropy import units as u\n", + "\n", + "from matplotlib.backends.backend_pdf import PdfPages\n", + "from matplotlib import patches\n", + "from matplotlib.patches import Ellipse\n", + "from decimal import Decimal" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.rcParams['font.size'] = 7" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "target = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET/1556208112_14:29:42.950_-62:40:46.200.fits'\n", + "target_w = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET_W/1556208112_14:29:42.950_-62:40:46.200.W.fits'\n", + "target_w2 = '~/Downloads/MSName_RandChunk0_DynSpecs_1556208112/TARGET_W/1556208112_14:29:42.950_-62:40:46.200.W2.fits'\n", + "\n", + "target_hdu = fits.open(target)[0]\n", + "target_w_hdu = fits.open(target_w)[0]\n", + "target_w2_hdu = fits.open(target_w2)[0]\n", + "\n", + "data = target_hdu.data\n", + "weights = target_w_hdu.data\n", + "weights2 = target_w2_hdu.data\n", + "\n", + "header = target_hdu.header" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Loop over the items and print each one\n", + "for keyword in header:\n", + "# print(f'{keyword} = {header[keyword]}')\n", + " frq_min = header['FRQ-MIN'] * 1e-6\n", + " frq_max = header['FRQ-MAX'] * 1e-6\n", + " tmin = 0\n", + " tmax = data.shape[2] * header['CDELT1']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "weights_data = np.where(weights==0,1, weights)\n", + "weights_data2 = np.where(weights2==0,1, weights2)\n", + "\n", + "sigma_inv = np.sqrt(np.power(weights_data, 2) / weights_data2)\n", + "\n", + "for i in range(data.shape[0]):\n", + " #data[i, :, :] = data[i, :, :] * norm_w * np.cbrt(data[i, :, :])\n", + " data[i, :, :] = data[i, :, :] * sigma_inv\n", + "\n", + "print('--- Data has shape ---')\n", + "print('(stokes, freq, time)')\n", + "print(data.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "stokes_idx = 0\n", + "d = data[stokes_idx, :, :]\n", + "extent=[tmin, tmax, frq_min, frq_max]\n", + "plt.imshow(d, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(d), vmax=np.max(d), extent=extent)\n", + "plt.title('Dynamic Spectrum of total intensity (Stokes I)')\n", + "plt.xlabel('time (s)')\n", + "plt.ylabel(r'$\\nu$ (MHz)')\n", + "cbar_t = plt.colorbar(pad=.07)\n", + "cbar_t.set_label(r'$mJy/B$', size=7)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Stokes I, Q, U, V Plotted in order below\n", + "```python\n", + "I = data[0, :, :]\n", + "Q = data[1, :, :]\n", + "U = data[2, :, :]\n", + "V = data[3, :, :]\n", + "```\n", + "\n", + "For now, we slice the time (t0, t1):\n", + "```sh\n", + "6325, 6345\n", + "12953, 12993\n", + "19803, 19843\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_spectral_features(s, t0, t1):\n", + " data_sl = data[s, :, t0:t1]\n", + " t = []\n", + " if not np.max(data_sl) == 0:\n", + " t.append(t0)\n", + " print(\"Flaring from {} to {} seconds\".format(t0, t1))\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_stokes(stokes_symbol, stokes_idx, t0, t1):\n", + " \"\"\"Plot dynamic spectra for given data array and Stokes param\n", + " \n", + " Parameters\n", + " ----------\n", + " stokes_symbol: str\n", + " Stokes parameter of interest mapped ('I', 'Q', 'U', 'V')\n", + " stokes_idx: int\n", + " Stokes parameter of interest mapped ('I', 'Q', 'U', 'V')->(0, 1, 2, 3)\n", + " t0: int\n", + " `start_time` to zoom-in on\n", + " t1: int\n", + " `end_time to` zoom-in on\n", + " \"\"\"\n", + " data_sl = data[stokes_idx, :, t0:t1]\n", + " print(\"--- Stats for Stokes {} ---\".format(stokes_symbol))\n", + " print(\"----Slicing throug time axis from {} to {} seconds after start of observation\".format(t0, t1))\n", + " print(\"Min: \", np.min(data_sl))\n", + " print(\"Max: \", np.max(data_sl))\n", + " print(\"Noise: \", np.std(data_sl))\n", + " print(\"Mean: \", np.mean(data_sl))\n", + " #print(np.argmax(d[stokes_idx, :, :]))\n", + "\n", + " # Plot\n", + " extent = [t0, t1, frq_min, frq_max]\n", + " plt.title(stokes_symbol)\n", + " plt.imshow(data_sl, origin='lower', cmap='viridis', aspect='auto', vmin=np.min(data_sl), vmax=np.max(data_sl), extent=extent)\n", + " plt.xlabel('time (s)')\n", + " plt.ylabel(r'$\\nu$ (MHz)')\n", + " cbar = plt.colorbar(pad=.07)\n", + " cbar.set_label(r'$mJy/B$', size=7)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(0, data.shape[2], 20):\n", + " t0, t1 = i, i+20\n", + " get_spectral_features(0, t0, t1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Flaring events every 35 minutes that repeat every 20 seconds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(19820-17440)/60" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_stokes('Stokes I', 0, 0, 40)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes I', 0, 12953, 12993)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes I', 0, 19803, 19843)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes I', 0, 19803, 19843)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes Q', 1, 19803, 19843)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes U', 2, 19803, 19843)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot_stokes(data, 'Stokes V', 3, 19803, 19843)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### During Fast Imaging meeting\n", + "Oleg advised I try this\n", + "``` sh\n", + "pip install owlcat\n", + "plot-evelation-tracks.py \n", + "```\n", + "This will help us determine if the features we see were during the duty cycle of the target source" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![](lst-elev.png)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lst-elev.png b/lst-elev.png new file mode 100644 index 0000000..4012fd4 Binary files /dev/null and b/lst-elev.png differ diff --git a/ms2dynspec.py b/ms2dynspec.py index 106468a..030bbdc 100755 --- a/ms2dynspec.py +++ b/ms2dynspec.py @@ -7,19 +7,21 @@ from past.builtins import cmp __author__ = "Cyril Tasse, and Alan Loh" __credits__ = ["Cyril Tasse", "Alan Loh"] -from DynSpecMS.dynspecms_version import version +from dynspecms_version import version from DDFacet.Other import logger log=logger.getLogger("ms2dynspec") from DDFacet.Other import ModColor __version__ = version() +import numpy as np SaveFile = "last_dynspec.obj" +from ClassDynSpecMS import ClassGiveCatalog """ ========================================================================= DESCRIPTION Blablabla - Modif - + Modif + Example: python ms2dynspec.py --msfile /usr/data/ --data CORRECTED --model PREDICT --sols killms.npz --srclist SRCPOS.txt --rad 10 @@ -52,10 +54,10 @@ import glob, os import pylab from DDFacet.Other import MyPickle -from DynSpecMS import logo +import logo logo.PrintLogo(__version__) -from DynSpecMS.ClassDynSpecMS import ClassDynSpecMS -from DynSpecMS import ClassSaveResults +from ClassDynSpecMS import ClassDynSpecMS +from ClassSaveResults import ClassSaveResults from DDFacet.Data.ClassMS import expandMSList from DDFacet.Other import ModColor from DDFacet.Other import progressbar @@ -86,6 +88,8 @@ def angSep(ra1, dec1, ra2, dec2): def main(args=None, messages=[]): if args is None: args = MyPickle.Load(SaveFile) + if args.UseRandomSeed!=0: + np.random.seed(args.UseRandomSeed) MSList=None if args.ms: @@ -103,52 +107,88 @@ def main(args=None, messages=[]): else: DT[T].append(MSName) t.close() + if len(DT)>1: log.print(ModColor.Str("FOUND %i time periods"%len(DT))) else: DT={0:MSList} - for k in DT.keys(): - MSList=DT[k] - D = ClassDynSpecMS(ListMSName=MSList, - ColName=args.data, ModelName=args.model, - SolsName=args.sols, - TChunkHours=args.TChunkHours, - ColWeights=args.WeightCol, - UVRange=args.uv, - FileCoords=args.srclist, - Radius=args.rad, - NOff=args.noff, - DicoFacet=args.DicoFacet, - ImageI=args.imageI, - ImageV=args.imageV, - SolsDir=args.SolsDir,NCPU=args.NCPU, - BaseDirSpecs=args.BaseDirSpecs, - BeamModel=args.BeamModel, - BeamNBand=args.BeamNBand, - SourceCatOff_FluxMean=args.SourceCatOff_FluxMean, - SourceCatOff_dFluxMean=args.SourceCatOff_dFluxMean, - SourceCatOff=args.SourceCatOff, - options=args) - - if D.NDirSelected==0: - return - - if D.Mode=="Spec": D.StackAll() - - SaveMachine=ClassSaveResults.ClassSaveResults(D,DIRNAME=args.OutDirName) - if D.Mode=="Spec": - SaveMachine.WriteFits() - if args.SavePDF: - SaveMachine.PlotSpec() - SaveMachine.SaveCatalog() - if args.DoTar: SaveMachine.tarDirectory() - else: - SaveMachine.SaveCatalog() - SaveMachine.PlotSpec(Prefix="_replot") - D.killWorkers() + L_radec=[] + for MSName in MSList: + tField = table("%s::FIELD"%MSName, ack=False) + ra0, dec0 = tField.getcol("PHASE_DIR").ravel() + if ra0<0.: ra0+=2.*np.pi + L_radec.append((ra0,dec0)) + tField.close() + L_radec=list(set(L_radec)) + if len(L_radec)>1: stop + ra0,dec0=L_radec[0] + + NChunk=1 + if args.NMaxTargets!=0: + CGC=ClassGiveCatalog.ClassGiveCatalog(args, + ra0,dec0, + args.rad, + args.srclist) + PosArray=CGC.giveCat() + NChunk=PosArray.size//args.NMaxTargets+1 + log.print(ModColor.Str("Will process the catalog in %i chunks"%NChunk,col="green")) + + + SubSet=None + for iChunk in range(NChunk): + log.print("====================================================================") + if NChunk>1: + SubSet=(iChunk,NChunk) + for ik,k in enumerate(sorted(list(DT.keys()))): + MSList=DT[k] + DIRNAME="%s"%(args.OutDirName) + if len(DT)>1: + DIRNAME+="_T%i"%ik + if NChunk>1: + DIRNAME+="_RandChunk%i"%iChunk + + + D = ClassDynSpecMS(ListMSName=MSList, + ColName=args.data, ModelName=args.model, + SolsName=args.sols, + TChunkHours=args.TChunkHours, + ColWeights=args.WeightCol, + UVRange=args.uv, + FileCoords=args.srclist, + Radius=args.rad, + NOff=args.noff, + DicoFacet=args.DicoFacet, + ImageI=args.imageI, + ImageV=args.imageV, + SolsDir=args.SolsDir,NCPU=args.NCPU, + BaseDirSpecs=args.BaseDirSpecs, + BeamModel=args.BeamModel, + BeamNBand=args.BeamNBand, + SourceCatOff_FluxMean=args.SourceCatOff_FluxMean, + SourceCatOff_dFluxMean=args.SourceCatOff_dFluxMean, + SourceCatOff=args.SourceCatOff, + options=args, + SubSet=SubSet) + + if D.NDirSelected==0: + return + + if D.Mode=="Spec": D.StackAll() + + SaveMachine=ClassSaveResults(D,DIRNAME=DIRNAME) + if D.Mode=="Spec": + SaveMachine.WriteFits() + if args.SavePDF: + SaveMachine.PlotSpec() + SaveMachine.SaveCatalog() + if args.DoTar: SaveMachine.tarDirectory() + else: + SaveMachine.SaveCatalog() + SaveMachine.PlotSpec(Prefix="_replot") + D.killWorkers() Multiprocessing.cleanupShm() - + # ========================================================================= # ========================================================================= @@ -157,7 +197,7 @@ def main(args=None, messages=[]): parser.add_argument("--ms", type=str, help="Name of MS file / directory", required=False) parser.add_argument("--data", type=str, default="CORRECTED", help="Name of DATA column", required=False) parser.add_argument("--TChunkHours", type=float, default=0., help="Chunk size in hours", required=False) - + parser.add_argument("--WeightCol", type=str, default=None, help="Name of weights column to be taken into account", required=False) parser.add_argument("--model", type=str, help="Name of MODEL column",default="")#, required=True) parser.add_argument("--sols", type=str, help="Jones solutions",default="") @@ -178,7 +218,8 @@ def main(args=None, messages=[]): parser.add_argument("--UseLoTSSDB", type=int, default=0, help="Use LoTSS DB for target list", required=False) parser.add_argument("--UseGaiaDB", type=str, default=None, help="Use Gaia DB for target list", required=False) parser.add_argument("--DoTar", type=int, default=1, help="Tar final products", required=False) - + parser.add_argument("--UseRandomSeed", type=int, default=0, help="Use random seed", required=False) + parser.add_argument("--NCPU", type=int, default=0, help="NCPU", required=False) parser.add_argument("--BeamModel", type=str, default=None, help="Beam Model to be used", required=False) parser.add_argument("--BeamNBand", type=int, default=1, help="Number of channels in the Beam Jones matrix", required=False) @@ -187,6 +228,7 @@ def main(args=None, messages=[]): parser.add_argument("--SourceCatOff", type=str, default="", help="Read the code", required=False) parser.add_argument("--SourceCatOff_FluxMean", type=float, default=0, help="Read the code", required=False) parser.add_argument("--SourceCatOff_dFluxMean", type=float, default=0, help="Read the code", required=False) + parser.add_argument("--NMaxTargets", type=int, default=0, help="Read the code", required=False) args = parser.parse_args()