{ "cells": [ { "cell_type": "markdown", "id": "f8fb7695-bc47-4c48-b98a-7b6a548ad05b", "metadata": {}, "source": [ "# Gaussian Fitting\n", "This notebook shows how to use \n", "[`astropy`]() and \n", "[`specutils`]() to fit a Gaussian line profile to a\n", "[Spectrum](https://dysh.readthedocs.io/en/latest/reference/modules/dysh.spectra.html#dysh.spectra.spectrum.Spectrum).\n", "The data comes from the \n", "[calseq tutorial](https://dysh.readthedocs.io/en/latest/tutorials/examples/calseq.html).\n", "\n", "## Dysh commands\n", "\n", "The following dysh commands are introduced (leaving out all the function arguments):\n", "\n", " filename = dysh_data()\n", " ta = Spectrum.read()\n", " p1 = ta.plot()\n", " p1.oshow(ta1)\n", "\n", "## Loading Modules\n", "We start by loading the modules we will use for the data reduction. \n" ] }, { "cell_type": "code", "execution_count": 1, "id": "ffc73528-c29b-4573-93e1-04abdffc49b6", "metadata": {}, "outputs": [], "source": [ "# These are the modules we will use for loading the data and fitting.\n", "from dysh.spectra import Spectrum\n", "from specutils.fitting import fit_lines\n", "from astropy.modeling import models\n", "from astropy import units as u\n", "from dysh.log import init_logging\n", "\n", "# These modules are used for file I/O\n", "from dysh.util.files import dysh_data\n", "from pathlib import Path" ] }, { "cell_type": "markdown", "id": "37ab6793-4d78-4657-aeef-253f30992286", "metadata": {}, "source": [ "## Setup\n", "We start the dysh logging, so we get more information about what is happening.\n", "This is only needed if working on a notebook.\n", "If using the CLI through the ``dysh`` command, then logging is setup for you." ] }, { "cell_type": "code", "execution_count": 2, "id": "2e5d9a79-ae71-4426-90f4-e9847baeae81", "metadata": {}, "outputs": [], "source": [ "init_logging(2)\n", "\n", "# also create a local \"output\" directory where temporary notebook files can be stored.\n", "output_dir = Path.cwd() / \"output\"\n", "output_dir.mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "id": "fc410cba-61d0-485a-9453-37a35c49e374", "metadata": {}, "source": [ "## Data Retrieval\n", "\n", "Download the example FITS spectrum, if necessary." ] }, { "cell_type": "code", "execution_count": 3, "id": "9b7b319b-2045-4e8f-8d60-06692046aa2a", "metadata": {}, "outputs": [], "source": [ "filename = dysh_data(example=\"nod-W/outputs/M82_ifnum_3_polavg.fits\")" ] }, { "cell_type": "markdown", "id": "b96f9f20-643c-4c0e-8293-d5597ccf85d2", "metadata": {}, "source": [ "## Data Loading\n", "\n", "We use \n", "[Spectrum.read](https://dysh.readthedocs.io/en/latest/reference/modules/dysh.spectra.html#dysh.spectra.spectrum.Spectrum.read)\n", "to load the spectrum." ] }, { "cell_type": "code", "execution_count": 4, "id": "a3a77e4c-cff6-4e23-b340-8dbd93800197", "metadata": {}, "outputs": [], "source": [ "spec = Spectrum.read(filename, format=\"fits\")" ] }, { "cell_type": "markdown", "id": "a73aecc6-d737-4a84-8f7a-34935c02ad25", "metadata": {}, "source": [ "The loaded spectrum is now a \n", "[Spectrum](https://dysh.readthedocs.io/en/latest/reference/modules/dysh.spectra.html#dysh.spectra.spectrum.Spectrum),\n", "with all its methods.\n", "We can plot it, using \n", "[plot](https://dysh.readthedocs.io/en/latest/reference/modules/dysh.spectra.html#dysh.spectra.spectrum.Spectrum.plot)." ] }, { "cell_type": "code", "execution_count": 5, "id": "e19cdac1-29a0-4350-a5e2-eff395a5f4e4", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a2a59394231945bb99145dcd0e53aba3", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sp_plt = spec.plot(xaxis_unit=\"GHz\")" ] }, { "cell_type": "markdown", "id": "c6fae124-45c1-4c85-b6ce-c3a8a0d1fad5", "metadata": {}, "source": [ "Before fitting, lets remove the edges of the spectrum.\n", "We only keep channels between 88.2 and 89.5 GHz." ] }, { "cell_type": "code", "execution_count": 6, "id": "8951cf5d-f2e8-403c-a81a-aadc64778009", "metadata": {}, "outputs": [], "source": [ "spec = spec[88.2*u.GHz:89.5*u.GHz]" ] }, { "cell_type": "markdown", "id": "d9c1047f-6180-4917-b0d0-942e6799f0ab", "metadata": {}, "source": [ "Plot again." ] }, { "cell_type": "code", "execution_count": 7, "id": "72c0c1e5-8842-4093-a374-d5d6b23ee1d9", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a97326cde91e46899eaf7befe244554a", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "spec_plt = spec.plot(xaxis_unit=\"GHz\")" ] }, { "cell_type": "markdown", "id": "11b34992-b667-466f-9d0f-a03ec08bdfc1", "metadata": {}, "source": [ "Now we are ready to start fitting.\n", "\n", "## Gaussian Fitting\n", "\n", "Here we show how to fit a Gaussian profile using [`astropy`]() and [`specutils`]().\n", "There are other options available, but we won't cover them all.\n", "\n", "### Fitting a Single Gaussian\n", "\n", "We start fitting a single Gaussian line profile to the line on the left, at about 88.55 GHz.\n", "We create a model, with starting values, and then use the [`specutils.fitting.fit_lines`](https://specutils.readthedocs.io/en/stable/api/specutils.fitting.fit_lines.html#specutils.fitting.fit_lines) function to fit the model to the spectrum.\n", "After the fit is done, we evaluate the best fit model over the spectral axis of the spectrum so we can plot it on top of the data." ] }, { "cell_type": "code", "execution_count": 8, "id": "0b135e3e-d9e0-434b-8fab-c050ac6a346f", "metadata": {}, "outputs": [], "source": [ "g_init = models.Gaussian1D(amplitude=0.5*u.K, mean=88.55*u.GHz, stddev=0.1*u.GHz)\n", "g_fit = fit_lines(spec, g_init)\n", "y_fit = g_fit(spec.spectral_axis)" ] }, { "cell_type": "markdown", "id": "bd3599c1-4ea8-401f-bcbe-30d87d1f91b5", "metadata": {}, "source": [ "Now plot it.\n", "We assign a \"group id\" (gid), so we can remove the line after." ] }, { "cell_type": "code", "execution_count": 9, "id": "04541b19-a46e-4594-aa5a-87cffee77e6e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spec_plt.axis.plot(spec.spectral_axis.to(\"GHz\"), y_fit, gid=\"1gauss\")\n", "spec_plt.figure # This will show the figure here." ] }, { "cell_type": "markdown", "id": "b00bc4f0-337e-4c60-bae8-454aa92abc07", "metadata": {}, "source": [ "#### Extracting Best Fit Parameters\n", "\n", "The best fit parameters are now properties of the output best fit model.\n", "We can access them through the parameter names, for the Gaussian, amplitude, mean and stddev.\n", "For example, the line amplitude" ] }, { "cell_type": "code", "execution_count": 10, "id": "9d0306ce-52ba-4c7c-bdc8-f6e8f982fb2f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Parameter('amplitude', value=0.3670704696308813, unit=K)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g_fit.amplitude" ] }, { "cell_type": "markdown", "id": "aa7481b3-6d68-46c9-a0f6-75a654b12896", "metadata": {}, "source": [ "and its formal error" ] }, { "cell_type": "code", "execution_count": 11, "id": "51219229-2e5f-4262-974b-c39e084f9bde", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(0.006236258674343534)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g_fit.amplitude.std" ] }, { "cell_type": "markdown", "id": "ac5067bc-94e5-4b56-854c-218960629975", "metadata": {}, "source": [ "#### Fitting Two Gaussians\n", "\n", "We can also combine models.\n", "In this case, we add a second Gaussian line profile to fit the second line at about 89.15 GHz." ] }, { "cell_type": "code", "execution_count": 12, "id": "3a4235fb-690f-4341-bba8-b3c7dc226779", "metadata": {}, "outputs": [], "source": [ "g_init2 = models.Gaussian1D(amplitude=0.6*u.K, mean=89.15*u.GHz, stddev=0.1*u.GHz)\n", "g_fit = fit_lines(spec, g_init + g_init2)\n", "y_fit = g_fit(spec.spectral_axis)" ] }, { "cell_type": "markdown", "id": "196d90e3-3898-4d13-a3f4-23008e41bac1", "metadata": {}, "source": [ "Before plotting the new model with two Gaussians, we remove the previous line." ] }, { "cell_type": "code", "execution_count": 13, "id": "9e96f4a3-f56a-418c-8df0-4e9f76fb28a2", "metadata": {}, "outputs": [], "source": [ "spec_plt.clear_lines(\"1gauss\")" ] }, { "cell_type": "markdown", "id": "72e75de1-a4e5-4fcc-b53f-e9c84ea2a11e", "metadata": {}, "source": [ "We repeat the plotting with the new model results.\n", "This time, we add a label to the best fit model line (`label=\"best fit\"`)." ] }, { "cell_type": "code", "execution_count": 14, "id": "071d8714-15ab-412b-b4f9-afc04ec3a61b", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a97326cde91e46899eaf7befe244554a", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "spec_plt.axis.plot(spec.spectral_axis.to(\"GHz\"), y_fit, label=\"best fit\")\n", "spec_plt.axis.legend()\n", "spec_plt.figure\n", "spec_plt.show()" ] }, { "cell_type": "markdown", "id": "40cb4d12-b356-49b3-881b-5d936c619eec", "metadata": {}, "source": [ "You can find more examples of how to use [``specutils.fitting.fit_lines``](https://specutils.readthedocs.io/en/stable/api/specutils.fitting.fit_lines.html#specutils.fitting.fit_lines) in the [`specutils` documentation](https://specutils.readthedocs.io/en/stable/fitting.html).\n", "\n", "In the case of a compound model, the model parameters get assigned an underscore and a number to differentiate them.\n", "So the amplitude of the first Gaussian is stored in `g_fit.amplitude_0` and that of the second Gaussian in `g_fit.amplitude_1`" ] }, { "cell_type": "code", "execution_count": 15, "id": "f9cfb0c9-f494-470d-98f9-10a3aa133cbc", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Parameter('amplitude', value=0.3670705202212676, unit=K)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g_fit.amplitude_0" ] }, { "cell_type": "code", "execution_count": 16, "id": "eb348b4b-02cb-4cc5-92c6-b362e3dd4e2e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Parameter('amplitude', value=0.5281505181498451, unit=K)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g_fit.amplitude_1" ] }, { "cell_type": "markdown", "id": "3f6c0844-f30d-4a7f-8c9c-bb70e8e5d924", "metadata": {}, "source": [ "#### Other Profiles and Further Examples\n", "Additional profiles and models are available through the [`astropy.modeling` submodule](https://docs.astropy.org/en/stable/modeling/predef_models1D.html).\n", "\n", "The [`astropy`]() documentation also provides examples of how to [fit models with constraints](https://docs.astropy.org/en/stable/modeling/example-fitting-constraints.html), and [examples of how to fit models](https://docs.astropy.org/en/stable/modeling/fitting.html), without using `specutils`.\n", "\n" ] }, { "cell_type": "markdown", "id": "be79ada3-6dc5-4c01-8342-f61a22b48e4e", "metadata": {}, "source": [ "## Final Stats\n", "Finally, at the end we compute some statistics over a spectrum, merely as a checksum if the notebook is reproducible.\n" ] }, { "cell_type": "code", "execution_count": 17, "id": "75ba811b-c13b-4a4b-9502-2374e21fe001", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "16:02:08.300 I Note: found 27 NaN (masked) values\n", "16:02:08.301 I rms is OK \n" ] } ], "source": [ "spec.check_stats(0.11606252 * u.K)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.9" } }, "nbformat": 4, "nbformat_minor": 5 }