# mapping demo script # prepared for HIFI PSP/SDP DP workshop April 2010 # # v0.0 CM: created # #------------------------------------------------------------------------------ # Spectral cubes are created at the end of the level 2 pipeline using the # doGridding task. # # This script demonstrates how to create a cube yourself, using options of the # doGridding task. After, some aspects of the CubeSpectrumAnalysisToolbox (CAT), # which is used to visualise and analyse spectral cubes, are demonstrated. # # The doGridding task is described in chapter 10 of the HIFI User manual. You # can also inspect the task parameters with # print doGridding # or # print doGridding.__doc__ # # The CAT is covered in Chapter 5.2.6 of the Data Analysis Guide. # # DATA: # obsid 1342178898 # HIFI Mapping mode DBS Raster - data used in DR 21 press release # # Data processed with HIPE 1.0.6, so level 2 data not split into USB and LSB, # and are in MHz rather than GHz. If you run doGridding on this data, you will # see many "WARNING CubeBuilder: could not add parameter sideband" messages # because of this, as data processed by the 2.6 pipeline will be split into LSB # and USB. # # Contains 1101 GHz 13CO(10-9) [WBS subband 3] # and 1113 GHz H2O(1_11-0_00) [WBS subband 2] #------------------------------------------------------------------------------- # ################################################################################ # 1. Get data into session ################################################################################ # obs = getObservation("1342178898", poolName="pool_1342178898") # # Like all pipeline steps, when using the doGridding task in interactive mode # you should supply it with an HTP (HifiTimelineProduct) rather than an # ObservationContext. # Furthermore, it should be a level 2 HTP so that the spectra have been # resampled to a linear frequency axis (doFreqGrid in the Level-2 pipeline) # htp = obs.refs["level2"].product.refs["WBS-H"].product # ################################################################################ # 2. run DoGridding with default parameters to make cubes for each subband ################################################################################ # cubes = doGridding(htp=htp) # # To do this in the GUI #select a HifiTimelineProduct variable in the 'Variables' view, # then in the 'Tasks' view go to : 'Applicable' -> 'doGridding' # and double click to popup the GUI of the task of course also under: # 'Tasks' -> 'All' -> 'doGridding' # 'Tasks' -> 'By Category' -> 'HIFI' -> 'doGridding' # # or open the GUI from the command line with # openTask("doGridding") # # If you opened the GUI after selecting an htp in the Variable view (otherwise # drag htp to the "htp" bullet), then clicking "Accept" will run the task # with default paramters. # # Nice GUI feature: if you hover on the parameter names you will see a # tip / descriptive hint for the parameter # # The output of doGridding is an array of cubes, one for each subband in the HTP. # Assuming there are 4 subbands in the data, each cube can be extracted by: # cube_1 = cubes[0] cube_2 = cubes[1] cube_3 = cubes[2] cube_4 = cubes[3] # # Please note that in HIPE 3.x that a secondary output will be a MapContext # containing all the cubes so that it can be browsed with the HIPE contexts browser. # # If you are creating scripts to do this, you might like to do this automatically # for subband in range(len(cubes)): cube = cubes[subband] subband = cube.meta['subband'].value cube_name = "cube_%d" % subband vars()[cube_name] = cube # Note that this code snippet just overwrites the cube_1 etc produced above, so # you only see the new variables "cube", "subband", "cube_name" # Or even automatically name each of the cubes something meaningful, such as by # obsid, spectrometer (apid), and subband def getCubeName(cube): """ Build a meaningful name for the cube """ subband = cube.meta['subband'].value obsid = cube.meta['obsid'].value apid = cube.meta['apid'].value apidName = herschel.hifi.pipeline.util.PipelineUtils.getApidName(apid) apidUpperCase = String.toUpperCase(apidName) cubeName = "cube_%d_%s_subband_%d" % (obsid, apidUpperCase, subband) return cubeName # "def" here will define getCubeName, which is used below. Highlight the # section of text above from "def" to "return cubeName" and then run it all with # the single green arrow before stepping through the next part, which will # get a separate variable for each cube computed for each subband for subband in range(len(cubes)): cube = cubes[subband] cube_name = getCubeName(cube) vars()[cube_name] = cube # ################################################################################ # 3. "Basic" doGridding # The "basic" usage of doGridding allows to define the: # subbands - useful to save memory if there is only a line in one subband # the beam size (half power beam width) # the weighting used in the convolution into a pixel # the filter used in the convolution - gaussian or box # the parameters of the filter ################################################################################ # # ---Make a cube for subbands 2 and 3 only--- cubes = doGridding(htp=htp, subbands=Int1d([2, 3])) # # If you prefer the GUI, make a variable, e.g, subbands = Int1d([2, 3]) and # drag tha variable into the "subbands" bullet. # # ---Specify a beam size (in arcsec)--- # At present, the beam size is calculated as 75.44726* wavelength [mm]. It is # not expected to be inaccurate but you may wish to change this cubes = doGridding(htp=htp, beam = Double1d([20.6])) # even to an elliptical beam (x,y) cubes = doGridding(htp=htp, beam = Double1d([10.0, 20.0])) # # Again, to use this in the GUI you must create the appropriate variable # and drag it to the "beam" bullet. # #---Change weighting--- # The spectra can be equally weighted during convolution, "equal", this is default. # or given the weighting assigned during pipelining with "selection" cubes = doGridding(htp=htp,weightMode="selection") # # In the GUI use the drop down menu # #---Change the convolution filter parameters--- # By default, a Gaussian ("gaussian") filter is used for convolution. This may # not be appropriate for a raster map and a "box" filter could be preferred. # It is also possible to modify the parameters that characterise the filters: # length - x, y - of the pixel, measured in pixels,for a box filter xFilterParameters = Double1d([0.5]) yFilterParameters = Double1d([1.5]) cubes = doGridding(htp=htp,weightMode="equal",filterType="box",xFilterParams=xFilterParameters, yFilterParams=yFilterParameters) # or, more compactly parameters = [Double1d([0.5]), Double1d([1.5])] cubes = doGridding(htp=htp,weightMode="equal",filterType="box",filterParams=parameters) # # The next example specifies the parameters length and sigma of the Gaussian # filter function, when using a Gaussian filter (default case). # The "influence area" is the area surrounding a grid point # where the algorithm must pick up all the available data points. influence_area = 1.95 # length in pixels # sigma of the gaussian function times SQRT(2)# sigma_sqrt2 = 0.3 # in pixels xFilterParameters = Double1d([influence_area, sigma_sqrt2]) # default case: influence_area = 1.8; sigma_sqrt2 = 0.36 yFilterParameters = Double1d([influence_area, sigma_sqrt2]) # cubes = doGridding(htp=htp,filterType="gaussian", xFilterParams=xFilterParameters, yFilterParams =yFilterParameters) # # it is again possible to pass both set of parameters in a single input: parameters = [Double1d([1.8,0.4]), Double1d([1.6,0.3])] cubes = doGridding(htp=htp,weightMode="equal",filterType="gaussian",filterParams=parameters) # # In the GUI, the filter type can be selected from the drop down menu, the # filter parameter variables should be defined as above and dragged into the # x- and yFilterParams bullets. # # ################################################################################ # 4. "Expert" doGridding # # There are many additional functions available to further refine the behaviour # of doGridding. # If you prefer to use the GUI then you must define a variable as below and # drag that variable into the appropriate bullet. # ################################################################################ # #---datasetTypes--- # By default the DoGridding task reads all the datasets of the HTP that are # flagged in the summary table with isLine to true ("on" spectra) and whose # sds_type correspond to "science" datasets. In the case that you wish to # change the default, you can specify yourself which type of dataset should be read. types = String1d(["science","other"]) cubes = doGridding(htp=htp, datasetTypes=types) # # By default, doGridding ignores any science data set that is not the line # (i.e. an OFF or reference), to change that behaviour cubes = doGridding(htp=htp, ignoreOFFs=False) # #---smoothFactor--- # The smoothFactor can be used to extend the area of influence of each pixel, # i.e. the kernel length. # This factor is given by means of a Double1d input with one or two elements, # the smoothing factor along each of the x and y axis. If both are equal, # a Double1d with a single element can be provided. The following example, where # the pixel size is set to 15"x 25", uses an smoothing factor equal to 1.5 # along the x-axis (horizontal axis of the grid, RA axis) and 2.0 along the # vertical axis of the grid. # pixelSize = Double1d([15,25]) cubes = doGridding(htp=htp,weightMode="selection",smoothFactor=Double1d([1.5,2.]),pixelSize=pixelSize) # #---Pixel size--- # Pixel size must be defined in arcsec. Here a pixel size of 20" square is defined: cubes = doGridding(htp=htp,weightMode="selection",filterType="gaussian",pixelSize=Double1d([15])) # # If the beam size is specified, and the pixel size is not specified, the pixel # size will be function of the beam size taking into account the Nyquist # criterion and the smooth factor (if any given). Usually, for nyquist sampling, # the default pixel size becomes half the beam size # #---mapSize--- # Define the map size in pixels (10 by 20) cubes = doGridding(htp=htp, mapSize=Int1d([10,20])) # #---ref pixel and refPixelCoordinates--- # Treat this with caution: it is advised to use the defaults unless you know # well what your data is and what you are doing with doGridding. # # Specify which pixel is the reference pixel (0,0) of the map. By default it is # the centre of the bottom-most, left-most pixel in the map. cubes = doGridding(htp=htp,refPixel = Double1d([3.5, 4.0])) # The coordinates of this pixel will be calculated, or you can provide them yourself refPixel = Double1d([0,0]) longitude = 307.9 latitude = 40.36 refPixelCoordinates = Double1d([longitude, latitude]) cubes = doGridding(htp=htp,refPixel=Double1d([0,0]), refPixelCoordinates=Double1d([longitude, latitude])) # or: cubes = doGridding(htp=htp,refPixel=refPixel, refPixelCoordinates=refPixelCoordinates) # #---detail--- # #The user maybe interested on getting more details about the process performed # by the gridding task. For that purpose, there is an optional, boolean input # parameter called "detail". If this is set to true, like in the next example, # the user will get an output table dataset called "convolutionTable". This # provides a table that shows which spectra have contributed to each pixel, # with which weight, etc. In addition the header of this table informs about # the filters used (type, parameters) and about the regular grid (pixel size, # beam size, reference pixel, grid dimensions, etc.). # cubes = doGridding(htp=htp,detail=True) convolutionDataset = doGridding.convolutionTable # or cubes = doGridding(htp=htp,detail=Boolean.TRUE) convolutionDataset = doGridding.convolutionTable # # or clicking on the "detail" checkbox of the "expert" GUI # #---offsetsTable--- # Other secondary outputs can be obtained from the task, e.g. the offsets # (measured in radians wrt the projection centre) of each point of the # regular grid along the x and y axis of the grid # xPoints = doGridding.xPoints yPoints = doGridding.yPoints print xPoints #