Analysis Tools: Fit models to data

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: ABC

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? (i.e., have 2 spatial axes)
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:

bool: True if the observational inputs are spatial maps, False otherwise.

property has_scalar#

Are the Measurements used scalars.

Returns:

bool: True if the observational inputs are scalars, False otherwise.

property has_vectors#

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

Returns:

bool: True if the observational inputs are a vector, False otherwise.

abstractmethod 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, visual extinction, and ortho-to-para ratio in \(H_2\) excitation diagrams. With a two temperature model assumed, the fit will find \(T_{hot}, T_{cold}, N_{hot}(H_2), N_{cold}(H_2),\) and optionally \(A_v\) or \(OPR\). The base class ExcitationFit can be used to create a tool to fit a different molecule. Single temperature fitting is also available.

Similarly, there are tools to fit \(^{12}CO, ^{13}CO, ^{12}C^{18}O, ^{13}C^{18}O, {\rm and}~ CH^{+}\).

class pdrtpy.tool.excitation.BaseExcitationFit(molecule: BaseMolecule, measurements: dict | Measurement = None)[source]#

Bases: ToolBase

Base class for creating excitation fitting tools for various species.

Parameters:

moleculeBaseMolecule: The molecule whose transitions will be fit.
measurementsMeasurement or dict, optional.: Input measurements to be fit. If input is a dictionary of measurements, the keys must Measurement identifiers. The default is None.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

add_measurement(m: Measurement)[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:

mMeasurement: A measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\).

property av#

The visual extinction

Returns:

Measurement: The fitted Av if it was determined in the fit, otherwise 0.

property av_fitted#

Was the visual extinction fitted?

Returns:

bool: True if Av was fitted, False if not.

average_column_density(position=None, size=None, norm=True, unit=Unit('1 / cm2'), line=True, clip=None)[source]#

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

Parameters:

positiontuple: 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.
sizeint or 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).
normbool: 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\).
unitstr or astropy.units.Unit: The units in which to return the column density. Default: \({\rm cm}^{-2}\).
linebool: If True, the returned dictionary index is the Line name, otherwise it is the upper state \(J\) number.
clipastropy.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:

dict: Dictionary of column density Measurements, with keys as \(J\) number or Line name.

property cold_colden#

The fitted cold gas total column density

Returns:

Measurement

colden(component)[source]#

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

Parameters:

componentstr: ‘hot’, ‘cold’, or ‘total’.

Returns:

Measurement

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

The computed upper state column densities of stored intensities

Parameters:

normbool: If True, normalize the column densities by the statistical weight of the upper state, \(g_u\). Default: False.
unitstr or astropy.units.Unit: The units in which to return the column density. Default: \({\rm cm}^{-2}\).
linebool: If True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: True.

Returns:

dict: Dictionary of column densities indexed by upper state \(J\) number or Line name.

energies(line=True)[source]#

Upper state energies of stored intensities, in K.

Parameters:

linebool: If True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: True.

Returns:

dict: Dictionary indexed by upper state \(J\) number or Line name.

property extinction_model#

The extinction law used when fitting for visual extinction, \(A_v\).

Returns:

modelBaseExtModel or astropy.modeling.Model: The model to be used to calculate dust extinction.

property fit_result#

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

Returns:

lmfit.model.ModelResult

gu(id, opr)[source]#

Get the upper state statistical weight \(g_u\) for the given transition identifier, 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:

idstr: The measurement identifier.
oprfloat: Ortho-to-para ratio.

Returns:

float

Raises:

KeyError: If id not in existing Measurements.

property hot_colden#

The fitted hot gas total column density

Returns:

Measurement

property intensities#

The stored intensities. See add_measurement()

Returns:

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:

coldenMeasurement: Upper state column density.

Returns:

Measurement: Optically thin intensity.

property molecule: BaseMolecule#

The molecule being fitted by this ExcitationFit

Returns:

Molecule: The molecule as represented by the Molecule class.

property numcomponents#

Number of temperature components in the fit

Returns:

int

property opr#

The ortho-to-para ratio (OPR)

Returns:

Measurement: The fitted OPR if it was determined in the fit, otherwise the canonical LTE OPR.

property opr_fitted#

Was the ortho-to-para ratio fitted?

Returns:

bool: True if OPR was fitted, False if canonical LTE value was used or this molecule’s OPR cannot vary.

remove_measurement(identifier: str)[source]#

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

Parameters:

identifierstr: The measurement identifier.

Raises:

KeyError: If identifier not in existing Measurements.

replace_measurement(m: Measurement)[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:

mMeasurement: A measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\).

