From d91df4e041e467e1245f3caf0ad35a5da4f826a1 Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Thu, 17 Mar 2022 10:21:25 -0400 Subject: [PATCH] support statistic='median' * defaults to 'average' * median does not support partial-pixel-weights --- .../jwst_boxcar/boxcar_extraction.ipynb | 48 +++++++++++++++++-- specreduce/background.py | 28 ++++++++++- 2 files changed, 72 insertions(+), 4 deletions(-) diff --git a/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb b/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb index f0d6eb7d..8f631913 100644 --- a/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb +++ b/notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb @@ -508,13 +508,55 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Trace" + "Note that when using median, partial pixel weights are not supported:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'slit[0] slice')" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAFLCAYAAABSuvQBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUk0lEQVR4nO3df7DldX3f8dc7u8sSfhhACV2BiBGiQ9KAnS2aaluFmKC1hTQ2lUlT0jGzyUxopWN/EDNtdabMaCdqkoljZw1UpjGiRRwZJyMipTXOtNRFiQKrBS0qW2AxgqwxIuC7f5wvkxu6y73cH3s+3PN4zNy55/vrnDczZ77Lc8/3fLe6OwAAAMzXD8x7AAAAAMQZAADAEMQZAADAAMQZAADAAMQZAADAAMQZAADAAMQZAADAAMQZAEOoqrdU1R9Mj3+kqr5dVVsOse9pVdXTPrtW+Pzvq6o/r6p7VjHbf6uqX5ke/2JVfeLpPgcALEecATCc7v5adx/T3Y8nfzmOnuS47t79xEJVnVdVX6yq71TVTVX1vCXP+ctJXr0Os72/u39mrc8DAE8mzgDYFKrqOUmuTfJvkpyQZE+SD851KAB4GsQZAIdVVf3rqtpXVQeq6ktVdd5B9nnissWtVXV5kr+Z5Pemyxh/7xBP/feT3N7d/6W7v5vkLUnOqqoXrXCuI6vqD6rqT6vqoar6TFWddJD9frmqPr1k+cer6oaq+mZV3V9Vb57W/0BVXVZVX56e80NVdcJKZgFgMYkzAA6bqnphkkuS/PXuPjbJzya5+6mO6e7fTPLHSS6ZLnW85BC7/niSP1ly3J8l+fK0fiUuTvJDSU5N8uwkv5bkz5/qgKo6Nsknk3w8yXOTnJ7kxmnzP01yYZK/PW17MMm7VzgLAAtInAFwOD2eZHuSM6tqW3ff3d1fXqfnPibJt5607ltJjl3h8Y9mFmWnd/fj3X1Ldz+8zDGvTXJfd7+ju7/b3Qe6++Zp268l+c3uvqe7H8nsk7zXVdXWFc4DwIIRZwAcNt19V5JLMwuV/VV1dVU9d52e/ttJnvWkdc9KcmCFx//nJNcnubqq/m9V/Yeq2rbMMadm9uncwTwvyUemSyQfSrI3szj9/y6VBIBEnAFwmHX3H3b3yzOLl07y9pUctoJ9bk9y1hMLVXV0khdM61cy16Pd/dbuPjPJ38jsU7F/vMxhX0/yo0+x7dXdfdySnyO7e99K5gFg8YgzAA6bqnphVZ1bVduTfDez73R9fwWH3p9DR9ATPpLkJ6rq56vqyCT/Nsnnu/uLK5ztlVX1V6d/W+3hzC5zXG62jyXZUVWXVtX2qjq2ql4ybfuPSS5/4nb+VXViVV2wklkAWEziDIDDaXuStyX5RpL7kvxwkt9YwXG/k9n3tR6sqt892A7d/UCSn09yeWY333hJktc/jdn+SpJrMguzvUn+e2aXOh5Sdx9I8qokfzez/547k7xyyczXJflEVR1I8j+nmQDgoKp7JVeKAMA4pk+jvpTZp2//srvfu4JjrkjyD5Ls7+7TN3hEAHjaxBkAAMAAXNYIAAAwAHEGAAAwgMP6D2EeUdv7yBx9OF8SAABgGAfy4De6+8SDbTuscXZkjs5L6rzD+ZIAAADD+GRf89VDbXNZIwAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwAAO6630f+wnv5Prr7/1cL4kAADAMLbsOPQ2n5wBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMQJwBAAAMYNk4q6ojq+p/VdWfVNXtVfXWaf3zq+rmqrqrqj5YVUds/LgAAACb00o+OXskybndfVaSs5OcX1UvTfL2JO/q7tOTPJjkDRs2JQAAwCa3bJz1zLenxW3TTyc5N8k10/qrkly4EQMCAAAsghV956yqtlTVrUn2J7khyZeTPNTdj0273JPk5A2ZEAAAYAGsKM66+/HuPjvJKUnOSfKilb5AVe2qqj1VteeBP318dVMCAABsck/rbo3d/VCSm5L8VJLjqmrrtOmUJPsOcczu7t7Z3TtPfPaWtcwKAACwaa3kbo0nVtVx0+MfTPKqJHszi7TXTbtdnOSjGzQjAADAprd1+V2yI8lVVbUls5j7UHd/rKruSHJ1Vf37JJ9LcsUGzgkAALCpVXcfthd7Vp3QL6nzDtvrAQAAjOSTfc0t3b3zYNue1nfOAAAA2BjiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADiDAAAYADLxllVnVpVN1XVHVV1e1W9cVr/lqraV1W3Tj+v2fhxAQAANqetK9jnsSRv6u7PVtWxSW6pqhumbe/q7t/auPEAAAAWw7Jx1t33Jrl3enygqvYmOXmjBwMAAFgkT+s7Z1V1WpIXJ7l5WnVJVX2+qq6squMPccyuqtpTVXsezSNrmxYAAGCTWnGcVdUxST6c5NLufjjJe5K8IMnZmX2y9o6DHdfdu7t7Z3fv3Jbta58YAABgE1pRnFXVtszC7P3dfW2SdPf93f14d38/yXuTnLNxYwIAAGxuK7lbYyW5Isne7n7nkvU7luz2c0luW//xAAAAFsNK7tb4siS/lOQLVXXrtO7NSS6qqrOTdJK7k/zqck/0Yz/5nVx//a3L7QYAALApbdlx6G0ruVvjp5PUQTb90epHAgAAYKmndbdGAAAANoY4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGIA4AwAAGMDWw/li//vzR+Vnn3v24XxJAACAgdx1yC0+OQMAABiAOAMAABiAOAMAABiAOAMAABiAOAMAABiAOAMAABiAOAMAABiAOAMAABjAsnFWVadW1U1VdUdV3V5Vb5zWn1BVN1TVndPv4zd+XAAAgM1pJZ+cPZbkTd19ZpKXJvn1qjozyWVJbuzuM5LcOC0DAACwCsvGWXff292fnR4fSLI3yclJLkhy1bTbVUku3KAZAQAANr2tT2fnqjotyYuT3JzkpO6+d9p0X5KTDnHMriS7kuTIHLXqQQEAADazFd8QpKqOSfLhJJd298NLt3V3J+mDHdfdu7t7Z3fv3JbtaxoWAABgs1pRnFXVtszC7P3dfe20+v6q2jFt35Fk/8aMCAAAsPmt5G6NleSKJHu7+51LNl2X5OLp8cVJPrr+4wEAACyGlXzn7GVJfinJF6rq1mndm5O8LcmHquoNSb6a5Bc2ZEIAAIAFsGycdfenk9QhNp+3vuMAAAAsphXfEAQAAICNI84AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGIM4AAAAGsGycVdWVVbW/qm5bsu4tVbWvqm6dfl6zsWMCAABsbiv55Ox9Sc4/yPp3dffZ088fre9YAAAAi2XZOOvuTyX55mGYBQAAYGGt5Ttnl1TV56fLHo8/1E5Vtauq9lTVnkfzyBpeDgAAYPNabZy9J8kLkpyd5N4k7zjUjt29u7t3dvfObdm+ypcDAADY3FYVZ919f3c/3t3fT/LeJOes71gAAACLZVVxVlU7liz+XJLbDrUvAAAAy9u63A5V9YEkr0jynKq6J8m/S/KKqjo7SSe5O8mvbtyIAAAAm9+ycdbdFx1k9RUbMAsAAMDCWsvdGgEAAFgn4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAA4gwAAGAAy8ZZVV1ZVfur6rYl606oqhuq6s7p9/EbOyYAAMDmtpJPzt6X5PwnrbssyY3dfUaSG6dlAAAAVmnZOOvuTyX55pNWX5DkqunxVUkuXN+xAAAAFsvWVR53UnffOz2+L8lJh9qxqnYl2ZUkR+aoVb4cAADA5rbmG4J0dyfpp9i+u7t3dvfObdm+1pcDAADYlFYbZ/dX1Y4kmX7vX7+RAAAAFs9q4+y6JBdPjy9O8tH1GQcAAGAxreRW+h9I8j+SvLCq7qmqNyR5W5JXVdWdSX56WgYAAGCVlr0hSHdfdIhN563zLAAAAAtrzTcEAQAAYO3EGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwADEGQAAwAC2ruXgqro7yYEkjyd5rLt3rsdQAAAAi2ZNcTZ5ZXd/Yx2eBwAAYGG5rBEAAGAAa42zTvKJqrqlqnYdbIeq2lVVe6pqz6N5ZI0vBwAAsDmt9bLGl3f3vqr64SQ3VNUXu/tTS3fo7t1JdifJs+qEXuPrAQAAbEpr+uSsu/dNv/cn+UiSc9ZjKAAAgEWz6jirqqOr6tgnHif5mSS3rddgAAAAi2QtlzWelOQjVfXE8/xhd398XaYCAABYMKuOs+7+SpKz1nEWAACAheVW+gAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAAMQZwAAAANYU5xV1flV9aWququqLluvoQAAABbNquOsqrYkeXeSVyc5M8lFVXXmeg0GAACwSNbyydk5Se7q7q909/eSXJ3kgvUZCwAAYLGsJc5OTvL1Jcv3TOv+kqraVVV7qmrPo3lkDS8HAACweW34DUG6e3d37+zunduyfaNfDgAA4BlpLXG2L8mpS5ZPmdYBAADwNK0lzj6T5Iyqen5VHZHk9UmuW5+xAAAAFkt19+oPrnpNkt9OsiXJld19+TL7P5Dkq0tWPSfJN1Y9AGws709G5b3JqLw3GZn3J6N4XnefeLANa4qztaqqPd29c24DwFPw/mRU3puMynuTkXl/8kyw4TcEAQAAYHniDAAAYADzjrPdc359eCren4zKe5NReW8yMu9PhjfX75wBAAAwM+9PzgAAAIg4AwAAGMLc4qyqzq+qL1XVXVV12bzmgKo6tapuqqo7qur2qnrjtP6Eqrqhqu6cfh8/71lZTFW1pao+V1Ufm5afX1U3T+fPD1bVEfOekcVUVcdV1TVV9cWq2ltVP+XcyQiq6p9Pf6bfVlUfqKojnTt5JphLnFXVliTvTvLqJGcmuaiqzpzHLJDksSRv6u4zk7w0ya9P78fLktzY3WckuXFahnl4Y5K9S5bfnuRd3X16kgeTvGEuU0HyO0k+3t0vSnJWZu9T507mqqpOTvLPkuzs7p9IsiXJ6+PcyTPAvD45OyfJXd39le7+XpKrk1wwp1lYcN19b3d/dnp8ILP/uTg5s/fkVdNuVyW5cC4DstCq6pQkfyfJ70/LleTcJNdMu3hvMhdV9UNJ/laSK5Kku7/X3Q/FuZMxbE3yg1W1NclRSe6NcyfPAPOKs5OTfH3J8j3TOpirqjotyYuT3JzkpO6+d9p0X5KT5jUXC+23k/yrJN+flp+d5KHufmxadv5kXp6f5IEk/2m67Pb3q+roOHcyZ929L8lvJflaZlH2rSS3xLmTZwA3BIFJVR2T5MNJLu3uh5du69m/OeHfneCwqqrXJtnf3bfMexY4iK1J/lqS93T3i5P8WZ50CaNzJ/Mwfc/xgsz+AuG5SY5Ocv5ch4IVmlec7Uty6pLlU6Z1MBdVtS2zMHt/d187rb6/qnZM23ck2T+v+VhYL0vy96rq7swu/z43s+/4HDddqpM4fzI/9yS5p7tvnpavySzWnDuZt59O8n+6+4HufjTJtZmdT507Gd684uwzSc6Y7ppzRGZf0rxuTrOw4Kbv8FyRZG93v3PJpuuSXDw9vjjJRw/3bCy27v6N7j6lu0/L7Dz5X7v7F5PclOR1027em8xFd9+X5OtV9cJp1XlJ7ohzJ/P3tSQvraqjpj/jn3hvOncyvJpdcTCHF656TWbfpdiS5Mruvnwug7DwqurlSf44yRfyF9/reXNm3zv7UJIfSfLVJL/Q3d+cy5AsvKp6RZJ/0d2vraofzeyTtBOSfC7JP+ruR+Y4Hguqqs7O7GY1RyT5SpJ/ktlf/Dp3MldV9dYk/zCzOzJ/LsmvZPYdM+dOhja3OAMAAOAvuCEIAADAAMQZAADAAMQZAADAAMQZAADAAMQZAADAAMQZAADAAMQZAADAAP4fayc9iLtgdBgAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "bg_med = Background.two_sided(image, ext_center, bkg_sep, width=bkg_width, statistic='median')\n", + "plt.figure(figsize=(15, 15))\n", + "plt.imshow(bg_med.bkg_wimage[::,0:100], origin=\"lower\")\n", + "plt.title(\"slit[0] slice\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Trace" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, "outputs": [], "source": [ "# optional: refine the trace on the background subtracted image (once PR#85 is merged)\n", @@ -530,7 +572,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ @@ -541,7 +583,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 20, "metadata": {}, "outputs": [ { diff --git a/specreduce/background.py b/specreduce/background.py index e18bc181..da6496df 100644 --- a/specreduce/background.py +++ b/specreduce/background.py @@ -25,6 +25,10 @@ class Background: extract the background width : float width of extraction aperture in pixels + statistic: string + statistic to use when computing the background. 'average' will + account for partial pixel weights, 'median' will include all partial + pixels. disp_axis : int dispersion axis crossdisp_axis : int @@ -37,6 +41,7 @@ class Background: image: CCDData traces: list = field(default_factory=list) width: float = 5 + statistic: str = 'average' disp_axis: int = 1 crossdisp_axis: int = 0 @@ -53,6 +58,10 @@ def __post_init__(self): extract the background width : float width of each background aperture in pixels + statistic: string + statistic to use when computing the background. 'average' will + account for partial pixel weights, 'median' will include all partial + pixels. disp_axis : int dispersion axis crossdisp_axis : int @@ -80,8 +89,17 @@ def _to_trace(trace): if np.any(bkg_wimage > 1): raise ValueError("background regions overlapped") + if self.statistic == 'median': + # make it clear in the expose image that partial pixels are fully-weighted + bkg_wimage[bkg_wimage > 0] = 1 + self.bkg_wimage = bkg_wimage - self.bkg_array = np.average(self.image, weights=self.bkg_wimage, axis=0) + if self.statistic == 'average': + self.bkg_array = np.average(self.image, weights=self.bkg_wimage, axis=0) + elif self.statistic == 'median': + self.bkg_array = np.median(self.image[self.bkg_wimage > 0], axis=0) + else: + raise ValueError("statistic must be 'average' or 'median'") @classmethod def two_sided(cls, image, trace_object, separation, **kwargs): @@ -99,6 +117,10 @@ def two_sided(cls, image, trace_object, separation, **kwargs): separation from ``trace_object`` for the background regions width : float width of each background aperture in pixels + statistic: string + statistic to use when computing the background. 'average' will + account for partial pixel weights, 'median' will include all partial + pixels. disp_axis : int dispersion axis crossdisp_axis : int @@ -124,6 +146,10 @@ def one_sided(cls, image, trace_object, separation, **kwargs): above the trace, negative below. width : float width of each background aperture in pixels + statistic: string + statistic to use when computing the background. 'average' will + account for partial pixel weights, 'median' will include all partial + pixels. disp_axis : int dispersion axis crossdisp_axis : int