"""This script shows how to use the line identification tool""" # Scenario # -------- # # You have a spectral scan and you'd like to identify the lines in it. # For the sake of simplicity, we'll assume you cleaned this spectral scan # (no fringe, no spurs,...) # # In this guide, we'll use a small region of the spectrum of Orion South in band # 1a between 485 and 490 GHz # Notes # --------------- # # This document consists of 3 main sections: # - Basic usage # This is where you'll find the minimum you need to identify the lines of a # spectrum. # # - Advanced usage # For people who want to get the maximum out of the line identification # # - A guided tour # The line identification comes with handy functionalities for handling, # filtering, visualizing your lines. This section will teach you how to # do that. # # Minimal version: HIPE 13.0 build 3876 # ----------------------------------------------------------------------------- # Initialization # ----------------------------------------------------------------------------- import os from herschel.hifi.dp.tools.linelist import Linelist # Set the correct path to line_identification_guide/data: INPUT_DATA_DIR = "/Users/sfbeaulieu/Desktop/data-“lines” # ----------------------------------------------------------------------------- # Basic usage # ----------------------------------------------------------------------------- # In order to identify the lines in a spectrum you need at least 4 things: # - a spectrum: Spectrum1d # - a linelist file: the path to the CASSIS linelist file used as the known lines # - a vlsr: the VLSR of the source # - identifyLines: the task for the line identification # The first 3 things are the data: linelistFile = os.path.join(INPUT_DATA_DIR, "16293_SIS.txt") spectrum = fitsReader(os.path.join(INPUT_DATA_DIR, "Orion_S_baselined_485_490.fits")) vlsr = 6.7 # The line identification is performed with the task identifyLines. Pass the data # defined above to the task: lines = identifyLines(data=spectrum, linelist=linelistFile, vlsr=vlsr) # You get the identified, unidentified lines and the baselined spectrum identified = lines["identified"] unidentified = lines["unidentified"] baselinedSpectrum = lines["baseline"] # Note: Whilst the command # # lines = identifyLines(data=spectrum, linelist=linelistFile, vlsr=vlsr) # # may be enough to perform the line identification, it makes many assumptions # and sets default values to the optional parameters. Some users may want to # have more control on what's actually happening, that's what we'll describe in # the section "Advanced usage". # ----------------------------------------------------------------------------- # Advanced usage # # 'identifyLines' allows you to determine how the different steps # of the identification must be performed through several optional # parameters. # Before introducing them it's necessary to understand what 'identifyLines' # does. # # The line identification is a 2-step process: # 1. the lines are detected # - In order to optimize the quality and the speed of the detection, the # spectrum is split into N smaller spectra ('nbOfSplits') with a given # 'overlap' # - The detection is performed by SmoothBaseline that needs a 'box' value # for the smoothing and a 'midcyle' value. (cf SmoothBaseline documentation) # - The detected regions of lines are then separated using a boxcar # smoothing ('smoothBox') and a slope detection algorithm. # - All the lines that are detected may not be real lines but noise, so # a filter using a minimal signal-to-noise ratio ('minSNR') is necessary to # keep only relevant lines. The signal-to-noise ratio is computed with an # estimation of the FWHM of the lines ('fwhm'), the 'calibration' value of # the instrument and the 'width' of the region around the line used to # compute a local RMS. # # 2. these detected lines are then compared with the known lines ('linelist') # - Basically, each detected line is compared with each known line and if # the frequency between the detected line and the known line is inferior # to a given 'tolerance' the detected line is identified otherwise it is # unidentified. # For the sake of completeness, here are the units of the parameters: # - midcycle: per inverse wavenumber in micron # - fwhm: km/s # - width: GHz # - tolerance: km/s # - overlap: GHz # # For more details about all the parameters, type: print identifyLines.__doc__ # or read the corresponding entry in the User Reference Manual. # ----------------------------------------------------------------------------- # Now that the optional parameters are clearer to you, let's rewrite the # previous command line with the default values of the optional parameters # so you have an overview of what you can do: lines = identifyLines(data=spectrum, linelist=linelistFile, vlsr=vlsr) # is actually lines = identifyLines(data=spectrum, linelist=linelistFile, vlsr=vlsr, nbOfSplits=4, overlap=0.25, # the spectrum is split box=200, midcycle=1.7E6, # the regions of lines are detected smoothBox=7 # the regions of lines are separated into lines fwhmEmission=3.0, fwhmAbsorption=1.0, minSNR=5.0, width=0.05, # the lines with a minimum signal-to-noise ratio are kept toleranceEmission=4.5, toleranceAbsorption=1.5, # the lines are compared to the lines in the linelist ) # You get the identified, unidentified lines and the baselined spectrum identified = lines["identified"] unidentified = lines["unidentified"] baselinedSpectrum = lines["baseline"] # That's it! The identification is finished. Identifying your lines is one # thing, analyzing what was identified is another story. Fortunately, # the line identification tool comes with a few functions that may be helpful. # ----------------------------------------------------------------------------- # A guided tour # ----------------------------------------------------------------------------- # From here you should read the comments and run the code line by line # to make sure you understand everything. # Getting to know your lines # -------------------------- # Let's look at the identified lines. The variable 'identified' is a Linelist # (a data structure containing Line objects) print identified # You should see the number and the species of the lines. If you want to print # more details for each line, use the inspect method: identified.inspect() # Note that 'inspect' doesn't return anything, it just prints the result. # Sometimes a table is more handy: table = identified.toTable() # If you open 'table' in the Dataset viewer you'll see the same result as with # 'inspect' but this time each line is in a row, and each attribute of a line # is in a column. If you want to get, for example, the sky frequency of all the # lines, you can do it with: print table["skyFrequency"].data # But you can also do it with: print identified.get("skyFrequency") # Working with the Linelist has other advantages. We can select the lines with # a specific attribute. For example, this is how we select all the A-CH3OH # lines: print identified.select(species="A-CH3OH") # This command returned a Linelist, which means you can apply all the methods # you know so far, including 'select'. ach3oh = identified.select(species="A-CH3OH").select(maxAmplitude=[0.7, 0.9]) print ach3oh # You've just selected all the A-CH3OH lines whose maximum amplitude is between 0.7 and # 0.9 K. You can check it's true with: print ach3oh.get("maxAmplitude") # If you want more details about 'select', type: print identified.select.__doc__ # Saving and plotting # ------------------- # If you're happy with the lines you have (I assume you are), let's save them # in a FITS file so that we can use them later. ach3oh.saveLinesTo("ach3oh_lines.fits") identified.saveLinesTo("identified_lines.fits") # 'saveLinesTo' saves all the lines in a TableDataset into a FITS file or a CSV # file (use .csv instead of .fits if you prefer CSV). # So far you just used numbers and tables. You can have a spectrum of the # lines by using the 'spectrum' method: ach3ohSpectrum = ach3oh.spectrum(data=baselinedSpectrum) # or just plot them with 'plot': ach3oh.plot(data=baselinedSpectrum) # Note: You need the baselined spectrum (or at least a spectrum) for 'spectrum' # and 'plot'. We'll explain why later as it involves some details on how the # Line was implemented. # The same way, we used 'saveLinesTo' to save our lines to a file, you can save # your spectrum to a file with 'saveSpectrumTo': ach3oh.saveSpectrumTo("ach3oh_spectrum.fits", spectrum=baselinedSpectrum) # You probably noticed that the spectrum created is just a continuum at 0 K # except where the lines were detected. Let's choose the first and see what we # can learn. # Looking closer # -------------- # I said earlier that a Linelist contains Line objects. # A Line just stores information about the line detected in your spectrum. # In order to get a line, specify its index (starting from 0) print ach3oh[0] # In fact, you can treat the Linelist like any other iterable you're used to # iterate on it with a for loop: for line in ach3oh: print line # Like for a Linelist, if you want to know in details the attributes of this # line, just use 'inspect': print ach3oh[0].inspect() # All these attributes are easily accessible. Let's say you want to get the # signal-to-noise ratio (snr) of this line, just specify this attribute: print ach3oh[0].snr # Linelist and Line share some methods like 'plot' and 'spectrum', # So you can plot this line with the 'plot' method. ach3oh[0].plot(data=baselinedSpectrum) # and get the corresponding spectrum with the 'spectrum' method: ach3ohSpectrum = ach3oh[0].spectrum(data=baselinedSpectrum) # One remark about these 2 methods and the philosophy behind the Line object: # As you saw a Line object has just information about the line in the spectrum, # not the spectrum itself. The attribute allowing to bind the line and the # spectrum is the frequency range. This is why the methods 'plot' and # 'spectrum' need the spectrum as argument in order to select only the part of # the spectrum within that frequency range. # Adding and removing lines # ------------------------- # Any convenient data structure allows you to add or remove elements. # So does the Linelist. # First, let's create an empty Linelist: myLinelist = Linelist() # You can add a Line or a Linelist into this Linewith with 'add': # Add all the A-CH3OH lines and the H2CS line: ach3ohAll = identified.select(species="A-CH3OH") h2cs = identified.select(species="H2CS") myLinelist.add(ach3ohAll) myLinelist.add(h2cs) # Check your Linelist, there should be 5 lines in it: print myLinelist print len(myLinelist) # We could have had the same result by removing the CS, v=0-4 line from our # initial Linelist with method 'remove': cs_v0_4 = identified.select(species="CS, v=0-4") identified.remove(cs_v0_4) # Note: you could also have removed it with its index (i.e 5): # identified.remove(5) # Be careful with 'remove' because it removes permanently the lines from the # Linelist. # If you want to make a copy of your Linelist, you can do it by creating a new # Linelist with the lines you instead of letting it empty: identifiedCopy = Linelist(identified) print identifiedCopy # We remove all the A-CH3OH lines identified.remove(ach3ohAll) # We print the result print "The original linelist has %d lines." % len(identified) print "The copied linelist has %d lines." % len(identifiedCopy) # identified doesn't have A-CH3OH lines anymore and identifiedCopy isn't modified. # Loading your lines # ------------------ # It's generally a good idea to save your lines so you can reload them if # you want to share them or if something goes wrong. Fortunately, that's what # we did. We can load all the lines we had at the beginning with 'loadLinesFrom'. # First we create an empty that will be populated by the lines that we'll load # from the file "identified_lines.fits": originalLines = Linelist() originalLines.loadLinesFrom("identified_lines.fits")