Analysis Tools: Fit models to data

The tool module contains the analysis tools in the PDR Toolbox. All tools are derived from ToolBase.

For examples how to use LineRatioFit, see the notebooks PDRT_Example_Find_n_G0_Single_Pixel.ipynb and PDRT_Example_Make_n_G0_maps.ipynb.

For an example how to use H2ExcitationFit and ExcitationPlot see the notebook PDRT_Example_H2_Excitation.ipynb.

ToolBase

The base class of all tools. Tools have a run() method both of which subclasses must define.

class pdrtpy.tool.toolbase.ToolBase[source]

Bases: object

Base class object for PDR Toolbox tools. This class implements a simple interface with a run method. Tools will generally do some set up such as reading in observational data before run() can be invoked.

Attributes:
has_maps

Are the Measurements used map-based?.

has_scalar

Are the Measurements used scalars.

has_vectors

Are the Measurements used a Nx1 vector, e.g.

Methods

run()

Runs the tool.
property has_maps

Are the Measurements used map-based?. (i.e., have 2 spatial axes)

Returns:

True, if the observational inputs are spatial maps, False otherwise

Return type:

bool

property has_scalar

Are the Measurements used scalars.

Returns:

True, if the observational inputs are scalars, False otherwise

Return type:

bool

property has_vectors

Are the Measurements used a Nx1 vector, e.g. read in from a table with from_table().

Returns:

True, if the observational inputs are a vector, False otherwise

Return type:

bool

run()[source]

Runs the tool. Each subclass Tool must implement its own run() method.

Excitation Diagram Fitting

H2ExcitationFit is a tool for fitting temperature, column density, and ortho-to-para ratio in \(H_2\) excitation diagrams. A two temperature model is assumed, and the fit will find \(T_{hot}, T_{cold}, N_{hot}(H_2), N_{cold}(H_2),\) and optionally OPR. The base class ExcitationFit can be used to create a tool to fit a different molecule.

class pdrtpy.tool.h2excitation.ExcitationFit(measurements=None, constantsfile=None)[source]

Bases: ToolBase

Base class for creating excitation fitting tools for various species.

Parameters:

measurements (list of Measurement.) – Input measurements to be fit.

Attributes:
has_maps

Are the Measurements used map-based?.

has_scalar

Are the Measurements used scalars.

has_vectors

Are the Measurements used a Nx1 vector, e.g.

Methods

add_measurement(m)

Add an intensity Measurement to internal dictionary used to compute the excitation diagram.

remove_measurement(identifier)

Delete a measurement from the internal dictionary used to compute column densities.

replace_measurement(m)

Safely replace an existing intensity Measurement.

run()

Runs the tool.
add_measurement(m)[source]

Add an intensity Measurement to internal dictionary used to compute the excitation diagram. This method can also be used to safely replace an existing intensity Measurement.

Parameters:

m – A Measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)

remove_measurement(identifier)[source]

Delete a measurement from the internal dictionary used to compute column densities. Any associated column density will also be removed.

Parameters:

identifier (str) – the measurement identifier

Raises:

KeyError – if identifier not in existing Measurements

replace_measurement(m)[source]

Safely replace an existing intensity Measurement. Do not change a Measurement in place, use this method. Otherwise, the column densities will be inconsistent.

Parameters:

m – A Measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)

class pdrtpy.tool.h2excitation.H2ExcitationFit(measurements=None, constantsfile='RoueffEtAl.tab')[source]

Bases: ExcitationFit

Tool for fitting temperatures, column densities, and ortho-to-para ratio(OPR) from an \(H_2\) excitation diagram. It takes as input a set of \(H_2\) rovibrational line observations with errors represented as Measurement.

Often, excitation diagrams show evidence of both “hot” and “cold” gas components, where the cold gas dominates the intensity in the low J transitions and the hot gas dominates in the high J transitions. Given data over several transitions, one can fit for \(T_{cold}, T_{hot}, N_{total} = N_{cold}+ N_{hot}\), and optionally OPR. One needs at least 5 points to fit the temperatures and column densities (slope and intercept \(\times 2\)), though one could compute (not fit) them with only 4 points. To additionally fit OPR, one should have 6 points (5 degrees of freedom).