run(position=None, size=None, fit_opr=False, fit_av=False, components=2, **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:

positiontuple

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.

sizeint or 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_oprbool

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.

fit_avbool

Whether to fit the visual extinction. If True, the Av will be varied to determine the best value. If False, the Av is fixed at zero.

workersint or None

Number of worker processes for parallel pixel fitting. None (default) runs serially. -1 uses all available CPUs. Any positive integer uses that many workers. Matches the LineRatioFit.run() API.

Performance note: each pixel is submitted as a separate task to ProcessPoolExecutor, so inter-process communication overhead is paid once per pixel. For the excitation fits in this package (n ≤ ~20 spectral lines, lightweight lmfit minimisation) the per-pixel compute time is short enough that parallel execution only outperforms serial on maps with roughly 5 000 or more valid (unmasked) pixels. On smaller maps the IPC overhead dominates and serial is faster.

emcee fitting is excluded from the parallel path regardless of this setting.

chunk_sizeint

Number of pixels batched into each parallel task (default: 32). Each worker process fits chunk_size pixels serially, so IPC serialisation overhead is paid once per chunk rather than once per pixel. Larger values reduce overhead further but coarsen progress-bar granularity and may cause load imbalance on the last chunk. Ignored when workers is None or when the fitting method is 'emcee'.

set_extinction_model(model)[source]#

Set the model to be used for fitting visual extinction, \(A_v\). This is typically a model from the dust_extinction package.

Parameters:

modelBaseExtModel or astropy.modeling.Model: The model to be used to calculate dust extinction.

Returns:

None.

property tcold#

The fitted cold gas excitation temperature

Returns:

Measurement

property temperature#

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

Returns:

dict

property thot#

The fitted hot gas excitation temperature

Returns:

Measurement

property total_colden#

The fitted total column density

Returns:

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:

intensityMeasurement: A measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\).
unitstr or astropy.units.Unit: The units in which to return the column density. Default: \({\rm cm}^{-2}\).

Returns:

Measurement: The column density.

wavelengths(line=True, units=False)[source]#

Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

Parameters:

linebool: If True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: True.
unitsbool: If True, values are returned with units as astropy Quantity.

Returns:

dict: Dictionary indexed by upper state \(J\) number or Line name.

class pdrtpy.tool.excitation.C13O18ExcitationFit(measurements: Measurement = None)[source]#

Bases: BaseExcitationFit

Tool for fitting temperatures, column densities, \(A_v\) from an \(^{13}C^{18}O\) excitation diagram. It takes as input a set of \(^{13}C^{18}O\) 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 \(A_v\). One needs at least 5 points to fit two temperatures and column densities (slope and intercept \(\times 2\)), though one could compute (not fit) them with only 4 points.

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

Parameters:

measurementslist of Measurement: Input \(^{13}C^{18}O\) measurements to be fit.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

class pdrtpy.tool.excitation.C13OExcitationFit(measurements: Measurement = None)[source]#

Bases: BaseExcitationFit

Tool for fitting temperatures, column densities, \(A_v\) from an \(^{13}C^{16}O\) excitation diagram. It takes as input a set of \(^{13}CO\) rovibrational line observations with errors represented as Measurement.

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

Parameters:

measurementslist of Measurement: Input \(^{13}CO\) measurements to be fit.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

class pdrtpy.tool.excitation.CHplusExcitationFit(measurements: Measurement = None)[source]#

Bases: BaseExcitationFit

Tool for fitting temperatures, column densities, \(A_v\), and ortho-to-para ratio(OPR) from an \(CH^{+}\) excitation diagram. It takes as input a set of \(CH^{+}\) 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}\). 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.

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

Parameters:

measurementslist of Measurement: Input \(CH^{+}\) measurements to be fit.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

class pdrtpy.tool.excitation.CO18ExcitationFit(measurements: Measurement = None)[source]#

Bases: BaseExcitationFit

Tool for fitting temperatures, column densities, \(A_v\) from an \(^{12}C^{18}O\) excitation diagram. It takes as input a set of \(^{12}C^{18}O\) rovibrational line observations with errors represented as Measurement.

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

Parameters:

measurementslist of Measurement: Input \(^{12}C^{18}O\) measurements to be fit.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

class pdrtpy.tool.excitation.COExcitationFit(measurements: Measurement = None)[source]#

Bases: BaseExcitationFit

Tool for fitting temperatures, column densities, \(A_v\) from an \(^{12}C^{16}O\) excitation diagram. It takes as input a set of \(^{12}CO\) 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 \(A_v\). One needs at least 5 points to fit two temperatures and column densities (slope and intercept \(\times 2\)), though one could compute (not fit) them with only 4 points. To additionally fit \(A_v\), one should have 6 points (5 degrees of freedom).

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

