15.3. Using doGridding...

15.3.1. ...to change beam, pixel, and map size

  • Mapping observations from OD 835 onwards, and processed with HIPE 8 onwards, are gridded using the scan and readout spacing used during the observation, whilst earlier observations are gridded using a better approach than in HIPE 7. Note that the data must be processed including the doUplink step from the Level 0 pipeline so you must either retrieve data from the HSA processed already with HIPE 8 onwards, or you should run the script given in the introduction to this chapter, where more information is given on this topic, Section 15.1.

    Despite the need to grid data according to your original specifications, you may wish to regrid data with a new pixel or beam size in order to compare with other observations. The beam width in arcseconds is automatically entered into the GUI when you pass the HTP to the task. It can be found in the beam box and can be edited there. If you modify the beam size then the filter parameters will also be automatically updated unless you specify parameter values yourself, see Section 15.3.6 below.

    For observations from OD 835 onwards, the beam width is calculated from the calibration tree using the frequency requested in HSpot. Note that this results in slightly different beam sizes for the WBS and the HRS as a consequence of how the spectrometers are set up. For older observations, the beam size is calculated from HPBW["] = 75.44726 * wavelength[mm] where wavelength is the LO Frequency.

    The beam width can be set in the command line with:

    # The extended version which can be adapted to the rest of this section
    cubesContext, cubes, cube, xPoints, yPoints, convolutionTable, grid = doGridding(htp=htp,beam=[15.4])
    
    # or the shorther version, where, here, we are only interested in obtaining the cubesContext
    cubesContext = doGridding(htp=htp, beam=[15.4])[0]
    

    The HIFI beam is symmetric but you can also specify an elliptical beam by specifying the dimensions of the x and y axes:

    # specify an elliptical beam. In this case the beam is wider along the vertical axis.
    cubesContext = doGridding(htp=htp, beam=[20.0, 40.0])[0]
    
  • For observations from OD 835 onwards, the doGridding task will compute the pixel size to be appropriate according to the spacing used in the observation - metadata values of mapReadoutSep and mapLineStep - these reflect the angular separation of the scan lines as carried on the sky at the sampling requested with the observation. For earlier maps, i.e. when mapReadoutSep and mapLineStep are not present in the product metadata, different estimates are used depending on what information can be found in the metadata, such as the structure and number of datasets within the HTP, and the metadata parameters supersampling, n_cycles for OTF maps, and crossStep for DBS raster/cross maps. In the most sparse information case, the task will assume the beam size is given by 75.44726 * wavelength[mm] (where wavelength is the LO Frequency) and take (beam width/2) for half-beam spacing maps and (beam width/2.4) for Nyquist sampled maps.

    You can specify the pixel size by editing the pixelSize boxes in the GUI, units are arcseconds. Note that if you do this then the filter parameters will also be updated to reflect the new pixel size. However, editing the filter parameters does not affect the pixelSize parameter. Editing the beam width parameter does not affect the pixelSize parameter but entering a non-zero flyAngle (see Section 15.3.5) will in order that flux is conserved when rotating the map.

    The pixelSize can be modified in the command line as:

    # For a square pixel
    cubesContext = doGridding(htp=htp, pixelSize=[16.36])[0]
    # For a rectangular pixel
    cubesContext = doGridding(htp=htp, pixelSize=[10.0, 20.0])[0]
    
    [Warning] Warning

    Flux Conservation in Spectral Cubes from Mapping Observations:

    The doGridding task in the HIFI pipeline is responsible for convolving the spectral datasets acquired in an OTF or DBS Raster mapping observation into a spectral cube with a specified pixel scale, which by default should match how the scan lines and readout points within each line were spaced during the observation. The scheme of signal filtering and interpolation to put the data on the specified grid may affect the overall flux conservation, at a level which is low but you should be aware of. For example, the total signal summed from a spectral cube produced using a Gaussian filter over the datasets of a Nyquist-sampled OTF map is generally < 1% lower than the sum of the signal taken directly from the input datasets (the Level 2 HTP datasets) before convolution. A part, or all of this slight mismatch may be on the assumed versus actual beam shape at the observed frequency. If you wish to put the map on a coarser grid, effectively reducing the spatial resolution to a wider beam in order to match another observation, then the flux losses become more noticeable. Doubling the pixels sizes from their default (native map point spacing) can reduce the total flux by as much as 10%, accompanied by an increase in baseline RMS noise. The effect is more prevalent in OTF maps than DBS Raster, and in addition to deviations from ideal beam shape characteristics becoming more important, the filtering and interpolation method, and parameter values can be influential. No changes to the doGridding algorithm are planned, and this issue applies to all HIPE versions.

    For convenience, the beam widths per HIFI band and beam spacing for half-beam and Nyquist sampled maps are given in the table below.

    Band Beam Width (") Nyquist spacing (") Half-beam spacing (")
    1a 43.5 18.4 22.0
    1b 37.7 15.9 19.0
    2a 33.3 13.9 17.0
    2b 29.8 12.5 15.0
    3a 27.3 11.4 14.0
    3b 24.9 10.4 13.0
    4a 22.5 9.4 11.0
    4b 20.8 8.7 10.0
    5a 19.6 8.0 9.0
    5b 18.6 7.8 9.0
    6a 15.2 6.3 8.0
    6b 14.0 5.8 7.0
    7a 12.8 5.3 7.0
    7b 12.2 5.2 6.0
  • The pixel size is multiplied by a smoothing factor, which is 1 by default. You can modify the smoothing factor by editing the smoothFactor box in the GUI or in the command line by:

    cubesContext = doGridding(htp=htp, smoothFactor=[2.0, 2.0])[0]
    

    The convolution filter parameters are dependent on the smoothing factor and modifications to it will cause the filter parameters to be automatically updated, see Section 15.3.6 below.

  • By default, the mapSize parameter is set to 0 in order to let the algorithm choose the optimal map dimension. The map size is automatically calculated by doGridding based on the spatial information in the input spectra. The map size can be specified in units of pixels by entering values in the mapSize boxes on the GUI. If you wish to fix the map size, you must set mapSize to your choice (say [8.0, 7.0]), and then you must set pixelSize =(0.0, 0.0) and (x,y)FilterParams =(0.0, 0.0).

    In the command line, this can be done as:

    cubesContext = doGridding(htp=htp, mapSize=[8.0, 7.0])[0]
    

    Note that if the map dimension in pixels are specified, then unless it is also specified, the pixel size will be given by the area observed divided by the number of pixels specified in the map size parameter.

Additional information about map size is available in the HIFI User's Reference Manual under doGriddingTask.

15.3.2. ...to make cubes of combined H- and V- polarisation

At the moment, the doGridding task only works with one HTP at a time. However, you may wish to create maps of combined H and V polarisation in order to increase signal to noise or to be able to compare with HSpot noise predictions, which assume that both polarisations are combined.

This can be done using the mergeHtps task, as follows:

# Get the HTPs of H and V data, here we use Level 2 WBS data from two polarisations of an observation (obs)
htpv=obs.refs["level2"].product.refs["WBS-V-USB"].product.copy()
htph=obs.refs["level2"].product.refs["WBS-H-USB"].product.copy()
#
# Merge the HTPS
htps = [htpv, htph]
mergedHtp = mergeHtps(htps=htps)
#
# Create the cube of merged data
cubesContext = doGridding(htp=mergedHtp)[0]

The resulting cube is labelled with the name of the first HTP passed to mergeHtps so in the example above, the resulting cube has the name cube_WBS_V_USB_1. Generally, the mergeHtps task attempts to merge metadata in a sensible fashion, and the details can be found in the mergeHtps URM entry.

Note that differences may be seen in H and V profiles, as a consequence of the beam separation between the H and V polarisations, and of structural and/or velocity variations in your source. If you are particularly interested in the spatial structure of your source you may prefer not to average the H and V polarisations together. More detailed information is provided in the HIFI AOT Observing Mode Release and Performance Notes v3.0 (24 September 2011), available from the HIFI instrument and calibration web pages.

The task mkRms allows you to assess the statistics of your new combined cube compared to the original cubes and HSpot noise predictions.

The mergeHtps task can be used to combine any number of HTPs and so can be used to create combined maps of data from different observations at the same frequency.

15.3.3. ...to make cubes for Solar System Objects

From HIPE 10 onwards, cubes of Solar System Objects (SSOs) produced by the standard pipeline are created with co-moving coordinates, i.e., the map centre follows the moving target. Cubes for non-SSO targets are automatically centred on the centre of the area observed. You can change these defaults using the comoving parameter, which has three possible options:

  • True: creates the map in comoving coordinates

  • False: creates the map centred on the centre of the area observed (appropriate for non-SSO objects)

  • None (default): for SSO, creates the map in comoving coordinates - for non-SSO, creates the map centred on the centre of the area observed

To, for example, force creation of a map centred on the centre of the area observed (this was previously the default approach) select comoving = False in the doGridding GUI. To do the same in the command line, use:

cubesContext = doGridding(htp=htp, comoving=False)[0]

You may wish to view the positions in your maps as offsets (or relative positions) rather than in the absolute coordinates that are offered in the Spectrum Explorer. This can be done using the doOffset task, see Chapter 7.

15.3.4. ...to make cubes more efficiently (limiting data input)

  • By default, doGridding will create cubes for each subband. However, if you know (from inspection of Level 2 HTPs) that you are only interested in a subset of subbands, you can specify these in the comma separated list in the subbands box of the GUI. This will reduce processing time and memory usage.

    A command line example to create cubes only for subbands 1 and 4:

    cubesContext = doGridding(htp=htp, subbands=Int1d([1, 4]))[0]
    #
    # Extract the cube for subband 1
    cube1 = cubesContext.refs["cube_WBS_V_USB_1"].product
    # Extract the cube for subband 4
    cube4 = cubesContext.refs["cube_WBS_V_USB_4"].product
    #
    

    The metadata of each cube will include a "subband" parameter stating the subband of the spectra which was used to compute the cube. This can be checked with:

    print cube1.meta['subband']
  • You can limit the spectral range of the cubes produced to cover only the region of interest by providing a range, per subband, for which cubes should be generated. The channel numbers must be input, and these can be read off the text to the bottom left of the plot in Spectrum Explorer while hovering the cursor over the spectrum.

    The channel ranges can be entered in to the channels boxes in the GUI, by default the full channel range is already entered. Channel ranges can only be entered for those subbands selected in the subbands box at the top of the GUI, the non-editable boxes will turn red to indicate this.

    The example below shows how to create a cube for the first and fourth subbands of a given spectrometer, reading just the channels 200 to 1200 in the first subband, and the channels 400 to 700 in the fourth:

    cubesContext = doGridding(htp=htp, subbands = Int1d([1,4]), channels=[[200,1200],[400,700]])[0]
    
  • Finally, you can also create cubes for only a selection of datasets in the HTP. You do this by specifying the index of the dataset. Here we select subbbands 2 and 4, and datasets 3, 4, and 5 from the HTP.

    cubesContext = doGridding(htp=htp, subbands=Int1d([2,4]), datasetIndices=Int1d([3,4,5]))[0]
    cube_subband_2 = cubesContext.refs["cube_WBS_V_USB_2"].product
    cube_subband_4 = cubesContext.refs["cube_WBS_V_USB_4"].product
    

    To use the GUI, create the variable datasetIndices using the format in the example above and drag it from the Variables View into the datasetIndices bullet.

15.3.5. ...to make a rotated map or use a different WCS

  • In the HTP metadata, you will find the parameter pattAngle (for pattern angle) which corresponds to the Position Angle was provided in HSpot. If the pattAngle is non-zero, then doGridding calculates a parameter called flyAngle = 180 - pattAngle (units of degrees).

    You can apply a rotation angle to the map yourself using the flyAngle parameter. For example,

    cubesContext = doGridding(htp=htp, flyAngle=50.0)[0]
    

    Note that there is also a HifiUplink parameter called flyAngle and this is not the same quantity as the flyAngle created by doGridding. It is, unfortunately, possible to find the parameter flyAngle with several values in your data. For example, in the case of the observation 1342201689, which is an OTF map carried out with a rotation angle of 65 degrees, you find the following:

    • HifiUplinkProduct: flyAngle = 65 deg

      pattAngle takes the same value as flyAngle in the HifiUplinkProduct

    • cube headers in cubesContext: flyAngle = 0 deg

      Recall that in the case a map was carried out with a non-zero rotation angle, then two sets of cubes are generated at Level 2.5 by the SPG; one set is not rotated and this is contained in cubesContext, while the rotated cubes can be found in cubesContextRotated. But if you use the doGridding task on an htp, you will only have one sets of cubes generated i.e. with the flyAngle applied. Therefore, if you wish to recreate the native orientation on the sky, you need to set flyAngle=0.0.

    • cube headers in cubesContextRotated: flyAngle = 115 deg

      doGridding calculates flyAngle = (180 - 65) degrees.

  • By default, doGridding applies a generalised least squared (GLS) projection of the input spectral coordinates into a flat map. You can also supply your own WCS to perform this projection. To create a WCS in HIPE, see Defining and using the World Coordinates System (WCS) in the Herschel Data Analysis Guide. Your WCS, for example called my_wcs, is passed to doGridding by:

    cubesContext = doGridding(htp=htp, wcs=my_wcs)[0]
    

    To use the GUI, you should pass the my_wcs variable to the wcs bullet.

15.3.6. ...with a different convolution

  • By default the convolution is performed with a Gaussian filter function for OTF maps and Nyquist or smaller sampled raster maps, while a box filter is used for greater than Nyquist sampled raster maps. A box filter is the most appropriate choice for a raster map. However, in the case of small (less than 4x4) Nyquist sampled raster maps, a box filter which has length BeamSize/PixelSize or 2.4 times the pixel size, will result in a map with pixels at all the same values - the average of all the pixels in the map. For larger Nyquist sampled maps you may like to experiment with a box filter.

    The filter type can be selected from the drop-down menu in the GUI or in the command line as:

    # box filter
    cubesContext = doGridding(htp=htp, filterType="box")[0]
    # or Gaussian filter
    cubesContext = doGridding(htp=htp, filterType="gaussian")[0]
    
  • The parameters characterising the convolution filter are automatically computed by the task, and depend on the beam size, the pixel size, and the smoothing factor. If you change any of these parameters, the optimum filter parameters are recalculated by the task, and you can see the update occuring in the GUI. Note, however, that the reverse is not true if the filter parameters are modified. In case of doubt, you can set the filter paramters to zero to force the task to auto-compute the optimal values.

    While we note that it is more intuitive to change the smoothing parameters to alter the convolution than the filter parameters, there may be cases, for example to compare a HIFI map to other data with an unusual sampling, that you would wish to do so. The filter parameters are given in pixel lengths and can be set using xFilterParams and yFilterParams. The parameters can be entered directly into the GUI, or set in the command line as demonstrated in the examples below.

    The box filter is defined by lengths (in pixels) in the x (RA) and y (Dec) directions. The Gaussian filter is defined by the length (or influence area), which is a +/- value rather than a total range, the sigma of the Gaussian function multiplied by the square root of 2, and the order of the exponent. The default value for the length of the Gaussian filter is given by 3 times the kernel beam size divided by the pixel size, i.e. 0.9 times the telescope beam size divided by the pixel size (Section 15.4.1).

    Example 1: a box filter with lengths x=0.5 pixels in the x (RA) and y=1.5 pixels in the y (Dec) directions:

    cubesContext = doGridding(htp=htp, filterType="box", xFilterParams=[0.5], yFilterParams=[1.5])[0]
    #
    # Or you can wrap the x and y parameters into one
    parameters = [Double1d([0.5]), Double1d([1.5])]
    cubesContext = doGridding(htp=htp, filterType="box", filterParams=parameters)[0]
    

    Example 2: specify the parameters length and sigma of the Gaussian filter function:

    # 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])
    yFilterParameters = Double1d([influence_area, sigma_sqrt2])
    # default case:   influence_area = 1.8; sigma_sqrt2 = 0.36
    
    cubesContext = doGridding(htp=htp,filterType="gaussian", xFilterParams= xFilterParameters, \
    yFilterParams = yFilterParameters)[0]
    
  • The doGridding task assumes an equal weighting during the convolution, as the integration time per point is equal. In the weightMode drop-down menu in the GUI this option is labelled as 'all same weight'. In the command line, this option is selected by:

    cubesContext = doGridding(htp=htp, weightMode="equal")[0]
    

    You can also choose to use weights already present in the data, which may have been set by the pipeline or by you, by choosing 'read weights column from dataset' from the drop-down weightMode menu in the GUI. To do the same in the command line, use:

    cubesContext = doGridding(htp=htp, weightMode="selection")[0]
    
  • If you wish to learn more about the details of the convolution, you can check the detail box in the GUI, or set

    cubesContext = doGridding(htp=htp, detail=True)[0]
    

    This creates an output called convolutionTable, which can be viewed with the Dataset Viewer. The metadata of the table gives information about the map, beam, and pixel sizes (but note that the x values are negative because the x-axis values increment from right to left, while the y-axis values increment from bottom to top), and also information about the filter used in the convolution and the position of the reference pixel.

    The convolutionTable itself contains the following information for each spectrum that contribute to a given pixel:

    • The pixel number, as ypixel and xpixel. The number of spectra contributing to a given pixel depends on the gridding parameters used and the spacing of the readout points on the sky.

    • The latitude and longitude of the pixel centre in decimal degrees.

    • vm, v, vp, are the offsets of the bottom, centre, and top of the y position of the pixel from the map centre, while um, u, up are the same for the x position of the pixel.

    • PointSpectrum_id gives a reference to the location of the spectrum in the HTP. For example, "row 0: box_1: container 13" refers to the first (labelled 0.0) PointSpectrum of SpectrumDataset 0014 in box 1 of the HTP. Unfortunately, that was not a typo: the ConvolutionTable labels SpectrumDatasets start from 0, while in the ContextViewer SpectrumDatasets are labelled from 1.

    • dv and du are the distance of the spectrum (in pixels) from the pixel centre in the y and x directions, while the Pythagorean distance is given by distance (in pixels) and dist_beams (in beam size).

    • If a spectrum is not included in the convolution as a consequence of a very low weighting then the blanked column will show "true".

    • du_beams and dv_beams give the distance of the spectrum in the x and y directions from the pixel centre in antenna beam sizes.

    • filter_du (filter_dv) is the value computed by the x-axis (y-axis) filter for the distance of the spectrum to the pixel centre, while filter_value is the value computed by the filter for this element, combining the filter functions along both axes.

    • spectrum_latitude and spectrum_longitude give the position of the spectrum in decimal degrees.

    • The spectrum_flux column contains the flux of the contributing spectrum.

    • The weights of each channel of the spectrum are given in the spectrum_weights column.

    • The spectra within the HTP are sorted according to their coordinates in horizontal (x) and vertical (y) directions prior to gridding. The final three columns in the ConvolutionTable refer to the indices after sorting of the first and last spectra in the current strip as well as the index of the spectrum.

15.3.7. ...to specify the map centre

The doGridding task computes the map centre as the centre of the coordinates of the input spectra. The pixel that the map centre falls in is called the reference pixel. It is possible to specify the reference pixel, in either pixel coordinates or in celestial coordinates, to doGridding. This may be useful in the case that you know the coordinates of the map centre to higher accuracy than the Herschel-HIFI Absolute Pointing Error (2").

  • When specifying the reference pixel in pixel coordinates it is important to bear in mind that the (0,0) pixel is the bottom-most, left-most pixel in the cube. In contrast, the header information of the cube follows normal FITS convention in which the bottom-most, left-most pixel is (1,1).

    For example, if you want to force the map centre to be the pixel (3.5 , 4.0), then the refPixel input will be Double1d([3.5, 4.0]). This means that in the cube header the values of CRPIX1 and CRPIX2 will be (4.5, 5.0) and the centre of the coordinates of the input spectra will be located at the pixel (3.5, 4.0).

    cubesContext = doGridding(htp=htp,refPixel = Double1d([3.5, 4.0]))[0]
    

    To use the GUI, you will need to create a variable refPixel as above and drag it to the refPixel bullet.

  • To specify the celestial coordinates of the map centre, you can enter the map position in decimal degrees directly into the refPixelCoordinates box in the GUI. Here is a command line example:

    longitude = 307.9
    latitude = 40.36
    cubesContext = doGridding(htp=htp,refPixel=Double1d([0,0]), refPixelCoordinates=Double1d([longitude, latitude]))[0]
    # or more compactly,
    refPixel = Double1d([0,0])
    refPixelCoordinates = Double1d([longitude, latitude])
    cubesContext = doGridding(htp=htp,refPixel=refPixel, refPixelCoordinates=refPixelCoordinates)[0]
    

    You can check that the cubes generated have the same crval1 (longitude) and crval2 (latitude) as the map centre you specified in the following way:

    # Value returned:  True or False
    print cubesContext.refs["cube_WBS_V_USB_2"].product.wcs.crval1 == longitude # True
    print cubesContext.refs["cube_WBS_V_USB_2"].product.wcs.crval2 == latitude  # True
    

Of course, it is possible to set the map centre such that the defined grid falls partially, or completely, outside of the observed region. In that case, doGridding will not fail but will fill the pixels in the unobserved regions with NaNs.

15.3.8. ...with selected data types

The default action of doGridding is to take the science datasets, i.e. those datasets that are on source, and with the LO setting required to observe the frequency of interest. You can also select data to be gridded according to any data type you can see in the summary table of the HTP. This is particularly important if you choose to process the reference spectra (OFF positions) found in the product Calibration -> pipeline-out to allow you, e.g., to correct for any emission contamination in the OFF positions.

For example to make a grid of the hot/cold load observations:

cubesContext = doGridding(htp=htp, datasetType="hc")[0]

To include the OFF positions (found in Calibration -> pipeline-out -> ReferenceSpectra), which are ignored by default, in the gridding with:

cubesContext = doGridding(htp=htp, ignoreOffs=False)[0]

# In the GUI, you just deselect the ignoreOffs button

In practice, the Level 2 HTP contain only science data and these options are not expected to be used by astronomers unless the doCleanUp step of the Level 2 pipeline is omitted. The option to select dataType = "other" may be of use for calibration scientists studying engineering mode observations.

15.3.9. ...to deal with NaNs

On occasion, you may see 'Not a Numbers' (NaNs) at the edge of the cube. This is not a problem with computation but a pixel that is flagged as containing spectra that have a very low weighting. Typically, this is found in maps displaying a strong zig-zag effect, and it is because the data arise from too far away from the map area. You can avoid this to some degree by using the extrapolate option:

cubesContext = doGridding(htp=htp,extrapolate=1)[0]

Furthermore, you can adjust the threshold at which pixels are flagged at, the default is 1e-5 but higher values may help remove the NaNs.

Configuration.setProperty("hcss.hifi.dp.otf.filter.threshold", "1e-1")
cubesContext = doGridding(htp=htp)[0]
cubes = cubesContext.refs["cube_WBS_V_USB_1"].product