Once the fit is done, ExcitationPlot can be used to view the results.

Parameters:

measurements (list of Measurement.) – Input \(H_2\) measurements to be fit.

Attributes:
cold_colden

The fitted cold gas total column density

fit_result

The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables.

has_maps

Are the Measurements used map-based?.

has_scalar

Are the Measurements used scalars.

has_vectors

Are the Measurements used a Nx1 vector, e.g.

hot_colden

The fitted hot gas total column density

intensities

The stored intensities.

numcomponents

Number of temperature components in the fit

opr

The ortho-to-para ratio (OPR)

opr_fitted

Was the ortho-to-para ratio fitted?

tcold

The fitted cold gas excitation temperature

temperature

The fitted gas temperatures, returned in a dictionary with keys ‘hot’ and ‘cold’.

thot

The fitted hot gas excitation temperature

total_colden

The fitted total column density

Methods

add_measurement(m)

Add an intensity Measurement to internal dictionary used to compute the excitation diagram.

average_column_density([position, size, ...])

Compute the average column density over a spatial box.

colden(component)

The column density of hot or cold gas component, or total column density.

column_densities([norm, unit, line])

The computed upper state column densities of stored intensities

energies([line])

Upper state energies of stored intensities, in K.

gu(id, opr)

Get the upper state statistical weight $g_u$ for the given transition identifer, and, if the transition is odd-$J$, scale the result by the given ortho-to-para ratio.

intensity(colden)

Given an upper state column density \(N_u\), compute the intensity \(I\).

remove_measurement(identifier)

Delete a measurement from the internal dictionary used to compute column densities.

replace_measurement(m)

Safely replace an existing intensity Measurement.

run([position, size, fit_opr])

Fit the \(log N_u-E\) diagram with two excitation temperatures, a hot \(T_{ex}\) and a cold \(T_{ex}\).

upper_colden(intensity, unit)

Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
average_column_density(position=None, size=None, norm=True, unit=Unit("1 / cm2"), line=True, clip=<Quantity -1.e+40 1 / cm2>)[source]

Compute the average column density over a spatial box. The box is created using astropy.nddata.utils.Cutout2D.

