support statistic='median'

* defaults to 'average' * median does not support partial-pixel-weights
astropy · Mar 17, 2022 · d91df4e · d91df4e
1 parent ba994c3
commit d91df4e
Show file tree

Hide file tree

Showing 2 changed files with 72 additions and 4 deletions.
diff --git 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": [
+       "<Figure size 1080x1080 with 1 Axes>"
+      ]
+     },
+     "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
@@ -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