该博客主要涉及利用Python处理CDMS数据库和Molfit文件，进行分子线的搜索和证认。通过编写查询字符串，筛选频率和能量范围内的跃迁，并对已证认的分子进行操作。使用XCLASS包进行拟合，提取超过3σ的拟合结果。此外，还介绍了如何将Molfit文件拆分为单个文件并分析识别的过渡。代码中包含了数据库读取、数据处理和文件操作的详细步骤。
有两个功能没写好。
XCLASS的分子线证认已经写得很好，但是后续还要自己做一些比较繁杂的工作，这里放了部分的操作步骤。
根据CDMS的database文件和已经证认好分子的molfit文件，这里可以搜寻一下，得出拟合结果大于3sigma的跃迁。
下面的代码是写好的package。
#!/usr/bin/python
# -*- coding: utf-8 -*-
##********************************************************************************************************************************************************
##********************************************************************* load packages ********************************************************************
from __future__ import print_function ## import print_function package
import numpy ## import numpy package
import os ## import os package
import sys ## import sys package
import string ## import string package
import random ## import random package
import datetime ## import datetime package
import sqlite3 ## import sqlite3 package
from pylab import * ## import pylab package
import task_MAGIX ## import MAGIX function from XCLASS package
import task_myXCLASS ## import package myXCLASS
import task_myXCLASSFit ## import package myXCLASSFit
from collections import Counter
##--------------------------------------------------------------------------------------------------------------------------------------------------------
##--------------------------------------------------------------------------------------------------------------------------------------------------------
##
##
def QueryString(FreqMin,FreqMax,ElowMin, ElowMax,ListOfSelectedMolecules):
## define frequency interval form mouse position
FreqIntMin = FreqMin
FreqIntMax = FreqMax
## initialize some parameter
NameOfRadTransTable = "transitions"
ColumnNameForNameTransitions = "T_Name"
ColumnNameForFreqTransitions = "T_Frequency"
ColumnNameForIntTransitions = "T_Intensity"
ColumnNameForEinsteinATransitions = "T_EinsteinA"
ColumnNameForFreqErrTransitions = "T_Uncertainty"
ColumnNameForELowTransitions = "T_EnergyLower"
ColumnNameForgUpTransitions = "T_UpperStateDegeneracy"
ColumnNameForQNUpLabelTransitions = "T_UpperStateQuantumNumbers"
ColumnNameForQNLowLabelTransitions = "T_LowerStateQuantumNumbers"
ColumnNameForURLTransitions = "T_URL"
## define query string
##
## get names of columns with
## import sqlite3
## connection = sqlite3.connect('cdms_sqlite.db')
## cursor = connection.execute('select * from partitionfunctions')
## names = list(map(lambda x: x[0], cursor.description))
##
query_string = "SELECT "
query_string += ColumnNameForNameTransitions + ", " ## position 1: molecule name
query_string += ColumnNameForFreqTransitions + ", " ## position 2: frequency
query_string += ColumnNameForIntTransitions + ", " ## position 3: intensity
query_string += ColumnNameForEinsteinATransitions + ", " ## position 4: Einstein A
query_string += ColumnNameForFreqErrTransitions + ", " ## position 5: Error frequency
query_string += ColumnNameForELowTransitions + ", " ## position 6: E_low
query_string += ColumnNameForgUpTransitions + ", " ## position 7: upper state degeneracy
query_string += ColumnNameForQNUpLabelTransitions + ", " ## position 8: quantum number label for upper state
query_string += ColumnNameForQNLowLabelTransitions ## position 9: quantum number label for lower state
# query_string += ColumnNameForURLTransitions ## position 10: URL
query_string += " FROM " + NameOfRadTransTable
query_string += " WHERE (" + ColumnNameForFreqTransitions + " >= " + str(FreqIntMin)
query_string += " and " + ColumnNameForFreqTransitions + " <= " + str(FreqIntMax)
if (ElowMin < ElowMax):
query_string += " and " + ColumnNameForELowTransitions + " >= " + str(ElowMin)
query_string += " and " + ColumnNameForELowTransitions + " <= " + str(ElowMax)
if (len(ListOfSelectedMolecules) > 0):
if (ListOfSelectedMolecules[0] != ["None"]):
query_string += " and ("
counter = 0
for molecule in ListOfSelectedMolecules:
counter += 1
if (counter > 1):
query_string += " or "
query_string += ColumnNameForNameTransitions + " = " + chr(34) + molecule.strip() + chr(34)
query_string += ")"
query_string += ")"
query_string += " ORDER by " + ColumnNameForFreqTransitions
# Debug:
# print " "
# print ("\n\nquery_string = ", query_string)
# print ("len(query_string) = ", len(query_string))
return (query_string)
def GetTrans(rows):
names = []
freqs = []
for row in rows: ## loop over all entries of the database
element_counter = 0
for elements in row:
element_counter += 1
if (element_counter == 1): ## get name of molecule
name = elements.strip()
name = '%25s' % name[:25] ## name of molecule has max. 25 characters
names.append(name)
elif (element_counter == 2): ## get frequency (in MHz)
freqVal = '%17s' % elements
if (freqVal.strip() != "None"):
freqVal = float(freqVal)
freq = '%17.5f' % freqVal
freqs.append(freq)
else:
freq = freqVal
freqs.append(freq)
return (names,freqs)
def MKMolFileOfMolecule(MolfitsFile,path,spw,molecule):
f = open(MolfitsFile)
lines = f.readlines()
f.close()
for i in range(len(lines)):
CurrentLine = lines[i].strip()
w = CurrentLine.find("%")
if (w == 0):
CurrentLine = ""
elif (w > 0):
CurrentLine = CurrentLine[:w]
if (CurrentLine != ""):
SplitLines = CurrentLine.split()
helpstring = SplitLines[-1].strip()
else:
helpstring = ''
if (helpstring.isdigit()):
if int(helpstring)>3:
print ("helpstring is %s" %(helpstring))
print ("If helpstring larger than 3, lease modify the code!")
if (SplitLines[0]==molecule):
if '#' in molecule:
molfitfilename = string.replace(molecule,';','_')
else:
molfitfilename = string.replace(molecule[:-1],';','_')
print(molfitfilename)
if int(helpstring)==1:
numpy.savetxt(path+"/molfitfilein%s/%s.molfit" %(spw,molfitfilename),[CurrentLine,lines[i+1].strip()],fmt='%s')
elif int(helpstring)==2:
numpy.savetxt(path+"/molfitfilein%s/%s.molfit" %(spw,molfitfilename),[CurrentLine,lines[i+1].strip(),lines[i+2].strip()],fmt='%s')
elif int(helpstring)==3:
numpy.savetxt(path+"/molfitfilein%s/%s.molfit" %(spw,molfitfilename),[CurrentLine,lines[i+1].strip(),lines[i+2].strip(),lines[i+3].strip()],fmt='%s')
def GetIdentifiedTransitions(names,freqs,FreqStep,modeldata,expdata,molecule,sigma):
name = []
freq = []
value = []
tag_molecule = []
newmolname = tagmoleculenames(molecule)
for i in range(len(names)):
if (names[i].encode('raw_unicode_escape').strip()==molecule):
model_position = numpy.argwhere(numpy.array(modeldata)[:,0]-FreqStep/2>float(freqs[i].strip()))[0]
exp_position = numpy.argwhere(numpy.array(expdata)[:,0]-FreqStep/2>float(freqs[i].strip()))[0]
model_x,model_y = modeldata[:,0],modeldata[:,1]
exp_x,exp_y = expdata[:,0],expdata[:,1]
if model_y[model_position]>3*sigma:
tag_molecule.append(newmolname)
__name = names[i].encode('raw_unicode_escape').strip()
w = __name.find('#')
if (w<0):
__name = __name[:-1]
name.append(__name)
freq.append(float(freqs[i].strip()))
value.append(float(exp_y[exp_position]))
return (tag_molecule,name,freq,value)
def tagmoleculenames(molname):
new_molname = ''
value = molname.find('#')
if (value>0):
molname = molname[:value]
molname = molname[:-1]
if ('-' in molname):
start = molname.find('-')
if ('-' in molname[start+1:]):
end = start+1+molname[start+1:].find('-')
superscr = molname[start+1:end]
molname_extracted = string.replace(molname,molname[start:end+1],'')
subscr = []
count = 0
a = 0
for i in molname_extracted[:molname_extracted.find(';')]:
count+=1
if (i.isdigit()):
subscr.append(i)
new_molname = molname_extracted[:molname_extracted.find(';')]
if (count<start):
a += 1
b = dict(Counter(subscr))
c = b.keys()
for num in range(len(c)):
new_molname = string.replace(new_molname,c[num],'$_'+c[num]+'$')
if num ==0:
new_molname = new_molname+molname_extracted[molname_extracted.find(';'):]
if (subscr==[]):
new_molname = molname_extracted
print(superscr)
new_molname = string.replace(new_molname,new_molname[start+a*3-1:],'$^{'+superscr+'}$'+new_molname[start+a*3-1:])
else:
subscr = []
for i in molname[:molname.find(';')]:
if (i.isdigit()):
subscr.append(i)
new_molname = molname[:molname.find(';')]
b = dict(Counter(subscr))
a = b.keys()
for num in range(len(a)):
new_molname = string.replace(new_molname,a[num],'$_'+a[num]+'$')
if num ==0:
new_molname = new_molname+molname[molname.find(';'):]
if (len(subscr)==0):
new_molname = molname
return (new_molname)
def ListSelectMolecule(SelectMolecule):
##====================================================================================================================================================
## test, if SelectMolecule is a list, otherwise try to use SelectMolecule as filename
print (" ")
print (" ")
print ("Analyze selected molecules ..")
sys.stdout.flush()
ListOfSelectedMolecules = []
if (len(SelectMolecule) > 0 and SelectMolecule != ["None"]):
for entry in SelectMolecule: ## loop over all entries of SelectMolecule list
entryStrip = entry.strip() ## remove leading and trailing blanks
## check, if current entry is a file
if (os.path.isfile(entryStrip)):
##----------------------------------------------------------------------------------------------------------------------------------------
## if file is a molfit file, get names of molecules from molfit file
if (entryStrip.endswith("molfit")):
# Debug:
# print "Selection string is a molfit file!"
## read in contents of molfit file
f = open(entryStrip)
MolfitContents = f.readlines() ## read in contents of molfit file
f.close()
## get names of molecules
for lines in MolfitContents: ## loop over all lines of molfit file
CurrentLine = lines.strip() ## remove leading and trailing blanks
# Debug:
# print "CurrentLine = sys.stdout.flush() Python的作用", CurrentLine
## check for comments
w = CurrentLine.find("%") ## are there comments in the current line ?
if (w == 0): ## ignore lines which contains only comments
CurrentLine = "" ## if line is only a comment clear line
elif (w > 0): ## is there a comment in the current line ?
CurrentLine = CurrentLine[:w] ## remove comments
# Debug:
# print "CurrentLine = ", CurrentLine
## analyze, if the current line contains the name of a molecule
if (CurrentLine != ""): ## ignore empty lines
SplitLines = CurrentLine.split() ## split current line into columns
# Debug:
# print "SplitLines = ", SplitLines
## contains the last column an integer number?
helpstring = SplitLines[-1].strip()
is_int = "false"
if (helpstring.isdigit()): ## continue, if string is an integer number
is_int = "true"
## If yes, then the current line contains the name of a molecule
if (is_int == "true"): ## if last entry in a line is an integer number, then
MolName = "" ## the current line contains the name of a molecule
for col in SplitLines[:-1]: ## get whole line without integer number at the end
MolName += col ## construct entry for SelectMolecule
# Debug:
# print "MolName = ", MolName
ListOfSelectedMolecules.append(MolName) ## construct final array containing all molecule names
##----------------------------------------------------------------------------------------------------------------------------------------
## read names of molecules from ASCII file
else:
# Debug:
# print "Selection string is a ASCII file!"
## read in whole ASCII file
f = open(entryStrip)
AllLines = f.readlines()
f.close()
## append lines in ASCII file to selection list
for line in AllLines:
ListOfSelectedMolecules.append(line)
##--------------------------------------------------------------------------------------------------------------------------------------------
## continue here, if current entry is not a filename
else:
ListOfSelectedMolecules.append(entryStrip)
# Debug:
# print ("ListOfSelectedMolecules = ", ListOfSelectedMolecules)
return (ListOfSelectedMolecules)
## read data from database
def ReadDatabase(query_string):
##====================================================================================================================================================
## get path of default database file
dbFilename = task_myXCLASS.GetDefaultDBFile()
##====================================================================================================================================================
## connect to the sqlite3 database
## connect to database
try:
conn = sqlite3.connect(dbFilename)
print (" ")
print ("Connection to sqlite3 database %s established." % dbFilename)
print (" ")
print (" ")
except sqlite3.Error, e:
print (" ")
print ("Can not connect to sqlite3 databse %s." % dbFilename)
print ("Error: %d: %s" % (e.args[0], e.args[1]))
sys.exit(1)
cursor = conn.cursor()
cursor.execute(query_string)
rows = cursor.fetchall()
return (rows)
这里是在casa里运行的调用上面package的代码。
#data
from task_SelectTransitions import *
import os
#from pandas import DataFrame
import numpy as np
##remove the old runing data
import shutil
path = '/home/zhang/software/XCLASS-Interface/run/myXCLASS/'
if os.path.exists(path)==False:
os.makedirs(path)
files = os.listdir(path)
if len(files)>0:
shutil.rmtree(path)
os.makedirs(path)
spw = "spw35"
sigma = 0.3
if spw=='spw35':
FreqMin = 9.753638179980000132e+04
FreqMax = 9.941000529416945938e+04
else:
FreqMin = 9.947395539869999629e+04
FreqMax = 1.013475792277295404e+05
#data off
#FreqMin = 219345
#FreqMax = 221323
FreqStep = 0.4883
#######################
TelescopeSize = 1.000
######################
Inter_Flag = True
t_back_flag = False
#tBack = 0
#tslope = 0.0000000000E+00
nH_flag = False
#N_H = 3.0000000000E+24
#beta_dust = 2.0
#kappa_1300 = 0.02
iso_flag = False
#IsoTableFileName = " "
RestFreq = 0.0
vLSR = 0.0
MolfitsFile = "/home/zhang/data/ATOMS/spec/I17233-3606/3/%s.molfit" %(spw)
FileName = "/home/zhang/data/ATOMS/spec/I17233-3606/3/I17233-3606_%s_17h26m42.4915s_-36d9m18.1152s.txt" %(spw)
NumHeaderLines = 0
expdata = LoadASCIIFile()
MinTransEnergy = 0.0
#data
#xLowerLimit = 9.753638179980000132e+04
#xUpperLimit = 9.941000529416945938e+04
#data off
#xLowerLimit = 219345
#xUpperLimit = 221323
SelectMolecule = ["/home/zhang/data/ATOMS/spec/I17233-3606/3/%s.molfit" %(spw)]
FrequencyWidth = 5.0
yLowerLimit = -8
yUpperLimit = 109
PlotTitle = "I17233_%s" %(spw)
LegendFlag = True
SaveFigureFile = " "
#OutputDevice = " "
ElowMin = 0.0
ElowMax = 1200.0
ListOfSelectedMolecules = ListSelectMolecule(SelectMolecule)
query_string = QueryString(FreqMin,FreqMax,ElowMin, ElowMax,ListOfSelectedMolecules)
rows = ReadDatabase(query_string)
names,freqs = GetTrans(rows)
path = os.getcwd()
direxists = os.path.exists(path+'/molfitfilein%s/' %(spw))
if direxists==False:
os.makedirs('molfitfilein%s' %(spw))
if (os.path.exists('tag_transitions_%s.txt' %(spw))==False):
tag_molecules = []
note_molecules = []
tag_x_positions = []
tag_y_positions = []
for molecule in ListOfSelectedMolecules:
MKMolFileOfMolecule(MolfitsFile,path,spw,molecule)
if '#' in molecule:
molfitfilename = string.replace(molecule,';','_')
else:
molfitfilename = string.replace(molecule[:-1],';','_')
MolfitsFileName = path+"/molfitfilein%s/%s.molfit" %(spw,molfitfilename)
modeldata,log,TransEnergies,IntOptical,jobDir = myXCLASS ()
tag_molecule,name,freq,value = GetIdentifiedTransitions(names,freqs,FreqStep,modeldata,expdata,molecule,sigma)
tag_molecules.append(tag_molecule)
note_molecules.append(name)
tag_x_positions.append(freq)
tag_y_positions.append(value)
tag_molecules.remove([])
note_molecules.remove([])
tag_x_positions.remove([])
tag_y_positions.remove([])
tag = []
for a in tag_molecules:
for b in a:
tag.append(b)
x = []
for a in tag_x_positions:
for b in a:
x.append(b)
y = []
for a in tag_y_positions:
for b in a:
y.append(b)
note = []
for a in note_molecules:
for b in a:
note.append(b)
########
# saved txt file
ind = np.lexsort((y,x))
data = [[tag[i],x[i],y[i],note[i]] for i in ind]
np.savetxt('tag_transitions_%s.txt' %(spw), data, delimiter=' ',fmt='%s')
else:
pass
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
