Source code for ctapipe.utils.table_interpolator

Table interpolation class, to allow interpolation of 2D images in any number of other
dimensions. Reads in an interim FITS table as defined below:

Tables organised as a standard FITS image file with each table as an image HDU.
First (primary) HDU must contain the following header entries:

CRPIXx, CRVALx, CRDELTAx: Number, position and pixel spacing of the reference pixel in
the image. Where x is the axis number (1 or 2).

GRIDVALS: Names of values describing the grid points
Each HDU must contain a header entry containing the grid points described by GRIDVALS
e.g. ALT =70
AZI = 0

It is also strongly recommended to include documentation for each interpolation
dimension containing at least the units. This entry ust begin with DOC, followed by the
dimension name
e.g. DOCALT = "Altitude of event (deg)"

    - Improve error handling
    - Add option for caching nd interpolated value
    - Better deal with edges of phase space
    - Allow non-linear interpolation

import numpy as np
from import fits
from scipy import interpolate

[docs]class TableInterpolator: """ This is a simple class for loading lookup tables from a fits file and interpolating between them """ def __init__(self, filename, verbose=1): """ Initialisation of class to load templates from a file and create the interpolation objects Parameters ---------- filename: string Location of Template file verbose: int Verbosity level, 0 = no logging 1 = File + interpolation point information 2 = Detailed description of interpolation points """ self.verbose = verbose if self.verbose: print("Loading lookup tables from", filename) grid, bins, template = self.parse_fits_table(filename) x_bins, y_bins = bins self.interpolator = interpolate.LinearNDInterpolator( grid, template, fill_value=0 ) self.nearest_interpolator = interpolate.NearestNDInterpolator(grid, template) self.grid_interp = interpolate.RegularGridInterpolator( (x_bins, y_bins), np.zeros([x_bins.shape[0], y_bins.shape[0]]), method="linear", bounds_error=False, fill_value=0, )
[docs] def parse_fits_table(self, filename): """ Function opens tables contained within fits files and parses them into a format recognisable by the interpolator. Parameters ---------- filename: str Name of table file Returns ------- tuple (grid points, bin centres, images) """ file = template = list() grid = list() primary_hdu = file[0].header # We require first HDU to be primary # Below definitions are standard ix, iy = primary_hdu["CRPIX2"], primary_hdu["CRPIX1"] val_x, val_y = primary_hdu["CRVAL2"], primary_hdu["CRVAL1"] print(val_x, val_y) delta_x, delta_y = primary_hdu["CRDELTA2"], primary_hdu["CRDELTA1"] nbins_x, nbins_y = primary_hdu["NAXIS2"], primary_hdu["NAXIS1"] ix *= delta_x iy *= delta_y x_bins = np.arange(val_x - ix, val_x + (delta_x * nbins_x) - ix, step=delta_x) y_bins = np.arange(val_y - iy, val_y + (delta_y * nbins_y) - iy, step=delta_y) grid_vals = primary_hdu["GRIDVALS"] points = grid_vals.split(",") if self.verbose: print("Interpolation point source be called in order", points) if self.verbose > 1: for p in points: print(p, ":", primary_hdu["DOC" + p]) for hdu in file: template.append( hdu_pt = list() for p in points: hdu_pt.append(hdu.header[p]) grid.append(np.array(hdu_pt)) print(np.array(grid)) bins = (x_bins, y_bins) return grid, bins, template
[docs] def interpolate(self, params, pixel_pos_x, pixel_pos_y): """ Parameters ---------- params: ndarray numpy array of interpolation parameters currently [energy, impact distance, xmax] pixel_pos_x: ndarray pixel position in degrees pixel_pos_y: ndarray pixel position in degrees Returns ------- ndarray of expected intensity for all pixel positions given """ image = self.interpolated_image(params) self.grid_interp.values = image points = np.array([pixel_pos_x, pixel_pos_y]) return self.grid_interp(points.T)
[docs] def interpolated_image(self, params): """ Function for creating a ful interpolated image template from the interpolation library Parameters ---------- params: ndarray numpy array of interpolation parameters currently [energy, impact distance, xmax] Returns ------- ndarray of a single image template """ image = self.interpolator(params)[0] if np.isnan(image).all(): print("Found a NaN", params) image = self.nearest_interpolator(params)[0] return image