From a1d242ec96bc182d159f6d7588dd3b402ed1b045 Mon Sep 17 00:00:00 2001 From: "Adam Ginsburg (keflavich)" Date: Wed, 25 Jun 2025 20:11:14 -0400 Subject: [PATCH] refactor tutorial to include (but not exclusively use) specreduce --- 1_SpectroscopicTraceTutorial.ipynb | 224 ++++++++++++++++++++++++++++- 1 file changed, 222 insertions(+), 2 deletions(-) diff --git a/1_SpectroscopicTraceTutorial.ipynb b/1_SpectroscopicTraceTutorial.ipynb index fbd8c4f..fb2c227 100755 --- a/1_SpectroscopicTraceTutorial.ipynb +++ b/1_SpectroscopicTraceTutorial.ipynb @@ -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", @@ -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": { @@ -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." ] }, { @@ -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": { @@ -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": {