Parameters:

  • position (tuple) – The position of the cutout array’s center with respect to the data array. The position can be specified either as a (x, y) tuple of pixel coordinates.

  • size (int, array_like`) – The size of the cutout array along each axis. If size is a scalar number or a scalar Quantity, then a square cutout of size will be created. If size has two elements, they should be in (nx,ny) order [this is the opposite of Cutout2D signature]. Scalar numbers in size are assumed to be in units of pixels. Default value of None means use all pixels (position is ignored)

  • norm (bool) – if True, normalize the column densities by the statistical weight of the upper state, \(g_u\). For ortho-$H_2$ $g_u = OPR times (2J+1)$, for para-$H_2$ $g_u=2J+1$. In LTE, $OPR = 3$.

  • unit (str or astropy.units.Unit) – The units in which to return the column density. Default: \({\rm cm}^{-2}\)

  • line (bool) – if True, the returned dictionary index is the Line name, otherwise it is the upper state \(J\) number.

  • clip (astropy.units.Quantity) – Column density value at which to clip pixels. Pixels with column densities below this value will not be used in the average. Default: a large negative number, which translates to no clipping.

Returns:

dictionary of column density Measurements, with keys as \(J\) number or Line name

Return type:

dict

property cold_colden

The fitted cold gas total column density

Return type:

Measurement

colden(component)[source]

The column density of hot or cold gas component, or total column density.

Parameters:

component (str) – ‘hot’, ‘cold’, or ‘total

Return type:

Measurement

column_densities(norm=False, unit=Unit('1 / cm2'), line=True)[source]

The computed upper state column densities of stored intensities

Parameters:

  • norm (bool) – if True, normalize the column densities by the statistical weight of the upper state, \(g_u\). Default: False

  • unit (str or astropy.units.Unit) – The units in which to return the column density. Default: \({\\rm }cm^{-2}\)

  • line (bool) – if True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: False

Returns:

dictionary of column densities indexed by upper state \(J\) number or Line name. Default: False means return indexed by \(J\).

Return type:

dict

energies(line=True)[source]

Upper state energies of stored intensities, in K.

Parameters:

line (bool) – if True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: False

Returns:

dictionary indexed by upper state \(J\) number or Line name. Default: False means return indexed by \(J\).

Return type:

dict

property fit_result

The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables.

Return type:

lmfit.model.ModelResult

gu(id, opr)[source]

Get the upper state statistical weight $g_u$ for the given transition identifer, and, if the transition is odd-$J$, scale the result by the given ortho-to-para ratio. If the transition is even-$J$, the LTE value is returned.

Parameters:

  • id (str) – the measurement identifier

  • opr (float) –

Raises:

KeyError – if id not in existing Measurements

Return type:

float

property hot_colden

The fitted hot gas total column density

Return type:

Measurement

property intensities

The stored intensities. See add_measurement()

Return type:

list of Measurement

intensity(colden)[source]

Given an upper state column density \(N_u\), compute the intensity \(I\).

\[I = {A \Delta E~N_u \over 4\pi}\]

where \(A\) is the Einstein A coefficient and \(\Delta E\) is the energy of the transition.

Parameters:

colden (Measurement) – upper state column density

Returns:

optically thin intensity

Return type:

Measurement

property numcomponents

Number of temperature components in the fit

Return type:

int

property opr

The ortho-to-para ratio (OPR)

Returns:

The fitted OPR is it was determined in the fit, otherwise the canonical LTE OPR

Return type:

Measurement

property opr_fitted

Was the ortho-to-para ratio fitted?

Returns:

True if OPR was fitted, False if canonical LTE value was used

Return type:

bool

run(position=None, size=None, fit_opr=False, **kwargs)[source]

Fit the \(log N_u-E\) diagram with two excitation temperatures, a hot \(T_{ex}\) and a cold \(T_{ex}\).

If position and size are given, the data will be averaged over a spatial box before fitting. The box is created using astropy.nddata.utils.Cutout2D. If position or size is None, the data are averaged over all pixels. If the Measurements are single values, these arguments are ignored.

Parameters:

  • position (tuple) – The position of the cutout array’s center with respect to the data array. The position can be specified either as a (x, y) tuple of pixel coordinates.

  • size (int, array_like`) – The size of the cutout array along each axis in pixels. If size is a scalar number or a scalar Quantity, then a square cutout of size will be created. If size has two elements, they should be in (nx, ny) order [this is the opposite of Cutout2D signature]. Scalar numbers in size are assumed to be in units of pixels. Default value of None means use all pixels (position is ignored)

  • fit_opr (bool) – Whether to fit the ortho-to-para ratio or not. If True, the OPR will be varied to determine the best value. If False, the OPR is fixed at the canonical LTE value of 3.

property tcold

The fitted cold gas excitation temperature

Return type:

Measurement

property temperature

The fitted gas temperatures, returned in a dictionary with keys ‘hot’ and ‘cold’. :rtype: dict

property thot

The fitted hot gas excitation temperature

Return type:

Measurement

property total_colden

The fitted total column density

Return type:

Measurement

upper_colden(intensity, unit)[source]

Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission. Units of \(I\) need to be equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\).

\[ \begin{align}\begin{aligned}I &= {A \Delta E~N_u \over 4\pi}\\N_u &= 4\pi {I\over A\Delta E}\end{aligned}\end{align} \]

where \(A\) is the Einstein A coefficient and \(\Delta E\) is the energy of the transition.

Parameters:

  • intensity (Measurement) – A Measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)

  • unit (str or astropy.units.Unit) – The units in which to return the column density. Default: \({\rm }cm^{-2}\)

Returns:

a Measurement of the column density.

Return type:

