{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get image from the PS1 image server\n",
    "Query the PS1 image server to get a list of images and retrieve some images.  This sample script demonstrates the use of the PS1 image services.  See the <a href=\"https://outerspace.stsci.edu/x/ioOc\">PS1 Image Cutout Service documentation</a> for details of the services being used.  This notebook is available for \n",
    "<a href=\"https://ps1images.stsci.edu/ps1image.ipynb\" download=\"\">download</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from astropy.table import Table\n",
    "import requests\n",
    "from PIL import Image\n",
    "from io import BytesIO\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper functions to query the list of images and to extract images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getimages(ra,dec,filters=\"grizy\"):\n",
    "    \n",
    "    \"\"\"Query ps1filenames.py service to get a list of images\n",
    "    \n",
    "    ra, dec = position in degrees\n",
    "    size = image size in pixels (0.25 arcsec/pixel)\n",
    "    filters = string with filters to include\n",
    "    Returns a table with the results\n",
    "    \"\"\"\n",
    "    \n",
    "    service = \"https://ps1images.stsci.edu/cgi-bin/ps1filenames.py\"\n",
    "    url = f\"{service}?ra={ra}&dec={dec}&filters={filters}\"\n",
    "    table = Table.read(url, format='ascii')\n",
    "    return table\n",
    "\n",
    "\n",
    "def geturl(ra, dec, size=240, output_size=None, filters=\"grizy\", format=\"jpg\", color=False):\n",
    "    \n",
    "    \"\"\"Get URL for images in the table\n",
    "    \n",
    "    ra, dec = position in degrees\n",
    "    size = extracted image size in pixels (0.25 arcsec/pixel)\n",
    "    output_size = output (display) image size in pixels (default = size).\n",
    "                  output_size has no effect for fits format images.\n",
    "    filters = string with filters to include\n",
    "    format = data format (options are \"jpg\", \"png\" or \"fits\")\n",
    "    color = if True, creates a color image (only for jpg or png format).\n",
    "            Default is return a list of URLs for single-filter grayscale images.\n",
    "    Returns a string with the URL\n",
    "    \"\"\"\n",
    "    \n",
    "    if color and format == \"fits\":\n",
    "        raise ValueError(\"color images are available only for jpg or png formats\")\n",
    "    if format not in (\"jpg\",\"png\",\"fits\"):\n",
    "        raise ValueError(\"format must be one of jpg, png, fits\")\n",
    "    table = getimages(ra,dec,filters=filters)\n",
    "    url = (f\"https://ps1images.stsci.edu/cgi-bin/fitscut.cgi?\"\n",
    "           f\"ra={ra}&dec={dec}&size={size}&format={format}\")\n",
    "    if output_size:\n",
    "        url = url + \"&output_size={}\".format(output_size)\n",
    "    # sort filters from red to blue\n",
    "    flist = [\"yzirg\".find(x) for x in table['filter']]\n",
    "    table = table[numpy.argsort(flist)]\n",
    "    if color:\n",
    "        if len(table) > 3:\n",
    "            # pick 3 filters\n",
    "            table = table[[0,len(table)//2,len(table)-1]]\n",
    "        for i, param in enumerate([\"red\",\"green\",\"blue\"]):\n",
    "            url = url + \"&{}={}\".format(param,table['filename'][i])\n",
    "    else:\n",
    "        urlbase = url + \"&red=\"\n",
    "        url = []\n",
    "        for filename in table['filename']:\n",
    "            url.append(urlbase+filename)\n",
    "    return url\n",
    "\n",
    "\n",
    "def getcolorim(ra, dec, size=240, output_size=None, filters=\"grizy\", format=\"jpg\"):\n",
    "    \n",
    "    \"\"\"Get color image at a sky position\n",
    "    \n",
    "    ra, dec = position in degrees\n",
    "    size = extracted image size in pixels (0.25 arcsec/pixel)\n",
    "    output_size = output (display) image size in pixels (default = size).\n",
    "                  output_size has no effect for fits format images.\n",
    "    filters = string with filters to include\n",
    "    format = data format (options are \"jpg\", \"png\")\n",
    "    Returns the image\n",
    "    \"\"\"\n",
    "    \n",
    "    if format not in (\"jpg\",\"png\"):\n",
    "        raise ValueError(\"format must be jpg or png\")\n",
    "    url = geturl(ra,dec,size=size,filters=filters,output_size=output_size,format=format,color=True)\n",
    "    r = requests.get(url)\n",
    "    im = Image.open(BytesIO(r.content))\n",
    "    return im\n",
    "\n",
    "\n",
    "def getgrayim(ra, dec, size=240, output_size=None, filter=\"g\", format=\"jpg\"):\n",
    "    \n",
    "    \"\"\"Get grayscale image at a sky position\n",
    "    \n",
    "    ra, dec = position in degrees\n",
    "    size = extracted image size in pixels (0.25 arcsec/pixel)\n",
    "    output_size = output (display) image size in pixels (default = size).\n",
    "                  output_size has no effect for fits format images.\n",
    "    filter = string with filter to extract (one of grizy)\n",
    "    format = data format (options are \"jpg\", \"png\")\n",
    "    Returns the image\n",
    "    \"\"\"\n",
    "    \n",
    "    if format not in (\"jpg\",\"png\"):\n",
    "        raise ValueError(\"format must be jpg or png\")\n",
    "    if filter not in list(\"grizy\"):\n",
    "        raise ValueError(\"filter must be one of grizy\")\n",
    "    url = geturl(ra,dec,size=size,filters=filter,output_size=output_size,format=format)\n",
    "    r = requests.get(url[0])\n",
    "    im = Image.open(BytesIO(r.content))\n",
    "    return im"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample JPEG images\n",
    "This gets single-band grayscale and color JPEG images at the position of the Crab Nebula.  The extracted region size is 1280 pixels = 320 arcsec."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crab Nebula position\n",
    "ra = 83.633210\n",
    "dec = 22.014460\n",
    "size = 1280\n",
    "\n",
    "# grayscale image\n",
    "gim = getgrayim(ra,dec,size=size,filter=\"i\")\n",
    "# color image\n",
    "cim = getcolorim(ra,dec,size=size,filters=\"grz\")\n",
    "\n",
    "plt.rcParams.update({'font.size':12})\n",
    "plt.figure(1,(12,6))\n",
    "plt.subplot(121)\n",
    "plt.imshow(gim,origin=\"upper\",cmap=\"gray\")\n",
    "plt.title('Crab Nebula PS1 i')\n",
    "plt.subplot(122)\n",
    "plt.title('Crab Nebula PS1 grz')\n",
    "plt.imshow(cim,origin=\"upper\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and display a FITS image\n",
    "\n",
    "Note that the $y$-axis is flipped in the JPEG image compared with the original FITS image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy.visualization import PercentileInterval, AsinhStretch\n",
    "\n",
    "fitsurl = geturl(ra, dec, size=size, filters=\"i\", format=\"fits\")\n",
    "fh = fits.open(fitsurl[0])\n",
    "fim = fh[0].data\n",
    "# replace NaN values with zero for display\n",
    "fim[numpy.isnan(fim)] = 0.0\n",
    "# set contrast to something reasonable\n",
    "transform = AsinhStretch() + PercentileInterval(99.5)\n",
    "bfim = transform(fim)\n",
    "\n",
    "plt.figure(1,(12,6))\n",
    "plt.subplot(121)\n",
    "plt.imshow(gim,cmap=\"gray\",origin=\"upper\")\n",
    "plt.title('Crab Nebula PS1 i (jpeg)')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.title('Crab Nebula PS1 i (fits)')\n",
    "plt.imshow(bfim,cmap=\"gray\",origin=\"lower\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}