# MapHandsOn.py: script containing examples of use of doGridding and # the CubeSpectrumAnalysisToolbox # # v0.0 Jan 2011: CM created for NHSC OT-1 workshop # #------------------------------------------------------------------------------- # Data available (see hand-out for more information) # 1342210097 in pool 1342210097_map_otf_1b_halfbeam, 8x6 half-beam sampled (19.5"x19.2") OTF map in band 1b # 1342210099 in pool 1342210099_map_otf_1b, 8x7 Nyquist sampled (16.0"x16.2") OTF map in band 1b # 1342205474 in pool 1342205474_map_dbs_1b, 3x3 DBS Raster map in band 1b, with 4.5'x4.0' spacing # 1342205481 in pool 1342205481_map_dbs_7b, 3x3 DBS Raster map in band 7b, with 4.5'x4.0' spacing # #------------------------------------------------------------------------------- # #------------------------ EG 1 ------------------------------------------------- # Example 1: re-grid 1342210099 with the pixel size corresponding to the spacings # used in the observation. # # # Useful documentation: # Chapter 16, "How to make a spectral cube", of the HIFI User Manual # # Load data obs_map_otf_1b_nyquist = getObservation("1342210099", poolName="1342210099_map_otf_1b") # # Gridding, like all pipeline tasks, is performed on the HifiTimelineProduct (HTP). # For doGridding, Level 2 HTP must be used. Extract the WBS-H-USB HTP # htp_otf_nyquist_wbs_h = obs_map_otf_1b_nyquist.refs["level2"].product.refs["WBS-H-USB"].product # # # By plotting the HTP in SpectrumExplorer one can see that the strongest line # falls in subband 1. To reduce computing time, create a cube for subband 1 only. # The subband is specified using "subbands" and the value must be a 1d integer # # The pixel size (here 16.0"x16.2", see crib sheet) is set using "pixelSize" # and the beam size (here 37.7", see crib sheet) is set using "beam". # Both are Double precision 1d values. # # The output of doGridding is an array of cubes (in this case the array contains # only one cube). # cubeArray_otf_nyquist_wbs_h = doGridding(htp=htp_otf_nyquist_wbs_h, \ subbands=Int1d([1]), pixelSize=Double1d([16.0,16.2]), beam = Double1d([37.7])) # # The array contents are numbered from zero, to get the first cube: # cube_otf_nyquist_wbs_h = cubeArray_otf_nyquist_wbs_h[0] # # You can view the regridded cube in the SpectrumExplorer or the # CubeSpectrumAnalysisToolbox by right clicking on the variable # cube_otf_nyquist_wbs_h in the Variables View and selecting "Open With # SpectrumExplorer" or "Open With CubeAnalysisToolbox" # #------------------------------------------------------------------------------- # #------------------------ EG 2 ------------------------------------------------- # Example 2: compare the signal in a fully Nyquist sampled map (1342210099) with # that in a half-beam sampled (1342210097) map. # # Useful documentation: # Chapter 16, "How to make a spectral cube", of the HIFI User Manual # Sections 5.7-5.10, of the Herschel Data Analysis Guide but especially # 5.9.4.3. "Cube Comparison" and 5.9.2.2. "Multiple Spaxel Display" # # Repeat exercise one for 1342210097 # # Load data and extract HTP obs_map_otf_1b_halfbeam = getObservation("1342210097", poolName="1342210097_map_otf_1b_halfbeam") htp_otf_halfbeam_wbs_h = obs_map_otf_1b_halfbeam.refs["level2"].product.refs["WBS-H-USB"].product # # With the exception of the beam spacing, this observation is set up identically # to 1342210099 so the line also falls in subband 1. # # The pixel size (here 19.5"x19.2", see crib sheet) # cubeArray_otf_halfbeam_wbs_h = doGridding(htp=htp_otf_halfbeam_wbs_h, \ subbands=Int1d([1]), pixelSize=Double1d([19.5,19.2]), beam = Double1d([37.7])) cube_otf_halfbeam_wbs_h = cubeArray_otf_halfbeam_wbs_h[0] # # Now to compare the cubes: # # This can be done visually point-by-point in the CSAT using the cube comparison tool. # Open cube_otf_nyquist_wbs_h with CSAT as above. # (You might find it easiest to undock the CSAT by clicking on the tab at the top # of the editor (avoid the x!) and dragging it outside of HIPE. # Alternatively, maximize HIPE on your screen.) # # Select "cube comparison" from the Analysis menu, the cube comparison working # area appears to the right. Select cube_otf_halfbeam_wbs_h from the drop-down # menu at the top of the GUI and click on the load Cube button. # (Note that you can also load fits files for comparison with the load Product # button below) # # In the region beneath, cube_otf_halfbeam_wbs_h appears. You will need to # zoom the image to fit the screen using the zoom in/out buttons beneath. There # will be a message warning that the two maps do not have the same spatial extent # because the half-beam map is slightly larger, however, the maps centres should # be coincident. # # Move the mouse cursor over the map on the right hand side of the GUI, in the # live spectrum display at the top left of the GUI the spectrum for that pixel # (or the closest aligned) from the nyquist sampled cube (on the left) is # displayed in dark blue and the spectrum from the half-beam cube (loaded in on # the right is green. # # Spectra from individual pixels of the two maps can be extracted for comparison # too using the Single Spaxel Extractions task under the Spectral Extraction menu. # However, there is a bug and HIPE freezes if you try to save extracted # spectra to HIPE (using Save as Product). # You can save the result to script (remember to specify .py at the end of the # file name in order to be able to read the script in HIPE again and this will # show you how the extraction can be done via the command line: # nyquist_pix_3_3 = extractSinglePixelSpectrum(simplecube=cube_otf_nyquist_wbs_h,\ posX=3,posY=3) # # The above command extracts pixel (3,3) and creates a spectrum that can be viewed # in SpectrumExplorer and analysed with the spectrum (arithmetic) toolbox. The # pixels in cubes are numbered from (0,0) with that pixel being the bottom-most, # left-most pixel in the cube. The pixel numbers are displayed immediately # beneath the cube. # # The global averages of the maps can also be extracted and compared. # # In the GUI, use the Area Extraction task found under the Spectral Extraction # menu. The Area Extraction tab appears in the right hand side of the GUI. # Select whole cube (other extractions are also possible) and press Show Spectrum. # The map average spectrum is shown below. This can be safely sent to HIPE as a # product "using save a Product" and can be found in the Variable View. In the # case of the nyquist cube created above, the automatically created name of the # new variable will be # cube_otf_nyquist_wbs_h_regionspectrum_HIFI_cube_product__globalspectrum_ # # Repeat the procedure for the other cube. To compare the averages of the two # maps plot one of the average spectra in SpectrumExplorer then drag the variable # name of the other average spectrum into SpectrumExplorer to overplot. # # This can also be done in the command line as # av_spectrum = extractRegionPixelSpectrum(simplecube=cube_otf_nyquist_wbs_h,wholeImg=True) spectrum_otf_nyquist_wbs_h_av = extractRegionPixelSpectrum.finalspectrum # The task now needs to be reinitialised. extractRegionPixelSpectrum = ExtractRegionPixelSpectrumTask() av_spectrum = extractRegionPixelSpectrum(simplecube=cube_otf_halfbeam_wbs_h,wholeImg=True) spectrum_otf_halfbeam_wbs_h_av = extractRegionPixelSpectrum.finalspectrum # # The two average spectra can be overplotted in SpectrumExplorer as above. # You can notice how much signal is lost in the half beam map. # #------------------------------------------------------------------------------- # #------------------------ EG 3 ------------------------------------------------- # # Exercise 3: Make an integrated map for 1342210097 and then overlay contours # on the resulting image. Compare the result obtained with that if baseline # correction is done prior to gridding. # # Open cube_otf_halfbeam_wbs_h with the CSAT. # Select the Integrated map task from the Cube Manipulation menu. # Define the integration region on the global spectrum that appears to the right # (Mac users should right click on the plot and use SelectRange from the Tool menu, # the region can then be selected with a click and drag) # Press Perform Integration # # Wait... # # The integrated map appears to the left of the global spectrum, it will need to # be zoomed to fit using the zoom in/out buttons below # Press Save a Product(s) to create the image variable in HIPE. In this case # the name will be similar to # cube_otf_halfbeam_wbs_h_integratedImage_freq_576_33250_Width_0_25100 # The numbers indicate the region you selected. If several ranges were selected # you would create an image per range. # # In the command line, an integrated map is created as follows: # sra = Double1d([1400.0]) # The channel number of the start of the integration region (read it off the plot) era = Double1d([2100.0]) # The channel number of the end of the integration region num_ranges = 1 # The number of selected regions you want integrated_map =integrateMapFromCube(cube= cube_otf_halfbeam_wbs_h, \ nbStack=num_ranges,startArray=sra,endArray=era) # This is an array of images (one per selected region) image = integrated_map[0] # Get the first image in the array # # Now plot contours over this image. The easiest way to do it is to use the # automaticContour task # Click on your image in the Variables view, automaticContour will appear in the # Applicable folder of the Tasks view. Double click on it to open. # enter the number of contour levels you would like into the GUI, say 10, and # click on accept. # # In the command line contours = automaticContour(image=image,levels=10,min=-6.085906215198587E-13, \ max=4.000805955319528E-12,distribution=1) # # The min and max above are just the defaults from the GUI. You have three # choices for the distribution: Linear = 0, Log=1, Ln=2 # # Now open your image with the Standard Image Viewer (an option on right clicking # on the variable name) and then drag your contour variable from the Variables # view onto the image. # Keep this image open in a tab in the Editor view. # # Now compare the results after baseline correction # # Correct the baselines of the datasets in the HTP, here automatic masking is # used, you can set domask=1 if you prefer to do the masking yourself. fit = fitBaseline(htp=htp_otf_halfbeam_wbs_h,order=2,domask=2,movemasks=2, \ doglue=False,doreuse=True,restart=True,basemode="subtract",smoothResolution=2.0,\ interactive=False,htpcopy=False) htp_otf_halfbeam_wbs_h = fitBaseline.result_htp # # Note that the example here is for one HTP but if you intend to create cubes # for both spectrometers and polarizations it is more efficient to pass the # ObservationContext to fitBaseline. # # Now regrid the data as before # cubeArray_otf_halfbeam_wbs_h = doGridding(htp=htp_otf_halfbeam_wbs_h, \ subbands=Int1d([1]), pixelSize=Double1d([19.5,19.2]), beam = Double1d([37.7])) cube_otf_halfbeam_wbs_h_fitted = cubeArray_otf_halfbeam_wbs_h[0] # # Make a new integrated map sra = Double1d([1400.0]) era = Double1d([2100.0]) num_ranges = 1 integrated_map =integrateMapFromCube(cube= cube_otf_halfbeam_wbs_h_fitted, \ nbStack=num_ranges,startArray=sra,endArray=era) image_baseline_fitted = integrated_map[0] # # Now plot contours over this new image. # contours_baseline_fitted = automaticContour(image=image_baseline_fitted,levels=10,min=-6.085906215198587E-13, \ max=4.000805955319528E-12,distribution=1) # # Open the new image and drag the new contours into it as before. # # # In the case that you wish to more quickly practice creating an integrated map, # you might prefer to use one of the smaller DBS raster maps. # Below the cube is regridded using appropriate parameters and then the # integrated map is created. obs_map_dbs_7b = getObservation("1342205481", poolName="1342205481_map_dbs_7b") htp_dbs_7b_wbs_h = obs_map_dbs_7b.refs["level2"].product.refs["WBS-H-USB"].product # # Chop modes suffer less from baseline drift than total power modes but it is # still a good idea to correct baselines. # fit = fitBaseline(htp=htp_dbs_7b_wbs_h,order=2,domask=2,movemasks=2,doglue=False,doreuse=True,restart=True,basemode="subtract",smoothResolution=2.0,interactive=False,htpcopy=False) htp_dbs_7b_wbs_h = fitBaseline.result_htp # # The pixel size is 4.5"x4.0", the beam size is 12.2", and the line falls in subband 4, see crib sheet) # cubeArray_dbs_7b_wbs_h = doGridding(htp=htp_dbs_7b_wbs_h, subbands=Int1d([4]), \ pixelSize=Double1d([4.5,4.0]), beam = Double1d([12.2])) cube_dbs_7b_wbs_h = cubeArray_dbs_7b_wbs_h[0] # sra = Double1d([1481.0]) era = Double1d([1982.0]) num_ranges = 1 integrated_map =integrateMapFromCube(cube= cube_dbs_7b_wbs_h, \ nbStack=num_ranges,startArray=sra,endArray=era) image = integrated_map[0] # #------------------------------------------------------------------------------- # #------------------------ EG 4 ------------------------------------------------- # Exercise 4: Make a line intensity map # # You can also use the Line Intensity Map task of the CSAT to make a continuum # subtracted map of the integrated flux of a line. This is useful if the # baselines in your data only require a simple correction. # # Open a cube in the CSAT and Select the Line Intensity Map task from the # Analysis menu. # Select the continuum subtract mode and degree of Polynomial, in the script # below two ranges either side of the line seen in cube_otf_nyquist_wbs_h are # selected to have the continuum fitted. # Select the degree of Polynomial that will be used to fit to the contiuum. # Select the model to be used to fit the the line, below a First order Polynomial # and a Gaussian are used. # Press Integrate Intensity # # In the command line: # startarray = Double1d([1100.0]) # the start of the first region used to fit to the continuum endarray = Double1d([2150.0]) # the end of the second region # These arrays are also given in terms of channel number, which can be found # by reading from the information that appears under the plot in the right hand # pane when the mouse is hovered over the spectrum. # LineIntMap= lineIntensityMap(cube=cube_otf_nyquist_wbs_h,continuumSubtraction=1,\ PolyDegree= 1,fitModel= 'Gaussian automatic',typeOfLine= 'emission', \ regionStartArray=startarray,regionEndArray=endarray) # #------------------------------------------------------------------------------- # # These exercises have focussed mostly on the WBS-H spectra of the OTF maps # in order to keep things uncluttered for the purposes of demonstration. If you # have time you can repeat these exercises for WBS-V and/or the HRS data.