Measurement

LineRatioFit

LineRatioFit is a tool for determining photodissociation region external radiation field and hydrogen nucleus density (commonly given as \(G_0\) and \(n\)) from measured spectral line intensity ratios.

class pdrtpy.tool.lineratiofit.LineRatioFit(modelset=<pdrtpy.modelset.ModelSet object>, measurements=None)[source]

Bases: ToolBase

LineRatioFit is a tool to fit observations of intensity ratios to a set of PDR models. It takes as input a set of observations with errors represented as Measurement and ModelSet for the models to which the data will be fitted. The observations should be spectral line or continuum intensities. They can be spatial maps or single pixel values. They should have the same spatial resolution.

The models to be fit are stored as intensity ratios. The input observations will be use to create ratios that correspond to models. From there a minimization fit is done to determine the density and radiation field that best fit the data. At least 3 observations are needed in order to make at least 2 ratios. With fewer ratios, no fitting can be done. More ratios generally means better determined density and radiation field, assuming the data are consistent with each other.

Once the fit is done, LineRatioPlot can be used to view the results.

Parameters:

  • modelset (ModelSet) – The set of PDR models to use for fitting.

  • measurements (list or dict of Measurement. If dict, the keys should be the Measurement identifiers.) – Input measurements to be fit.

Attributes:
density

The computed hydrogen nucleus density value(s).

fit_result

The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables.

has_maps

Are the Measurements used map-based?.

has_scalar

Are the Measurements used scalars.

has_vectors

Are the Measurements used a Nx1 vector, e.g.

measurementIDs

The stored measurement IDs, which are strings.

measurements

The stored measurements as a dictionary with Measurement IDs as keys

modelset

The underlying ModelSet

observed_ratios

The list of the observed line ratios that have been input so far.

radiation_field

The computed radiation field value(s).

ratiocount

The number of ratios that match models available in the current ModelSet given the current set of measurements

table

Construct the table of input Measurements, and if the fit has been run, the density, radiation field, and \(\chi^2\) values

Methods

add_measurement(m)

Add a Measurement to internal dictionary used to compute ratios.

chisq([min])

The computed chisquare value(s).

read_models([unit])

Given a list of measurement IDs, find and open the FITS files that have matching ratios and populate the _modelratios dictionary.

reduced_chisq([min])

The computed reduced chisquare value(s).

remove_measurement(id)

Delete a measurement from the internal dictionary used to compute ratios.

run(**kwargs)

Run the full computation using all the observations added.

write_chisq([chi, rchi, overwrite])

Write the chisq and reduced-chisq data to a file
add_measurement(m)[source]

Add a Measurement to internal dictionary used to compute ratios. This measurement may be intensity units (erg \({\rm s}^{-1}\) \({\rm cm}^{-2}\)) or integrated intensity (K km/s).

Parameters:

m (Measurement.) – a Measurement instance to be added to this tool

chisq(min=False)[source]

The computed chisquare value(s).

Parameters:

min (bool) – If True return the minimum reduced \(\chi^2\). In the case of map inputs this will be a spatial map of mininum \(\chi^2\). If False with map inputs the entire \(\chi^2\) hypercube is returned. If True with single pixel inputs, a single value is returned. If False with single pixel inputs, \(\chi^2\) as a function of density and radiation field is returned.

Return type:

Measurement

property density

The computed hydrogen nucleus density value(s).

Return type:

Measurement

property fit_result

The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables. For each pixel, there will be one instance of lmfit.minimizer.MinimizerResult.

Return type:

FitMap

property measurementIDs

The stored measurement IDs, which are strings.

Return type:

dict_keys

property measurements

The stored measurements as a dictionary with Measurement IDs as keys

Return type:

dict of Measurement

property modelset

The underlying ModelSet

property observed_ratios

The list of the observed line ratios that have been input so far.

Return type:

list of str

property radiation_field

The computed radiation field value(s).

Return type:

Measurement

property ratiocount

The number of ratios that match models available in the current ModelSet given the current set of measurements

Return type:

int

