Cygrid quick tour

The following (in-complete) code snippet demonstrates, how one would use cygrid to grid raw data, sampled at (possibly irregular) positions, (lon, lat), with values rawdata into a FITS/WCS image:

from astropy.io import fits
import matplotlib.pyplot as plt
import cygrid

# read-in data; assuming rawdata is a 2D array (i.e. a list of spectra)
lon, lat, rawdata = get_data()

# define target FITS/WCS header
header = create_fits_header()  # assume a 3D image, need only spatial part

# prepare gridder
kernelsize_sigma = 0.2

kernel_type = 'gauss1d'
kernel_params = (kernelsize_sigma, )  # must be a tuple
kernel_support = 3 * kernelsize_sigma
hpx_maxres = kernelsize_sigma / 2

mygridder = cygrid.WcsGrid(header)
mygridder.set_kernel(
    kernel_type,
    kernel_params,
    kernel_support,
    hpx_maxres,
    )

# do the actual gridding
mygridder.grid(lon, lat, rawdata)

# query result and store to disk ...
data_cube = mygridder.get_datacube()
fits.writeto('example.fits', header=header, data=data_cube)

# ... or plot
fig = plt.figure()
ax = fig.add_subplot(111, projection=WCS(header).celestial)
# plot first plane of data cube
ax.imshow(
    data_cube[0], origin='lower', interpolation='nearest'
    )
plt.show()

Here, we have omitted the code to read-in the raw data samples, lon, lat, and rawdata. The lon and lat need to be one-dimensional arrays (or lists), while in this example rawdata has two dimensions (first dimension needs to have the same length as lon and lat, second dimension is the number of spectral channels). It would also be possible that rawdata is a 1D array, in this case the output data_cube would be 2D, i.e., a map. Likewise, the code to setup the target FITS/WCS header is not explicitly shown. Anything that WCS would accept to create a spatial WCS is sufficient.

The next part in the little program is about setting up the gridder class (WCSGrid). The constructor only needs the aforementioned WCS header. But before any gridding job would work, it is necessary to define the weighting kernel parameters. Details about a proper choice of these values is contained in the user manual (Kernel parameters). In most cases you will probably use a radial-symmetric Gaussian (therefore, kernel_type = 'gauss1d') with a certain width (kernelsize_sigma). The kernel_support parameter defines out to which angular distance the kernel is computed. For Gaussians \(3\ldots4\sigma\) is sufficient. The hpx_maxres is a very technical parameter. It tells cygrid, how fine the internal HEALPix cache needs to be. We recommend to set this parameter to half the kernel size.

The gridding itself is done by calling the grid method. Depending on the size of the target FITS image and the number of raw data samples, this might take a while. Once the gridding is finished, you can obtain the gridded data array with the get_datacube method and do whatever you want with it, e.g., store it to a FITS file.

Note

This “minimal example” is also contained in a fully-working Jupyter tutorial notebook.

Gridding large data sets

If you have a huge amount of raw data to grid, which might not fit into your computer’s memory, cygrid allows you to process the data in smaller chunks:

mygridder = cygrid.WcsGrid(header)
mygridder.set_kernel(...)

for chunk in chunks:
    lon, lat, rawdata = get_data(chunk)
    mygridder.grid(lon, lat, rawdata)

data_cube = mygridder.get_datacube()

More information can be found in the user manual (Decrease memory footprint).

Sight-line gridding

Sometimes, it can be useful to re-sample the input raw data set for a given list of coordinate pairs, e.g., if you want to grid to a funky coordinate system or projection, which is not supported by FITS/WCS. For this, one can use the SlGrid class:

mygridder = cygrid.SlGrid(target_lon, target_lat)
mygridder.set_kernel(...)

mygridder.grid(input_lon, input_lat, input_signal)

target_signal = gridder.get_datacube()

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(target_lon, target_lat, c=target_signal)
plt.show()