Skip to content

Refactor Tutorial 1 to include instructions for specreduce #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Refactor Tutorial 1 to include instructions for specreduce #3

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
224 changes: 222 additions & 2 deletions 1_SpectroscopicTraceTutorial.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -361,7 +361,7 @@
"We can use [moments](https://en.wikipedia.org/wiki/Moment_(mathematics)) to provide a different, possibly better, estimate of where the trace's center is.\n",
"The advantage of moment analysis is that we're using all of the data to estimate the vertical position, not just the single brightest value, which is what we used above.\n",
"\n",
"Note that we need to subtract off the background to avoid a bias toward the center, so we use the median of the whole image as our background estimate.\n",
"Note that we need to subtract off the background to avoid a bias toward the center, so we use the median of the whole image as our background estimate. Specreduce provides more sophisticated background estimation, but we adopt a simple constant estimate here.\n",
"\n",
"(the first-order moment is the intensity-weighted mean position, \n",
"$$m_1 = \\frac{\\Sigma_i x_i f(x_i)}{\\Sigma_i f(x_i)}$$\n",
Expand Down Expand Up @@ -784,6 +784,69 @@
"plt.ylabel(\"Residual (data-model)\");"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3.1 Tracing - with specreduce\n",
"\n",
"The [``specreduce``](https://github.com/astropy/specreduce) package provides wrapper utilities to make the whole process faster and simpler.\n",
"\n",
"All of the above tools are integrated into the [``FitTrace``](https://specreduce.readthedocs.io/en/latest/api/specreduce.tracing.FitTrace.html) class.\n",
"\n",
"In the example below, we skip some of the above data examination steps and we solve the noise problems differently:\n",
"\n",
" * Instead of rejecting pixels at x>1200, we bin the data into 20 bins along the X-axis. This approach provides enough independent bins to fit a 3rd-order polynomial. Binning increases the signal-to-noise so that we get a good fit.\n",
" * We already know that a 3rd order polynomial is a good fit, so we specify that immediately."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from specreduce.tracing import FitTrace"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"trace = FitTrace(image_array, bins=20, trace_model=Polynomial1D(3))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Show the trace overlaid on the image cutout around the spectrum\n",
"plt.imshow(\n",
" image_array[425:475, :], extent=[0, image_array.shape[1], 425, 475], cmap=\"gray\"\n",
")\n",
"plt.plot(xvals, trace.trace, \"r-\", linewidth=2, alpha=0.8, label=\"Fitted trace\")\n",
"plt.gca().set_aspect(10) # stretch in y-direction for better visibility\n",
"plt.xlabel(\"X position (pixels)\")\n",
"plt.ylabel(\"Y position (pixels)\")\n",
"plt.title(\"Fitted trace overlaid on spectrum\")\n",
"plt.legend()\n",
"\n",
"plt.figure(figsize=(8, 4))\n",
"plt.plot(xvals[~bad_moments], weighted_yaxis_values[~bad_moments], \"x\", alpha=0.5)\n",
"plt.plot(xvals, trace.trace, color=\"r\");"
]
},
{
"cell_type": "markdown",
"metadata": {
Expand All @@ -801,7 +864,9 @@
"source": [
"Now we can extract the data along that trace.\n",
"\n",
"We want to take a \"profile\" of the trace to see how many pixels on either side of the line we should include."
"We want to take a \"profile\" of the trace to see how many pixels on either side of the line we should include.\n",
"\n",
"The 'trace profile' is effectively the point spread function (PSF) of the telescope."
]
},
{
Expand Down Expand Up @@ -1081,6 +1146,115 @@
"Note that the Gaussian model and the direct trace yield nearly identical results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5.1 Spectral Extraction with specreduce\n",
"\n",
"We can extract the 1D spectrum using specreduce's extraction methods. We'll demonstrate both boxcar and optimal (Horne) extraction:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from specreduce.extract import BoxcarExtract, HorneExtract"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5.1.1. Boxcar Extraction\n",
"\n",
"Boxcar extraction sums all the pixels within a rectangular aperture centered on the trace. This is the simplest extraction method:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Perform boxcar extraction on the background-subtracted image\n",
"bg_subtracted = image_array - background\n",
"boxcar_extract = BoxcarExtract(bg_subtracted, trace, width=15)\n",
"boxcar_spectrum = boxcar_extract.spectrum"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(boxcar_spectrum.flux, label=\"Boxcar extracted spectrum\")\n",
"plt.xlabel(\"Pixel position\")\n",
"plt.ylabel(\"Flux\")\n",
"plt.title(\"1D Spectrum from Boxcar Extraction\")\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5.1.2: Optimal (Horne) Extraction\n",
"\n",
"Optimal (aka Horne) extraction weights pixels based on their signal-to-noise ratio, which provides better results for noisy data. As above, it uses a weighting profile, but it also _requires_ specifying an uncertainty per pixel so it can compute proper weights.\n",
"\n",
"Ideally, we would compute the uncertainty per pixel and attach that to our measurement, but we will skip that process here and revisit it in another notebook **TODO** and for the moment assume a constant uncertainty per pixel."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# For optimal extraction, we can use the background-subtracted image directly\n",
"# Note: For real optimal extraction, variance estimates should be provided\n",
"\n",
"from astropy.nddata import NDData, VarianceUncertainty\n",
"\n",
"bg_subtracted_nddata = NDData(\n",
" bg_subtracted,\n",
" uncertainty=VarianceUncertainty(np.ones_like(bg_subtracted)),\n",
" unit=u.adu,\n",
")\n",
"\n",
"optimal_extract = HorneExtract(bg_subtracted_nddata, trace)\n",
"optimal_spectrum = optimal_extract.spectrum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5.1.3 Compare Extraction Methods\n",
"\n",
"Let's compare the boxcar and optimal extraction results:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(12, 6))\n",
"plt.plot(boxcar_spectrum.flux, label=\"Boxcar extraction\", alpha=0.8)\n",
"plt.plot(optimal_spectrum.flux, label=\"Optimal (Horne) extraction\", alpha=0.8)\n",
"plt.xlabel(\"Pixel position\")\n",
"plt.ylabel(\"Flux\")\n",
"plt.title(\"Comparison of Extraction Methods\")\n",
"plt.legend()\n",
"plt.grid(True, alpha=0.3);"
]
},
{
"cell_type": "markdown",
"metadata": {
Expand Down Expand Up @@ -1207,6 +1381,52 @@
"plt.plot(spectrum2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6.1 Specreduce version"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"trace2 = FitTrace(image_array_2, bins=20, trace_model=Polynomial1D(3))\n",
"spectrum2boxcar = BoxcarExtract(image_array_2, trace2, width=15).spectrum"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Show the trace overlaid on the image cutout around the spectrum\n",
"plt.imshow(image_array_2[470:520, :], extent=[0, 1600, 470, 520], cmap=\"gray\")\n",
"plt.plot(xvals, trace2.trace, \"r-\", linewidth=2, alpha=0.8, label=\"Fitted trace\")\n",
"plt.gca().set_aspect(10) # stretch in y-direction for better visibility\n",
"plt.xlabel(\"X position (pixels)\")\n",
"plt.ylabel(\"Y position (pixels)\")\n",
"plt.title(\"Fitted trace overlaid on spectrum\")\n",
"plt.legend()\n",
"\n",
"plt.figure(figsize=(8, 4))\n",
"plt.plot(xvals[~bad_moments], weighted_yaxis_values2[~bad_moments], \"x\", alpha=0.5)\n",
"plt.plot(xvals, trace2.trace, color=\"r\");"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(spectrum2boxcar.flux)"
]
},
{
"cell_type": "markdown",
"metadata": {
Expand Down