Parameters:

measurementslist of Measurement: Input \(^{12}CO\) measurements to be fit.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

class pdrtpy.tool.excitation.H2ExcitationFit(measurements: Measurement = None)[source]#

Bases: BaseExcitationFit

Tool for fitting temperatures, column densities, \(A_v\), 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 \(A_v\) or \(OPR\). One needs at least 5 points to fit two temperatures and column densities (slope and intercept \(\times 2\)), though one could compute (not fit) them with only 4 points. To additionally fit \(A_v\) or \(OPR\), one should have 6 points (5 degrees of freedom).

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

Parameters:

measurementslist of Measurement: Input \(H_2\) measurements to be fit.

Attributes:

av: The visual extinction
av_fitted: Was the visual extinction fitted?
cold_colden: The fitted cold gas total column density
extinction_model: The extinction law used when fitting for visual extinction, \(A_v\).
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? (i.e., have 2 spatial axes)
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.
molecule: The molecule being fitted by this ExcitationFit
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 identifier, 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_av, ...])	Fit the \(log N_u-E\) diagram with two excitation temperatures, a `hot` \(T_{ex}\) and a `cold` \(T_{ex}\).
`set_extinction_model`(model)	Set the model to be used for fitting visual extinction, \(A_v\).
`upper_colden`(intensity, unit)	Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
`wavelengths`([line, units])	Wavelengths of transitions, in micron (assumed unit using Roueff et al table)

Fitting Line Ratios For Density and Radiation Field#

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, measurements=None)[source]#

Bases: ToolBase

Tool to fit observations of intensity ratios to a set of PDR models.

Takes as input a set of observations with errors represented as Measurement and a ModelSet for the models to fit. Observations should be spectral line or continuum intensities, either spatial maps or single pixel values at the same spatial resolution.

At least 3 observations are needed to make at least 2 ratios. Once the fit is done, LineRatioPlot can be used to view the results.

Parameters:

modelsetModelSet: The set of PDR models to use for fitting.
measurementslist or dict of Measurement, optional: Input measurements to be fit. If dict, keys should be Measurement identifiers.

Attributes:

density: The computed hydrogen nucleus density value(s).
fit_result: The result of the fitting procedure, including fit statistics, variable values and uncertainties, and correlations.
has_maps: Are the Measurements used map-based? (i.e., have 2 spatial axes)
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 the 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 with matching ratios and populate `_modelratios`.
`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 added observations.
`write_chisq`([chi, rchi, overwrite])	Write the chisq and reduced-chisq data to a file.

add_measurement(m)[source]#

Add a Measurement to the internal dictionary used to compute ratios.

The measurement may be in intensity units (\({\rm erg~s}^{-1}\) \({\rm cm}^{-2}\)) or integrated intensity (K km/s).

Parameters:

mMeasurement: A Measurement instance to be added to this tool.

chisq(min=False)[source]#

The computed chisquare value(s).

Parameters:

minbool, optional: If True, return the minimum \(\chi^2\). For map inputs, returns a spatial map of the minimum; for single-pixel inputs, returns a scalar. If False, returns the full \(\chi^2\) as a function of density and radiation field. Default: False.

Returns:

Measurement

property density#

The computed hydrogen nucleus density value(s).

Returns:

Measurement

property fit_result#

The result of the fitting procedure, including fit statistics, variable values and uncertainties, and correlations. One lmfit.minimizer.MinimizerResult per pixel.

Returns:

FitMap

property measurementIDs#

The stored measurement IDs, which are strings.

Returns:

dict_keys

property measurements#

The stored measurements as a dictionary with Measurement IDs as keys.

Returns:

dict of Measurement

property modelset#: The underlying ModelSet

property observed_ratios#

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

Returns:

list of str

property radiation_field#

The computed radiation field value(s).

Returns:

Measurement

property ratiocount#

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

Returns:

int

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

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

Uses pdrtpy.measurement.Measurement as a storage mechanism.

Parameters:

unitstr or astropy.units.Unit, optional: Units of the data.

reduced_chisq(min=False)[source]#

The computed reduced chisquare value(s).

Parameters:

minbool, optional: If True, return the minimum \(\chi_\nu^2\). For map inputs, returns a spatial map of the minimum; for single-pixel inputs, returns a scalar. If False, returns the full \(\chi_\nu^2\) as a function of density and radiation field. Default: False.

Returns:

Measurement

remove_measurement(id)[source]#

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

Parameters:

idstr: The measurement identifier.

Raises:

KeyError: If id not in existing Measurements.

run(**kwargs)[source]#

Run the full computation using all the added observations.

Checks compatibility of input observations (beam parameters, coordinate types, axes lengths) and raises exceptions if they don’t match.

Parameters:

masklist or None, optional

Indicate how to mask image observations before computing density and radiation field. Possible values:

[‘mad’, multiplier] — mask values between +/- multiplier*mad_std
[‘data’, (low, high)] — mask data values between low and high
[‘clip’, (low, high)] — mask data values outside [low, high]
[‘error’, (low, high)] — mask where error pixel is below low or above high
None — no masking (default)

radiation_field_rangeQuantity or sequence, optional

Restrict the fit to radiation field values within [lower, upper] (inclusive), excluding unphysical regimes. Accepts either a length-2 Quantity (e.g. [100, 1000]*habing_unit) or a length-2 sequence whose elements are each None or a scalar Quantity (e.g. [None, 1000*habing_unit] for an upper limit only). None on a side means no limit there. Note the unit must be inside the brackets for a one-sided limit, since None*unit is not allowed. Default: None (no restriction).

density_rangeQuantity or sequence, optional

Restrict the fit to density values within [lower, upper] (inclusive). Same accepted forms as radiation_field_range, e.g. [1e3, 1e4]/u.cm**3 or [None, 1e4/u.cm**3]. None on a side means no limit there. Default: None (no restriction).

methodstr, optional

Fitting method. Default: ’leastsq’ (Levenberg-Marquardt). See https://lmfit-py.readthedocs.io/en/latest/fitting.html#fit-methods-table.

nan_policystr, optional

Action if fit returns NaN. One of ’raise’ (default), ’propagate’, ’omit’.

workersint or None, optional

Worker processes for parallel pixel fitting. None uses serial fitting. -1 uses all CPUs. Ignored for single-pixel fits, emcee, and when joint_fit is not None. Default: None.

joint_fitstr or None, optional

Controls joint pixel fitting:

None (default): serial or parallel per-pixel fitting.
’hybrid’: joint scipy fit + single-pixel re-fits for boundary pixels. ~3× faster than workers=-1.
’fast’: joint scipy fit only, no post-processing. ~11× faster but ~7–9% accuracy loss near boundaries.

All joint-fit modes use TRF. Ignored for single-pixel fits and emcee.

Raises:

Exception: If no models match the input observations, observations are incompatible, parameters are unrecognized, or NaN is encountered.
ValueError: If radiation_field_range or density_range is malformed, has the wrong units, or specifies a window that does not overlap the model grid.

property table#

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

Returns:

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:

chistr, optional: FITS file to write the chisq map to. Default: "chisq.fits".
rchistr, optional: FITS file to write the reduced chisq map to. Default: "rchisq.fits".

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

A class that can store fit objects in a data array but has all the nice WCS properties of NDData.

Parameters:

dataarray-like: The data set, an array of lmfit.model.ModelResult or lmfit.minimizer.MinimizerResult.
namestr, optional: An identifying name for this object.

Attributes:

data: ndarray-like : The stored dataset.
mask: any type : Mask for the dataset, if any.
meta
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: 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 x,y.
`get_pixel_from_coord`(coord)	Return the nearest pixel coordinates to the input world coordinates.
`get_skycoord`(x, y)	Return the Sky Coordinate corresponding to the input pixel coordinates.
`get_world`(x, y)	Return the world coordinates corresponding to the input pixel coordinates.

get_pixel(world_x, world_y)[source]#

Return the nearest pixel coordinates to the input world coordinates x,y.

The pixel values will be rounded to the nearest integer.

Parameters:

world_xfloat: The horizontal world coordinate.
world_yfloat: The vertical world coordinate.

get_pixel_from_coord(coord)[source]#

Return the nearest pixel coordinates to the input world coordinates.

The pixel values will be rounded to the nearest integer.

Parameters:

coordSkyCoord: The world coordinate.

get_skycoord(x, y)[source]#

Return the Sky Coordinate corresponding to the input pixel coordinates.

Parameters:

xfloat: The horizontal pixel coordinate.
yfloat: The vertical pixel coordinate.

get_world(x, y)[source]#

Return the world coordinates corresponding to the input pixel coordinates.

Parameters:

xfloat: The horizontal pixel coordinate.
yfloat: The vertical pixel coordinate.

property name#

The name of this FitMap.

Returns:

str

Analysis Tools: Fit models to data

Contents

Analysis Tools: Fit models to data#

ToolBase#

Excitation Diagram Fitting#

Fitting Line Ratios For Density and Radiation Field#

FitMap#