read_models(unit=Unit(dimensionless))[source]

Given a list of measurement IDs, find and open the FITS files that have matching ratios and populate the _modelratios dictionary. Uses pdrtpy.measurement.Measurement as a storage mechanism.

param m:

list of measurement IDS (string)

type m:

list

param unit:

units of the data

type unit:

string or astropy.Unit

reduced_chisq(min=False)[source]

The computed reduced chisquare value(s).

Parameters:

min (bool) – If True return the minimum reduced \(\chi_\nu^2\). In the case of map inputs this will be a spatial map of mininum \(\chi_\nu^2\). If False with map inputs the entire \(\chi_\nu^2\) hypercube is returned. If True with single pixel inputs, a single value is returned. If False with single pixel inputs, \(\chi_\nu^2\) as a function of density and radiation field is returned.

Return type:

Measurement

remove_measurement(id)[source]

Delete a measurement from the internal dictionary used to compute ratios.

Parameters:

id (str) – the measurement identifier

Raises:

KeyError – if id not in existing Measurements

run(**kwargs)[source]

Run the full computation using all the observations added. This will check compatibility of input observations (e.g., beam parameters, coordinate types, axes lengths) and raise exceptions if the observations don’t match each other.

param mask:
Indicate how to mask image observations (Measurements) before computing the density

and radiation field. Possible values are:

[‘mad’, multiplier] - compute standard deviation using median absolute deviation (astropy.mad_std),

and mask out values between +/- multiplier*mad_std. Example: [‘mad’,1.0]

[‘data’, (low,high)] - mask based on data values, mask out data between low and high

[‘clip’, (low,high)] - mask based on data values, mask out data below low and above high

[‘error’, (low,high)] - mask based on uncertainty plane, mask out data where the corresponding error pixel value

is below low or above high

None - Don’t mask data. This is the default.

type mask:

list or None

param method:

the fitting method to be used. The default is ‘leastsq’, which is Levenberg-Marquardt least squares. For other options see https://lmfit-py.readthedocs.io/en/latest/fitting.html#fit-methods-table

type method:

str

param nan_policy:

Specifies action if fit returns NaN values. One of: * ’raise’ : a ValueError is raised [Default] * ’propagate’ : the values returned from userfcn are un-altered * ’omit’ : non-finite values are filtered

type nan_policy:

str

raises Exception:

if no models match the input observations, observations are not compatible, or on unrecognized parameters, or NaN encountered.

property table

Construct the table of input Measurements, and if the fit has been run, the density, radiation field, and \(\chi^2\) values

Return type:

astropy.table.Table

write_chisq(chi='chisq.fits', rchi='rchisq.fits', overwrite=True)[source]

Write the chisq and reduced-chisq data to a file

Parameters:

  • chi (str) – FITS file to write the chisq map to.

  • rchi (str) – FITS file to write the reduced chisq map to.

FitMap

When fitting either single pixels or spatial maps, the fit results are stored per pixel in an NDData object that contains ModelResult objects for H2ExcitationFit or MinimizerResult objects for LineRatioFit. The user can thus examine in detail the fit at any pixel.

class pdrtpy.tool.fitmap.FitMap(data, *args, **kwargs)[source]

Bases: NDData

Attributes:
data

~numpy.ndarray-like : The stored dataset.

mask

any type : Mask for the dataset, if any.

meta

dict-like : Additional meta information about the dataset.

name

The name of this FitMap

psf

Image representation of the PSF for the dataset.

uncertainty

any type : Uncertainty in the dataset, if any.

unit

~astropy.units.Unit : Unit for the dataset, if any.

wcs

any type : A world coordinate system (WCS) for the dataset, if any.

Methods

get_pixel(world_x, world_y)

Return the nearest pixel coordinates to the input world coordinates
get_pixel(world_x, world_y)[source]

Return the nearest pixel coordinates to the input world coordinates

Parameters:

  • world_x (float) – The horizontal world coordinate

  • world_y (float) – The vertical world coordinate

property name

The name of this FitMap :rtype: str