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
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
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])
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])
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
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
pixelSize parameter but entering a non-zero
Section 15.3.5) will in order that flux is conserved when
rotating the map.
pixelSize can be modified in the command line as:
# For a square pixel cubesContext = doGridding(htp=htp, pixelSize=[16.36]) # For a rectangular pixel cubesContext = doGridding(htp=htp, pixelSize=[10.0, 20.0])
Flux Conservation in Spectral Cubes from Mapping Observations:
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 (")|
The pixel size is multiplied by a smoothing factor, which is 1 by default. You can modify the smoothing factor by editing
smoothFactor box in the GUI or in the command line by:
cubesContext = doGridding(htp=htp, smoothFactor=[2.0, 2.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
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])
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.
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)
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
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.
mkRms allows you to assess the statistics of your new combined cube compared to the original
cubes and HSpot noise predictions.
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.
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)
comoving = False in the
doGridding GUI. To do the same in the command
cubesContext = doGridding(htp=htp, comoving=False)
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
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])) # # 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:
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]])
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])) 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
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,
doGridding calculates a parameter called
flyAngle = 180 - pattAngle (units
You can apply a rotation angle to the map yourself using the
flyAngle parameter. For example,
cubesContext = doGridding(htp=htp, flyAngle=50.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:
flyAngle = 65 deg
pattAngle takes the same value as
flyAngle in the HifiUplinkProduct
cube headers in
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
by the SPG; one set is not rotated and this is contained in
cubesContext, while the rotated cubes can be
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
the sky, you need to set flyAngle=0.0.
cube headers in
flyAngle = 115 deg
flyAngle = (180 - 65) degrees.
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
cubesContext = doGridding(htp=htp, wcs=my_wcs)
To use the GUI, you should pass the
my_wcs variable to the wcs bullet.
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") # or Gaussian filter cubesContext = doGridding(htp=htp, filterType="gaussian")
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.
filter parameters are given in pixel lengths and can be set using
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]) # # 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)
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)
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")
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")
If you wish to learn more about the details of the convolution, you can check the
detail box in the GUI,
cubesContext = doGridding(htp=htp, detail=True)
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
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.
longitude of the pixel centre in decimal degrees.
vp, are the offsets of the bottom, centre, and top of the
y position of the pixel from the map centre, while
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.
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
(in beam size).
If a spectrum is not included in the convolution as a consequence of a very low weighting then the
column will show "true".
dv_beams give the distance of the spectrum in the x and y directions from
the pixel centre in antenna beam sizes.
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_longitude give the position of the spectrum in decimal degrees.
spectrum_flux column contains the flux of the contributing spectrum.
The weights of each channel of the spectrum are given in the
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.
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
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]))
To use the GUI, you will need to create a variable
refPixel as above and drag it to the
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])) # or more compactly, refPixel = Double1d([0,0]) refPixelCoordinates = Double1d([longitude, latitude]) cubesContext = doGridding(htp=htp,refPixel=refPixel, refPixelCoordinates=refPixelCoordinates)
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.
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")
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) # 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
doCleanUp step of the Level 2 pipeline is omitted. The option to select
= "other" may be of use for calibration scientists studying engineering mode observations.
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)
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) cubes = cubesContext.refs["cube_WBS_V_USB_1"].product