Mapping observations from OD 835 onwards and processed with HIPE 8 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 or you should run the script given in the introduction to this chapter, where more information is given on this topic, Section 12.1.
Despite the reduced need to regrid 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 12.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:
cubes = 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. cubes = griddingTask(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. For earlier maps, different estimates are used depending on what information can be found in the metadata. 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 now does not affect the pixelSize parameter but entering a non-zero flyAngle (see Section 12.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 cubes=doGridding(htp=htp, pixelSize=[16.36]) # For a rectangular pixel cubes = doGridding(htp=htp, pixelSize=[10.0, 20.0])
For convenience, the beam widths per HIFI band and beam spacing for half-beam and Nyquist sampled maps are given in the table below.
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:
cubes=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 12.3.6 below.
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, by default these will appear as zeros to ensure the autocomputation. In the command line that can be done as:
cubes=doGridding(htp=htp, mapSize=[8.0, 7.0])
Note that if the map dimensions 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.
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. In the script below we combine Level 2 H and V polarisation data to create an HTP containing both before using the doGridding task to create a combined polarisation map. This method can also be used to create combined maps of data from different observations at the same frequency.
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.
# Get the level 2 WBS data from two polarizations of our observation htpv=obs.refs["level2"].product.refs["WBS-V-USB"].product htph=obs.refs["level2"].product.refs["WBS-H-USB"].product # # Define the reference frequency grid from one of the HTPs, and crop both to a # common range since they may be on a slightly different scale with respect to # one another (they are frequency calibrated from their own comb measurements). newGrid = mkFreqGrid(htph) # htph=doFreqGrid(htp=htph, grid=newGrid) htpv=doFreqGrid(htp=htpv, grid=newGrid) # # The next step now combines both HTPs index=htpv.getLastIndex() for i in range (1,htpv.count+1): htpv.set(htph[i] ,index+i) # # Update the new HTP's summary info to take the new frequency grid into # account, otherwise you just get the same HTP as before from herschel.hifi.pipeline.util.data import HifiSpectrumDataSupport HifiSpectrumDataSupport.adjustSummaryTable(htpv); # cubes = doGridding(htp=htpv) cubesContext = doGridding.cubesContext #
Cubes produced by the automatic pipeline are created assuming a fixed target. Therefore, you need to regrid the data such that the map centre follows the moving target. This is done by checking the comoving box in the doGridding GUI or in the command line as:
cubes = doGridding(htp=htp, comoving=True)
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:
cubes= doGridding(htp=htp, subbands=Int1d([1, 4]) # # Extract the cube for subband 1 cube1 = cubes # Extract the cube for subband 4 cube 4 = cubes #
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 SpectrumExplorer 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:
cubes = doGridding(htp=htp, subband = 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.
cubes = doGridding(htp=htp, subbands=Int1d([2,4]), datasetIndices=([3,4,5])) cube_subband_2 = cubes cube_subband_4 = cubes
To use the GUI, create the variable
datasetIndices using the format in the example above and
drag it from the Variables View into the datssetIndices bullet.
By default, doGridding applies a generalized 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 chapter 4.3.1 "Working with the World Coordinate System" in the Herschel Data Analysis Guide. Your WCS, for example called my_wcs, is passed to doGridding by:
cubes = doGridding(htp=htp, wcs=my_wcs)
To use the GUI you should pass the my_wcs variable to the wcs bullet.
The doGridding task assumes that maps are not rotated. The spatial parameters of the input WCS are copied to the output cube and this can be used to rotate the cube. In the HTP metadata, you will find the parameter pattAngle (for pattern angle) which corresponds to the Position Angle you have provided in HSpot. If the pattAngle is not equal to 0 then you can set the doGridding parameter flyAngle to (180deg - pattAngle), e.g. for pattAngle=15, flyAngle is then set to 165.
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 cubes = doGridding(htp=htp, filterType="box") # or Gaussian filter cubes = 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 recalulated 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 intuititive 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 multipled 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 beam size divided by the pixel size.
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:
cubes = 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])] cubes = 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]) # 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)
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:
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:
If you wish to learn more about the details of the convolution, you can check the detail box in the GUI, or set
cubes = doGridding(htp=htp, detail=True)
This creates an output called convolutionTable, which can be viewed with the Dataset Viewer. The meta data 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
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,
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
dist_beams (in beam size).
If a spectrum is not included in the convolution as a consequence of a very low weighting then the
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.
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).
cubes = 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. A command line example:
longitude = 307.9 latitude = 40.36 cubes = doGridding(htp=htp,refPixel=Double1d([0,0]), refPixelCoordinates=Double1d([longitude, latitude])) # or more compactly, refPixel = Double1d([0,0]) refPixelCoordinates = Double1d([longitude, latitude]) cubes = 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:
# 1 is returned if True, 0 if False print cubes.wcs.crval1 == longitude # 1, True print cubes.wcs.crval2 == latitude # 1, 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, that is 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. For example to make a grid of the hot/cold load observations:
cubes = doGridding(htp=htp, datasetType="hc")
You can also choose to include the OFF positions, which are ignored by default, in the gridding with:
cubes = doGridding(htp=htp, datasetType="science", ignoreOffs=false)
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.
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:
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") cubes=doGridding(htp=htp) cub=cubes