utils - Utilities¶
Introduction¶
gammapy.utils is a collection of utility functions that are used in many
places or don’t fit in one of the other packages.
Since the various sub-modules of gammapy.utils are mostly unrelated, they
are not imported into the top-level namespace. Here are some examples of how to
import functionality from the gammapy.utils sub-modules:
from gammapy.utils.random import sample_sphere
sample_sphere(size=10)
from gammapy.utils import random
random.sample_sphere(size=10)
Time handling in Gammapy¶
See gammapy.utils.time.
Time format and scale¶
In Gammapy, astropy.time.Time objects are used to represent times:
>>> from astropy.time import Time
>>> Time(['1999-01-01T00:00:00.123456789', '2010-01-01T00:00:00'])
<Time object: scale='utc' format='isot' value=['1999-01-01T00:00:00.123' '2010-01-01T00:00:00.000']>
Note that Astropy chose format='isot' and scale='utc' as default and in
Gammapy these are also the recommended format and time scale.
Warning
Actually what’s written here is not true. In CTA it hasn’t been decided if
times will be in utc or tt (terrestial time) format.
Here’s a reminder that this needs to be settled / updated: https://github.com/gammapy/gammapy/issues/284
When other time formats are needed it’s easy to convert, see the time format section and table in the Astropy docs.
E.g. sometimes in astronomy the modified Julian date mjd is used and for
passing times to matplotlib for plotting the plot_date format should be
used:
>>> from astropy.time import Time
>>> time = Time(['1999-01-01T00:00:00.123456789', '2010-01-01T00:00:00'])
>>> time.mjd
array([ 51179.00000143, 55197. ])
>>> time.plot_date
array([ 729755.00000143, 733773. ])
Converting to other time scales is also easy, see the time scale section, diagram and table in the Astropy docs.
E.g. when converting celestial (RA/DEC) to horizontal (ALT/AZ) coordinates, the
sidereal time is needed. This
is done automatically by astropy.coordinates.AltAz when the
astropy.coordinates.AltAz.obstime is set with a Time object in
any scale, no need for explicit time scale transformations in Gammapy (although
if you do want to explicitly compute it, it’s easy, see here).
The Fermi-LAT time systems in a nutshell page gives a good, brief explanation
of the differences between the relevant time scales UT1, UTC and TT.
Mission elapsed times (MET)¶
[MET] time references are times representing UTC seconds after a specific origin. Each experiment may have a different MET origin that should be included in the header of the corresponding data files. For more details see Fermi-LAT time systems in a nutshell.
It’s not clear yet how to best implement METs in Gammapy, it’s one of the tasks here: https://github.com/gammapy/gammapy/issues/284
For now, we use the gammapy.utils.time.time_ref_from_dict,
gammapy.utils.time.time_relative_to_ref and gammapy.utils.time.absolute_time functions
to convert MET floats to Time objects via the reference times
stored in FITS headers.
Energy handling in Gammapy¶
Basics¶
Most objects in Astronomy require an energy axis, e.g. counts spectra or effective area tables. In general, this axis can be defined in two ways.
- As an array of energy values. E.g. the Fermi-LAT diffuse flux is given at certain energies and those are stored in an ENERGY FITS table extension. In Gammalib this is represented by GEnergy.
- As an array of energy bin edges. This is usually stored in EBOUNDS tables, e.g. for Fermi-LAT counts cubes. In Gammalib this is represented by GEbounds.
In Gammapy both the use cases are handled by two seperate classes:
gammapy.utils.energy.Energy for energy values and
gammapy.utils.energy.EnergyBounds for energy bin edges
Energy¶
The Energy class is a subclass of Quantity and thus has the
same functionality plus some convenience functions for fits I/O
>>> from gammapy.utils.energy import Energy
>>> energy = Energy([1,2,3], 'TeV')
>>> hdu = energy.to_fits()
>>> type(hdu)
<class 'astropy.io.fits.hdu.table.BinTableHDU'>
EnergyBounds¶
The EnergyBounds class is a subclass of Energy. Additional functions are available e.g. to compute the bin centers
>>> from gammapy.utils.energy import EnergyBounds
>>> ebounds = EnergyBounds.equal_log_spacing(1, 10, 8, 'GeV')
>>> ebounds.size
9
>>> ebounds.nbins
8
>>> center = ebounds.log_centers
>>> center
<Energy [ 1.15478198, 1.53992653, 2.05352503, 2.73841963, 3.65174127,
4.86967525, 6.49381632, 8.65964323] GeV>
Reference/API¶
gammapy.utils.energy Module¶
Classes¶
Energy |
Energy quantity scalar or array. |
EnergyBounds |
EnergyBounds array. |
gammapy.utils.units Module¶
Units and Quantity related helper functions
Functions¶
standardise_unit(unit) |
Standardise unit. |
unit_from_fits_image_hdu(header) |
Read unit from a FITS image HDU. |
gammapy.utils.coordinates Package¶
Astronomical coordinate calculation utility functions.
Functions¶
angle_to_radius(angle, distance) |
Radius (pc), distance(kpc), angle(deg) |
cartesian(r, theta) |
Convert polar coordinates to cartesian coordinates. |
flux_to_luminosity(flux, distance) |
Distance is assumed to be in kpc |
fov_to_sky(lon, lat, lon_pnt, lat_pnt) |
Transform field-of-view coordinates to sky coordinates. |
galactic(x, y, z[, obs_pos]) |
Compute galactic coordinates lon, lat (deg) and distance (kpc) for given position in cartesian coordinates (kpc). |
luminosity_to_flux(luminosity, distance) |
Distance is assumed to be in kpc |
minimum_separation(lon1, lat1, lon2, lat2) |
Compute minimum distance of each (lon1, lat1) to any (lon2, lat2). |
motion_since_birth(v, age, theta, phi) |
Compute motion of a object with given velocity, direction and age. |
pair_correlation(lon, lat, theta_bins) |
Compute pair correlation function for points on the sphere. |
polar(x, y) |
Convert cartesian coordinates to polar coordinates. |
radius_to_angle(radius, distance) |
Radius (pc), distance(kpc), angle(deg) |
sky_to_fov(lon, lat, lon_pnt, lat_pnt) |
Transform sky coordinates to field-of-view coordinates. |
velocity_glon_glat(x, y, z, vx, vy, vz) |
Compute projected angular velocity in galactic coordinates. |
Variables¶
D_SUN_TO_GALACTIC_CENTER |
A Quantity represents a number with some associated unit. |
gammapy.utils.interpolation Module¶
Interpolation utilities
Functions¶
interpolation_scale([scale]) |
Interpolation scaling. |
Classes¶
ScaledRegularGridInterpolator(points, values) |
Thin wrapper around scipy.interpolate.RegularGridInterpolator. |
gammapy.utils.table Module¶
Table helper utilities.
Functions¶
table_standardise_units_copy(table) |
Standardise units for all columns in a table in a copy. |
table_standardise_units_inplace(table) |
Standardise units for all columns in a table in place. |
table_row_to_dict(row[, make_quantity]) |
Make one source data dict. |
table_from_row_data(rows, **kwargs) |
Helper function to create table objects from row data. |
gammapy.utils.fits Module¶
Gammapy FITS utilities¶
FITS tables¶
In Gammapy we use the nice astropy.table.Table class a lot to represent all
kinds of data (e.g. event lists, spectral points, light curves, source catalogs).
The most common format to store tables is FITS. In this section we show examples
and mention some limitations of Table FITS I/O.
Also, note that if you have the choice, you might want to use a better format
than FITS to store tables. All of these are nice and have very good support
in Astropy: ECSV, HDF5, ASDF.
In Astropy, there is the Table class with a nice data model
and API. Let’s make an example table object that has some metadata on the
table and columns of different types:
>>> from astropy.table import Table, Column
>>> table = Table(meta={'version': 42})
>>> table['a'] = [1, 2]
>>> table['b'] = Column([1, 2], unit='m', description='Velocity')
>>> table['c'] = ['x', 'yy']
>>> table
<Table length=2>
a b c
m
int64 int64 str2
----- ----- ----
1 1 x
2 2 yy
>>> table.info()
<Table length=2>
name dtype unit description
---- ----- ---- -----------
a int64
b int64 m Velocity
c str2
Writing and reading the table to FITS is easy:
>>> table.write('table.fits')
>>> table2 = Table.read('table.fits')
and works very nicely, column units and description round-trip:
>>> table2
<Table length=2>
a b c
m
int64 float64 bytes2
----- ------- ------
1 1.0 x
2 2.0 yy
>>> table2.info()
<Table length=2>
name dtype unit description
---- ------- ---- -----------
a int64
b float64 m Velocity
c bytes2
This is with Astropy 3.0. In older versions of Astropy this didn’t use to work, namely column description was lost.
Looking at the FITS header and table2.meta, one can see that
they are cheating a bit, storing table meta in COMMENT:
>>> fits.open('table.fits')[1].header
XTENSION= 'BINTABLE' / binary table extension
BITPIX = 8 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 18 / length of dimension 1
NAXIS2 = 2 / length of dimension 2
PCOUNT = 0 / number of group parameters
GCOUNT = 1 / number of groups
TFIELDS = 3 / number of table fields
TTYPE1 = 'a '
TFORM1 = 'K '
TTYPE2 = 'b '
TFORM2 = 'K '
TUNIT2 = 'm '
TTYPE3 = 'c '
TFORM3 = '2A '
VERSION = 42
COMMENT --BEGIN-ASTROPY-SERIALIZED-COLUMNS--
COMMENT datatype:
COMMENT - {name: a, datatype: int64}
COMMENT - {name: b, unit: m, datatype: int64, description: Velocity}
COMMENT - {name: c, datatype: string}
COMMENT meta:
COMMENT __serialized_columns__: {}
COMMENT --END-ASTROPY-SERIALIZED-COLUMNS--
>>> table2.meta
OrderedDict([('VERSION', 42),
('comments',
['--BEGIN-ASTROPY-SERIALIZED-COLUMNS--',
'datatype:',
'- {name: a, datatype: int64}',
'- {name: b, unit: m, datatype: int64, description: Velocity}',
'- {name: c, datatype: string}',
'meta:',
' __serialized_columns__: {}',
'--END-ASTROPY-SERIALIZED-COLUMNS--'])])
TODO: we’ll have to see how to handle this, i.e. if we want that behaviour or not, and how to get consistent output accross Astropy versions. See https://github.com/astropy/astropy/issues/7364
Let’s make sure for the following examples we have a clean table.meta
like we did at the start:
>>> table.meta.pop('comments', None)
If you want to avoid writing to disk, the way to directly convert between
Table and BinTableHDU is like this:
>>> hdu = fits.BinTableHDU(table)
This calls astropy.io.fits.table_to_hdu in BinTableHDU.__init__,
i.e. if you don’t pass extra options, this is equivalent to
>>> hdu = fits.table_to_hdu(table)
However, in this case, the column metadata that is serialised is
doesn’t include the column description.
TODO: how to get consistent behaviour and FITS headers?
>>> hdu.header
XTENSION= 'BINTABLE' / binary table extension
BITPIX = 8 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 18 / length of dimension 1
NAXIS2 = 2 / length of dimension 2
PCOUNT = 0 / number of group parameters
GCOUNT = 1 / number of groups
TFIELDS = 3 / number of table fields
VERSION = 42
TTYPE1 = 'a '
TFORM1 = 'K '
TTYPE2 = 'b '
TFORM2 = 'K '
TUNIT2 = 'm '
TTYPE3 = 'c '
TFORM3 = '2A '
Somewhat surprisingly, Table(hdu) doesn’t work and there is no
hdu_to_table function; instead you have to call Table.read
if you want to convert in the other direction:
>>> table2 = Table.read(hdu)
>>> table2.info()
<Table length=2>
name dtype unit
---- ----- ----
a int64
b int64 m
c str2
Another trick worth knowing about is how to read and write multiple tables
to one FITS file. There is support in the Table API to read any HDU
from a FITS file with multiple HDUs via the hdu option to Table.read;
you can pass an integer HDU index or an HDU extension name string
(see FITS).
For writing (or if you prefer also for reading) multiple tables, you should
use the in-memory conversion to HDU objects and the HDUList
like this:
hdu_list = fits.HDUList([
fits.PrimaryHDU(),
fits.BinTableHDU(table, name='spam'),
fits.BinTableHDU(table, name='ham'),
])
hdu_list.info()
hdu_list.writeto('tables.fits')
For further information on Astropy, see the Astropy docs at Data Tables (astropy.table) and FITS.
We will have to see if / what we need here in gammapy.utils.fits
as a stable and nice interface on top of what Astropy provides.
Functions¶
energy_axis_to_ebounds(energy) |
Convert EnergyBounds to OGIP EBOUNDS extension. |
earth_location_from_dict(meta) |
Create EarthLocation from FITS header dict. |
Classes¶
SmartHDUList(hdu_list) |
A FITS HDU list wrapper with some sugar. |
gammapy.utils.random Module¶
Random sampling for some common distributions
Functions¶
get_random_state(init) |
Get a numpy.random.RandomState instance. |
sample_sphere(size[, lon_range, lat_range, …]) |
Sample random points on the sphere. |
sample_sphere_distance([distance_min, …]) |
Sample random distances if the 3-dim space density is constant. |
sample_powerlaw(x_min, x_max, gamma[, size, …]) |
Sample random values from a power law distribution. |
gammapy.utils.distributions Package¶
Utility functions / classes for working with distributions (e.g. probability density functions)
Functions¶
density(func) |
Radial surface density of a given one dimensional PDF. |
draw(low, high, size, dist[, random_state]) |
Allows drawing of random numbers from any distribution. |
normalize(func, x_min, x_max) |
Normalize a 1D function over a given range. |
pdf(func) |
One-dimensional PDF of a given radial surface density. |
Classes¶
GeneralRandom(pdf, min_range, max_range[, …]) |
Fast random number generation with an arbitrary pdf of a continuous variable x. |
GeneralRandomArray(pdf) |
Draw random indices from a discrete probability distribution given by a numpy array. |
gammapy.utils.scripts Module¶
Utils to create scripts and command-line tools
Functions¶
read_yaml(filename[, logger]) |
Read YAML file |
write_yaml(dictionary, filename[, logger]) |
Write YAML file. |
make_path(path) |
Expand environment variables on Path construction. |
recursive_merge_dicts(a, b) |
Recursively merge two dictionaries. |
gammapy.utils.testing Module¶
Utilities for testing
Functions¶
requires_dependency(name) |
Decorator to declare required dependencies for tests. |
requires_data(name) |
Decorator to declare required data for tests. |
assert_quantity_allclose(actual, desired[, …]) |
Assert all-close for astropy.units.Quantity objects. |
assert_wcs_allclose(wcs1, wcs2) |
Assert all-close for astropy.wcs.WCS objects. |
assert_skycoord_allclose(actual, desired) |
Assert all-close for astropy.coordinates.SkyCoord objects. |
assert_time_allclose(actual, desired[, atol]) |
Assert all-close for astropy.time.Time objects. |
gammapy.utils.wcs Module¶
WCS related utility functions.
Functions¶
linear_wcs_to_arrays(wcs, nbins_x, nbins_y) |
Make a 2D linear binning from a WCS object. |
linear_arrays_to_wcs(name_x, name_y, …) |
Make a 2D linear WCS object from arrays of bin edges. |
get_wcs_ctype(wcs) |
Get celestial coordinate type of WCS instance. |
get_resampled_wcs(wcs, factor, downsampled) |
Get resampled WCS object. |
gammapy.utils.nddata Module¶
Utility functions and classes for n-dimensional data and axes.
Functions¶
sqrt_space(start, stop, num) |
Return numbers spaced evenly on a square root scale. |
Classes¶
NDDataArray(axes[, data, meta, interp_kwargs]) |
ND Data Array Base class |
DataAxis(nodes[, name, interpolation_mode]) |
Data axis to be used with NDDataArray |
BinnedDataAxis(lo, hi, **kwargs) |
Data axis for binned data |
gammapy.utils.time Module¶
Time related utility functions.
Functions¶
time_ref_from_dict(meta[, format, scale]) |
Calculate the time reference from metadata. |
time_ref_to_dict(time[, scale]) |
TODO: document and test. |
time_relative_to_ref(time, meta) |
Convert a time using an existing reference. |
absolute_time(time_delta, meta) |
Convert a MET into human readable date and time. |
gammapy.utils.fitting Package¶
Fitting framework and backends.
Functions¶
confidence_iminuit(minuit, parameters, …) |
|
covariance_iminuit(minuit) |
|
covariance_scipy(parameters, function) |
|
covariance_sherpa() |
|
mncontour(minuit, parameters, x, y, …) |
|
optimize_iminuit(parameters, function, **kwargs) |
iminuit optimization |
optimize_scipy(parameters, function, **kwargs) |
|
optimize_sherpa(parameters, function, **kwargs) |
Sherpa optimization wrapper method. |
Classes¶
Dataset |
Dataset abstract base class. |
Datasets(datasets[, mask]) |
Join multiple datasets. |
Fit(datasets) |
Fit class. |
Likelihood(function, parameters) |
Wrapper of the likelihood function used by the optimiser. |
Model([params]) |
Model base class. |
Parameter(name, factor[, unit, scale, min, …]) |
A model parameter. |
Parameters([parameters, covariance, …]) |
List of Parameter. |