Lec011_astropy-example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a0106a1-fcaa-4761-8da4-92a4bee4b4bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# \n",
    "# Notebook illustrating a few astropy ( https://www.astropy.org ) examples.\n",
    "# Although astropy was designed for the astronomical community many of its\n",
    "# builtin functions are general and can be used for a wide array of \n",
    "# image and time series processing for all sorts of data.\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11d66b0e-cc00-4052-bae4-cfe5ee833bd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Functions like Gauss2d (see below) provide simple ways to generate \n",
    "# Gaussian distributions with various parameters. \n",
    "#\n",
    "import numpy as np\n",
    "from astropy.modeling.models import Gaussian2D\n",
    "y, x = np.mgrid[0:500, 0:600]\n",
    "data = (Gaussian2D(1, 150, 100, 20, 10, theta=0.5)(x, y) +\n",
    "        Gaussian2D(0.5, 400, 300, 8, 12, theta=1.2)(x,y) +\n",
    "        Gaussian2D(0.75, 250, 400, 5, 7, theta=0.23)(x,y) +\n",
    "        Gaussian2D(0.9, 525, 150, 3, 3)(x,y) +\n",
    "        Gaussian2D(0.6, 200, 225, 3, 3)(x,y))\n",
    "data += 0.01 * np.random.randn(500, 600)\n",
    "#\n",
    "# This is an example of some \"noise\"!\n",
    "#\n",
    "cosmic_ray_value = 0.997\n",
    "data[100, 300:310] = cosmic_ray_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc8a8af8-ac58-4d1a-bfd1-ef6637b833c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.imshow(data, origin='lower');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64e24a17-385c-4cfd-93b3-b32b94d345c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.nddata import CCDData\n",
    "ccd = CCDData(data, \n",
    "              meta={'object': 'fake galaxy', 'filter': 'R'},\n",
    "              unit='adu')\n",
    "\n",
    "from astropy.nddata import Cutout2D\n",
    "position = (149.7, 100.1)\n",
    "size = (81, 101)     # pixels\n",
    "cutout = Cutout2D(ccd, position, size)\n",
    "plt.imshow(cutout.data, origin='lower') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80d87f99-e786-4918-abde-738e44d7d12e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(ccd, origin='lower')  \n",
    "cutout.plot_on_original(color='white') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63ccb805-b541-4ca4-89a3-7f81be7b37da",
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils.segmentation import detect_threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd9df95-729a-457c-9817-e8569c705307",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = detect_threshold(ccd.data, nsigma=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd0c9ec-7c65-42cf-ad88-ee70661140a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.convolution import Gaussian2DKernel\n",
    "from astropy.stats import gaussian_fwhm_to_sigma\n",
    "from photutils.segmentation import detect_sources\n",
    "sigma = 3.0 * gaussian_fwhm_to_sigma  # FWHM = 3.\n",
    "kernel = Gaussian2DKernel(sigma, x_size=3, y_size=3)\n",
    "kernel.normalize()\n",
    "segm = detect_sources(ccd.data, threshold, npixels=5, kernel=kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21ba9a94-678f-468a-9bba-32d0a1d10b66",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.visualization import SqrtStretch\n",
    "from astropy.visualization.mpl_normalize import ImageNormalize\n",
    "norm = ImageNormalize(stretch=SqrtStretch())\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12.5))\n",
    "ax1.imshow(ccd.data, origin='lower', cmap='Greys_r', norm=norm)\n",
    "ax1.set_title('Data')\n",
    "cmap = segm.make_cmap(seed=123)\n",
    "ax2.imshow(segm, origin='lower', cmap=cmap, interpolation='nearest')\n",
    "ax2.set_title('Segmentation Image')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4f16619-0ace-4542-b67b-ed8edeab768b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Scipy Curve fit equaivalent\n",
    "#\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.modeling import models, fitting\n",
    "\n",
    "\n",
    "mu, sigma, amplitude = 0.0, 10.0, 10.0\n",
    "N2 = 100\n",
    "x2 = np.linspace(-30, 30, N2)\n",
    "y2 = amplitude * np.exp(-(x2-mu)**2 / (2*sigma**2))\n",
    "y2 = np.array([y_point + np.random.normal(0, 1) for y_point in y2])   #Another way to add random gaussian noise\n",
    "sigma = 1\n",
    "y2_err = np.ones(N2)*sigma\n",
    "plt.errorbar(x2, y2, yerr=y2_err, fmt='k.')\n",
    "plt.xlabel('$x_2$')\n",
    "plt.ylabel('$y_2$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20687e7e-7396-4ffc-87ab-0182215cf74d",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_gauss = models.Gaussian1D()\n",
    "fitter_gauss = fitting.LevMarLSQFitter()\n",
    "best_fit_gauss = fitter_gauss(model_gauss, x2, y2, weights=1/y2_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ad79fb-48a1-4c70-9a67-0ce4853b9c7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x2, y2, yerr=y2_err, fmt='k.')\n",
    "plt.plot(x2, best_fit_gauss(x2), 'g-', linewidth=6, label='astropy.modeling')\n",
    "plt.xlabel('$x_2$')\n",
    "plt.ylabel('$y_2$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f00ad287-a628-49ae-86ae-a55eec86480d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "import astropy\n",
    "?astropy.modeling.models.Legendre1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b03606b-aa2f-41af-bed2-2fd8f7f52b6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "M1=astropy.modeling.models.Legendre1D(1,c1=1)\n",
    "M2=astropy.modeling.models.Legendre1D(2,c2=1)\n",
    "M3=astropy.modeling.models.Legendre1D(3,c3=1)\n",
    "M4=astropy.modeling.models.Legendre1D(4,c4=1)\n",
    "M5=astropy.modeling.models.Legendre1D(5,c5=1)\n",
    "M6=astropy.modeling.models.Legendre1D(6,c6=1)\n",
    "M7=astropy.modeling.models.Legendre1D(7,c7=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ff89501-7bf5-4273-9960-e96158573b6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.linspace(-1,1,200);\n",
    "plt.plot(x, M1(x) )\n",
    "plt.plot(x, M2(x) )\n",
    "plt.plot(x, M3(x) )\n",
    "plt.plot(x, M4(x) )\n",
    "plt.plot(x, M5(x) )\n",
    "plt.plot(x, M6(x) )\n",
    "plt.plot(x, M7(x) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acffc4aa-3ab1-434c-a004-8590927102a0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}