MeerKAT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img style=\"float: right;\" src=\"http://montage.ipac.caltech.edu/docs/apjac449af1_hr.jpeg\" width =600 height= 600 />\n",
    "\n",
    "\n",
    "# <font color=\"#880000\"> Creating a Merger of a MeerKAT map of the Galactic Center and a 2MASS mosaic created with the Montage Image Mosaic Engine\n",
    "\n",
    "### John Good, Rosemary Moseley and G. Bruce Berriman (Caltech/IPAC-NExScI)\n",
    "\n",
    "<p style=\"color:#c00000;\"> This Notebook demonstrates how to create a three-color Montage mosaic of images in the 2MASS all-sky survey and overlay a MeerKAT (1.28 GHz) map of the Galactic Center, itself a mosaic created with Montage (Heywood, I., Rammala, I., Camilo, F., et al.,  2022, \"The 1.28 GHz MeerKAT Galactic Center Mosaic\", The Ap J, 925, 165) and shown on the right. The Notebook will run on all all recent Python versions installed on *nix or Mac OS (Intel or Apple silicon chips). It was tested with Python 3.11; we recommend using only Python 3.8 or later.  Processing will require about 15 minutes on an Apple M2 chip with 24 GB memory (longer on a slow network connection). </p>\n",
    "\n",
    "Montage is a toolkit that assembles Flexible Image Transport System (FITS) images into science-grade custom mosaics that preserve the astrometric and calibration fidelity of the original images. See http://montage.ipac.caltech.edu/. Download the code from https://github.com/Caltech-IPAC/Montage. \n",
    "\n",
    "Read about the architecture in Jacob, Joseph C., Katz, Daniel S., Berriman, G. Bruce, et al.,  2010, \"Montage: a grid portal and software toolkit for science-grade astronomical image mosaicking\", https://doi.org/10.48550/arXiv.1005.4454. \n",
    "\n",
    "Read about Montage and image visualization in Berriman, G. Bruce & Good, J. C.,  2017, \"The Application of the Montage Image Mosaic Engine To The Visualization Of Astronomical Images\", Publications of the Astronomical Society of the Pacific, 129, 058006. https:/doi.org/10.1088/1538-3873/aa5456.\n",
    "\n",
    "\n",
    "## <font color=\"#880000\">Montage and Radio Astronomy\n",
    "\n",
    "In recent years, Montage has found growing applicability in radio astronomy. Here are some papers published in 2024 that cite Montage:\n",
    "\n",
    "* Rajohnson, Sambatriniaina H. A., Kraan-Korteweg, Renée C., Chen, Hao, et al.,  2024, \"H I galaxy signatures in the SARAO MeerKAT galactic plane survey - III. Unveiling the obscured part of the Vela Supercluster\", Monthly Notices of the Royal Astronomical Society, 531, 3486\n",
    "\n",
    "* Raycheva, N., Haverkorn, M., Ideguchi, S., et al.,  2024, \"Faraday moments of the Southern Twenty-centimeter All-sky Polarization Survey (STAPS)\", arXiv e-prints, arXiv:2406.06166\n",
    "\n",
    "* Goedhart, S., Cotton, W. D., Camilo, F., et al.,  2024, \"The SARAO MeerKAT 1.3 GHz Galactic Plane Survey\", Monthly Notices of the Royal Astronomical Society, 531, 649\n",
    "\n",
    "* Vanderwoude, S., West, J. L., Gaensler, B. M., et al.,  2024, \"Prototype Faraday Rotation Measure Catalogs from the Polarisation Sky Survey of the Universe's Magnetism (POSSUM) Pilot Observations\", The Astronomical Journal, 167, 226\n",
    "\n",
    "* Kurapati, Sushma, Kraan-Korteweg, Renée C., Pisano, D. J., et al.,  2024, \"H I galaxy signatures in the SARAO MeerKAT galactic plane survey - II. The Local Void and its substructure\", Monthly Notices of the Royal Astronomical Society, 528, 542\n",
    "\n",
    "* Riggi, S., Umana, G., Trigilio, C., et al.,  2024, \"Classification of compact radio sources in the Galactic plane with supervised machine learning\", Publications of the Astronomical Society of Australia, 41, e029\n",
    "\n",
    "* Khabibullin, I. I., Churazov, E. M., Bykov, A. M., et al.,  2024, \"Discovery of a one-sided radio filament of PSR J0538+2817 in S147: escape of relativistic PWN leptons into surrounding supernova remnant?\", Monthly Notices of the Royal Astronomical Society, 527, 5683\n",
    "\n",
    "* Rajohnson, Sambatriniaina H. A., Kraan-Korteweg, Renée C., Chen, Hao, et al.,  2024, \"H I galaxy signatures in the SARAO MeerKAT galactic plane survey - III. Unveiling the obscured part of the Vela Supercluster\", Monthly Notices of the Royal Astronomical Society, 531, 3486\n",
    "\n",
    "* Raycheva, N., Haverkorn, M., Ideguchi, S., et al.,  2024, \"Faraday moments of the Southern Twenty-centimeter All-sky Polarization Survey (STAPS)\", arXiv e-prints, arXiv:2406.06166\n",
    "\n",
    "* Goedhart, S., Cotton, W. D., Camilo, F., et al.,  2024, \"The SARAO MeerKAT 1.3 GHz Galactic Plane Survey\", Monthly Notices of the Royal Astronomical Society, 531, 649\n",
    "\n",
    "* Kurapati, Sushma, Kraan-Korteweg, Renée C., Pisano, D. J., et al.,  2024, \"H I galaxy signatures in the SARAO MeerKAT galactic plane survey - II. The Local Void and its substructure\", Monthly Notices of the Royal Astronomical Society, 528, 542\n",
    "\n",
    "## <font color=\"#880000\">Notebook Setup\n",
    "\n",
    "We recommend using the Anaconda distribution. All of our development and testing has been done using Python 3.11 but any of the current versions (3.8 to 3.12) will probably work.\n",
    "    \n",
    "You will need to install the following packages if you do not alread have them:\n",
    "    \n",
    "\n",
    "* OpenCV (<tt style=\"color: #c00000;\">pip install opencv-python</tt>).  An extensive toolkit of image manipulation tools (along the lines of Photoshop).\n",
    "    \n",
    "* build (<tt style=\"color: #c00000;\">pip install build</tt>).  A Python package building front-end used when building Montage as a Python extension.\n",
    "    \n",
    "* importlib_resources (<tt style=\"color: #c00000;\">pip install importlib_resources</tt>). A library which provides for access to resources in Python packages. \n",
    "   \n",
    "    \n",
    "You also need the Montage Python package (MontagePy) but this you need to build from scratch since the \"pip installable\" version doesn't include some of the newest capabilities.  The Montage Python package is essentially a Python binary extension built from Montage itself (code written in C; http://montage.ipac,caltech.edu). Contains all tools needed to compute a mosaic of images written in FITS format.\n",
    "    \n",
    "We recommend you build the wheel locally for this tutorial, to avoid problems with supporting different Python versions through PyPI. The process will take less than 5 minutes:\n",
    "    \n",
    "\n",
    "* Download the Montage code from GitHub (\"git clone https://github.com/Caltech-IPAC/Montage.git\") \n",
    "      \n",
    "* Make Montage (\"cd Montage; make\").\n",
    "\n",
    "* Make the the wheel file for the host machine (\"cd python/MontagePy; ./make_local.sh\").\n",
    "\n",
    "* Install the wheel under Python (\"cd dist; pip install --force-reinstall *.whl\").\n",
    "\n",
    "* Download the Montage Notebooks, as this is where this MeerKAT.ipynb notebook resides (\"git clone https://github.com/Caltech-IPAC/MontageNotebooks.git\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Startup.  The Montage modules are pretty much self-contained\n",
    "# but this script needs a few extra utilities.  If you don't have\n",
    "# these, they also can easily be \"pip installed\".\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import shutil\n",
    "import cv2 as cv\n",
    "import importlib_resources\n",
    "\n",
    "from MontagePy.main    import *\n",
    "from MontagePy.archive import *\n",
    "\n",
    "from IPython.display import Image\n",
    "\n",
    "# Capture and print the current date and time for logging purposes\n",
    "import datetime\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "\n",
    "# Note: we create and move into several subdirectories in this notebook \n",
    "# but we want to come back to the original startup directory \n",
    "# whenever we restart the processing. \n",
    "\n",
    "# Here we use the directory where we started the Notebook as that\n",
    "# home directory.  Feel free to explicitly point somewhere else.\n",
    "\n",
    "try:\n",
    "    home\n",
    "except:\n",
    "    home = os.getcwd() # Get the current working directory\n",
    "\n",
    "# Change to the home directory\n",
    "os.chdir(home)\n",
    "\n",
    "print(\"Startup folder: \" + home)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Downloading the MeerKAT Data\n",
    "\n",
    "Download the MeerKAT mosaic as a starting point:\n",
    "\n",
    "https://archive-gw-1.kat.ac.za/public/repository/10.48479/fyst-hj47/data/MeerKAT_Galactic_Centre_1284MHz-StokesI.fits\n",
    "\n",
    "Or visit:\n",
    "https://archive-gw-1.kat.ac.za/public/repository/10.48479/fyst-hj47/index.html \n",
    "\n",
    "and download:\n",
    "\n",
    "`MeerKAT_Galactic_Centre_1284MHz-StokesI.fits`\n",
    "\n",
    "\n",
    "This file is nearly 1 GB, so we recommend downloading it once and saving it to the directory designated as 'home' above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Steps in Computing the 2MASS Mosaic and Displaying the MeerKAT Overlay ..\n",
    "\n",
    "Our goal is to build a mosaic of the 2MASS data that we can overlay on the MeerKAT image. The steps are:\n",
    "\n",
    "* Define a FITS header for the output 2MASS mosaic file.\n",
    "* Set up a working environment for the processing.\n",
    "* Retrieving 2MASS data.\n",
    "* Reproject all the images to the same scale and projection as MeerKAT.\n",
    "* Coadd the reprojected images.\n",
    "* Visualize the image.\n",
    "* As an additional step, rectify the sky background in the 2MASS images to a common level, and create a new background corrected mosaic.\n",
    " \n",
    "Repeat these steps for all three 2MASS wavelengths and build a color composite PNG of the three.  Merge with the MeerKAT image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Define a FITS Header for the Output Mosaic File.\n",
    "\n",
    "Since we want our 2MASS mosaics to match the MeerKAT image, the simplest thing is to use the MeerKAT FITS header as the template for building them.  mGetHdr extracts the header from an image and writes it to a text file.  The MeerKAT header does contain keywords specific to MeerKAT, but we will only make use of the  projection information:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the FITS header from the MeerKAT image:\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "# Ensure we are working from the home directory\n",
    "os.chdir(home)\n",
    "\n",
    "\n",
    "# Define the path to the MeerKAT FITS file\n",
    "meerkat = home + \"/MeerKAT_Galactic_Centre_1284MHz-StokesI.fits\"\n",
    "\n",
    "# Use the Montage utility 'mGetHdr' to extract the FITS header and save it to 'region.hdr'\n",
    "rtn = mGetHdr(meerkat, \"region.hdr\")\n",
    "\n",
    "# Print the result of the mGetHdr operation for verification\n",
    "print(\"mGetHdr:          \" + str(rtn), flush=True)\n",
    "\n",
    "# Use the Montage utility 'mExamine' to extract key information about the FITS image\n",
    "rtn = mExamine(meerkat)\n",
    "\n",
    "# Extract the center location (longitude and latitude) of the image\n",
    "location = str(rtn[\"lonc\"]) + \" \" + str(rtn[\"latc\"])\n",
    "\n",
    "# Get the size of the image in pixels\n",
    "size = rtn[\"ximgsize\"]\n",
    "if rtn[\"yimgsize\"] > size:\n",
    "    size = rtn[\"yimgsize\"] # Set 'size' to the larger dimension (x or y)\n",
    "    \n",
    "# Calculate a ratio to scale the image to a standard size (1000 pixels)\n",
    "ratio = 1000./size\n",
    "\n",
    "# Calculate the new width and height based on the scaling ratio\n",
    "width  = rtn[\"ximgsize\"] * ratio\n",
    "height = rtn[\"yimgsize\"] * ratio\n",
    "\n",
    "# Define the shape as a tuple with the new dimensions (width, height)\n",
    "shape = (int(width), int(height))\n",
    "\n",
    "# Print the extracted and calculated details: location, size, and shape of the image\n",
    "print(\"location:         \" + location, flush=True)\n",
    "print(\"size:             \" + str(size))\n",
    "print(\"shape:            \" + str(shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <font color=\"#880000\">Set up a Working Environment for the Processoing\n",
    "\n",
    "We need to set up our working environment for building the mosaic.  We need to set up some subdirectories to manage the data.  This will all be under an instance-specific directory specified below (\"workdir\").  It is best not to use directory names with embedded spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean out any old copy of the work tree, then remake it\n",
    "# and the set of the subdirectories we will need.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "# Define the name of the working directory\n",
    "workdir = \"2MASSJ\"\n",
    "\n",
    "# Ensure we are working from the home directory\n",
    "os.chdir(home)\n",
    "\n",
    "# Attempt to delete the existing work directory, if it exists\n",
    "try:\n",
    "    shutil.rmtree(workdir)\n",
    "except:\n",
    "    print(\"Can't delete work tree; probably doesn't exist yet.\\n\", flush=True)\n",
    "\n",
    "# Print the name of the work directory for verification\n",
    "print(\"Work directory: \" + workdir, flush=True)\n",
    "\n",
    "# Create the main work directory where all subsequent processing will occur\n",
    "os.makedirs(workdir)  \n",
    "\n",
    "# Change into the newly created work directory\n",
    "os.chdir(workdir)\n",
    "\n",
    "# Create subdirectories within the work directory\n",
    "os.makedirs(\"raw\")\n",
    "os.makedirs(\"projected\")\n",
    "os.makedirs(\"diffs\")\n",
    "os.makedirs(\"corrected\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Retrieving Data from an Archive\n",
    "\n",
    "mArchiveDownload provides access to four wide-area sky surveys: 2MASS, WISE, SDSS and DSS.  Under the hood, mArchiveDownload uses mSearch to search for image metadata , along with a very fast R-Tree / memory-mapped utility. It returns URLS for all files within the search area. Once we have the data, scan the images with mImgtbl to retrieve their WCS sky coverage metadata\n",
    "\n",
    "Note 1: There are other ways to find images.  e.g. The International Virtual Astronomy Alliance (IVOA) has developed astandards for querying metadata (Simple Image Access: SIAP and Table Access Protocol: TAP) which many data providers support. \n",
    "\n",
    "Note 2: Other datasets may require more care.  Downloading all pointed observations of a specific region for a non-survey instrument will, for example, include a wide range of integration times and therefore noise levels.  The mosaicking should involve user-specified weighting of the images: Montage supports this but does not define weightings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve archive images covering the region then scan \n",
    "# the images for their coverage metadata.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "# Change to the working directory where the images will be stored and processed\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "dataset = '2MASS J'\n",
    "\n",
    "# Use the Montage utility 'mArchiveDownload' to retrieve images from the specified dataset\n",
    "# covering the region defined by 'location' and 'size'. The images will be stored in the 'raw' directory.\n",
    "rtn = mArchiveDownload(dataset, location, size, \"raw\")\n",
    "\n",
    "# Print the result of the mArchiveDownload operation for verification\n",
    "print(\"mArchiveDownload: \" + str(rtn), flush=True)\n",
    "\n",
    "# Use the Montage utility 'mImgtbl' to scan the downloaded images in the 'raw' directory\n",
    "# and generate a metadata table ('rimages.tbl') that describes the coverage of each image.\n",
    "rtn = mImgtbl(\"raw\", \"rimages.tbl\")\n",
    "\n",
    "print(\"mImgtbl (raw):    \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Reproject the Images\n",
    "\n",
    "We now have all the information to reproject all the images to a common frame.\n",
    "\n",
    "Montage has four reprojection modules:\n",
    "\n",
    "* mProject    - general purpose, based on spherical trigonometry and applies to all WCS projections, but can be slow. Conserves flux\n",
    "* mProjectPP  - fast custom plane-to-plane algorithm developed by the Spitzer Space Telescope and applies to small images with tangent plane projections. Conserves flux.\n",
    "* mProjectCube - reprojects three- and four-dimensional datacubes. Conserves flux.\n",
    "* mProjectQL - Uses Lanczos interpolation for speed: it is 20x times fadter Not flux conserving by itself. Montage does however take advantage of distortion parameters in FITS headers,  to mimic an arbitrary projection over a small region with a distorted gnomonic projection.  This allows us to use mProjectPP over a wider range of cases and still have flux conservation. \n",
    "\n",
    "mProjExec, used below, is a wrapper around the three main reprojection routines that determines whether mProjectPP or mProject should be used for each image, or calls mProjectQL to override them, as is the case here to provide speed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Reproject the original images to the  frame of the \n",
    "# output FITS header we created\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "# Use the Montage utility 'mProjExec' to reproject the images listed in 'rimages.tbl'\n",
    "# to the coordinate system defined in the 'region.hdr' FITS header.\n",
    "rtn = mProjExec(\"raw\", \"rimages.tbl\", \"../region.hdr\", projdir=\"projected\", quickMode=True)\n",
    "\n",
    "# Print the result of the mProjExec operation for verification\n",
    "print(\"mProjExec:           \" + str(rtn), flush=True)\n",
    "\n",
    "# Use the Montage utility 'mImgtbl' to scan the reprojected images in the 'projected' directory\n",
    "# and generate a metadata table ('pimages.tbl') that describes the coverage of each reprojected image.\n",
    "mImgtbl(\"projected\", \"pimages.tbl\")\n",
    "\n",
    "print(\"mImgtbl (projected): \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Coadd the Reprojected Images\n",
    "\n",
    "Now that we have a set of images all reprojected to a common frame, we can coadd them into a mosaic.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Coadd the projected images without backgound correction.\n",
    "# This step is just to illustrate the need for background correction\n",
    "# and can be omitted.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "# Use the Montage utility 'mAdd' to coadd the reprojected images listed in 'pimages.tbl'.\n",
    "# The combined image will be created according to the frame defined in 'region.hdr'\n",
    "# and saved as 'uncorrected.fits'.\n",
    "rtn = mAdd(\"projected\", \"pimages.tbl\", \"../region.hdr\", \"uncorrected.fits\")\n",
    "\n",
    "print(\"mAdd:    \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Visualize the Image\n",
    "\n",
    "Visualizing the image requires us to render it in JPEG or PNG format.  This involves choosing image parameters such as a stretch and a color table.  Montage provides a general visualization tool, mViewer, which processes single or multiple images for presentation.    mViewer supports  a custom stretch algorithm based on gaussian-transformed histogram equalization, optionally with an extra logarithmic transforms for images with very bright point-like sources.  Many astronomical images share the general characteristics of having a lot of pixels with  a gaussian-like distribution at a low flux level, either from background noise of low-level sky structure,  coupled with a long histogram tail of very bright point-like sources.  If we apply our algorithm to an image such as this, stretching from the -2 or -1 \"sigma\" value of the low-level distribution to the image brightness maximum we usually find a good balance of seeing the low-level structure while still seeing the structure and brightness variations of the bright sources.\n",
    "\n",
    "It supports overlays of various sorts, such as coordinate overlays and labels.\n",
    "\n",
    "mViewer specifications can become lengthy so the module provides three entry mechanisms.  The most terse (used here) is a \"parameter string\" based on the command-line arguments of the original stand-alone C program.  For more complicated descriptions the user can define a JSON string or JSON file.  See examples in the <a style=\"text-decoration: none; color: #c00000\" href=\"mViewer.ipynb\"> Sky Visualization </a> notebook example.\n",
    "\n",
    "We use the built-in IPython.display utility to show the resultant image, which shrinks it to fit the screen. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a PNG rendering of the mosaic data and display it.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "# Use the Montage utility 'mViewer' to create a PNG image of the mosaic\n",
    "rtn = mViewer(\"-ct 1 -gray uncorrected.fits -2s max gaussian-log -out uncorrected.png\", \"\", mode=2)\n",
    "\n",
    "print(\"mViewer: \" + str(rtn), flush=True)\n",
    "\n",
    "# Jupyter IPython display gets unhappy with big images, so shrink the PNG\n",
    "\n",
    "# Read the generated PNG image into memory using OpenCV\n",
    "img = cv.imread(\"uncorrected.png\")\n",
    "\n",
    "# Resize the image to the 'shape' dimensions calculated earlier\n",
    "img = cv.resize(img, shape)\n",
    "\n",
    "# Save the resized image as 'uncorrected_small.png'\n",
    "cv.imwrite('uncorrected_small.png', img)\n",
    "\n",
    "# Display the smaller image in the Jupyter notebook\n",
    "Image(filename='uncorrected_small.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Rectify the Sky Background in the 2MASS Images to a Common Level\n",
    "\n",
    "In the above image you can see vertical stripes.  Even though the images were accurately flux-calibrated, the background levels in the individual image varied due to changes in the sky background brightness of the sky.  This is a common problem in processing ground-based images; differential photometry is easier than absolute. AS there is no physical model for calculating the sky background,  Montage computes the set of minimum adjustments we can make to the individual image backgrounds to bring them all in-line with each other.  Often, this is just a constant offset, or it may include a slight slope.  Anything more and we are starting to fit the sky structure rather than the background differences.\n",
    "\n",
    "The first step is determining the corrections is to analyze the overlap areas between adjacent images.  We determine from the sky coverage image metadata where there are overlaps between images.  Then for each pairwise overlap, we compute the image difference.  There is an explicit assumption here that a such a pair the sources and other real-sky structure match, including flux scales, so the difference should have nothing in it but background differences.  We then fit each difference with a plane (ignoring large excursions just to be safe).\n",
    "\n",
    "Finally, given this set of difference fits, we determine iteratively a global mimimum difference which results in a set of corrections to each image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine the overlaps between images (for background modeling).\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "# Use the Montage utility 'mOverlaps' to determine the overlapping regions between the images\n",
    "# listed in 'pimages.tbl'.\n",
    "rtn = mOverlaps(\"pimages.tbl\", \"diffs.tbl\")\n",
    "\n",
    "print(\"mOverlaps:    \" + str(rtn), flush=True)\n",
    "\n",
    "\n",
    "# Generate difference images and fit them.\n",
    "# This step creates difference images for each pair of overlapping images,\n",
    "# which are then used to model and correct for background variations.\n",
    "rtn = mDiffFitExec(path=\"projected\", tblfile=\"diffs.tbl\", template=\"../region.hdr\", \n",
    "                   diffdir=\"diffs\", fitfile=\"fits.tbl\", archive=\"\", keepAll=0, \n",
    "                   levelOnly=0, noAreas=0, debug=0)\n",
    "\n",
    "print(\"mDiffFitExec: \" + str(rtn), flush=True)\n",
    "\n",
    "\n",
    "# Model the background corrections.\n",
    "# This step uses the fit parameters obtained from the difference images \n",
    "# to generate a model that corrects the background variations across the mosaic.\n",
    "rtn = mBgModel(\"pimages.tbl\", \"fits.tbl\", \"corrections.tbl\")\n",
    "\n",
    "print(\"mBgModel:     \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Background Correcting and Re-Mosaicking\n",
    "\n",
    "Now  we  apply the background corrections to the individual images and regenerate the mosaic.  While we don't attempt to maintain the global total flux (this would be meaningless in any case given the source of the offsets), in our final mosaic is usually lcose to this level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Background correct the projected images.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "# Use the Montage utility 'mBgExec' to apply the background corrections to the reprojected images.\n",
    "rtn = mBgExec(\"projected\", \"pimages.tbl\", \"corrections.tbl\", \"corrected\")\n",
    "\n",
    "print(\"mBgExec:             \" + str(rtn), flush=True)\n",
    "\n",
    "# Create a metadata table for the background-corrected images.\n",
    "rtn = mImgtbl(\"corrected\", \"cimages.tbl\")\n",
    "\n",
    "print(\"mImgtbl (corrected): \" + str(rtn), flush=True)\n",
    "\n",
    "\n",
    "# Coadd the background-corrected, projected images.\n",
    "\n",
    "# Use the Montage utility 'mAdd' to coadd the corrected images listed in 'cimages.tbl'.\n",
    "rtn = mAdd(\"corrected\", \"cimages.tbl\", \"../region.hdr\", \"2massj.fits\")\n",
    "           \n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "print(\"mAdd:                \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Final Image\n",
    "\n",
    "Now when we regenerate and display a PNG for the mosaic. The artifacts have been removed and the low-level structure is preserved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a PNG rendering of the data and display it.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSJ\")\n",
    "\n",
    "# Use the Montage utility 'mViewer' to create a PNG image of the final mosaic\n",
    "rtn = mViewer(\"-ct 1 -gray 2massj.fits -2s max gaussian-log -out 2massj.png\", \"\", mode=2)\n",
    "\n",
    "# Read the generated PNG image into memory using OpenCV\n",
    "img = cv.imread(\"2massj.png\")\n",
    "# Resize the image to the 'shape' dimensions calculated earlier to make it smaller\n",
    "img = cv.resize(img, shape)\n",
    "# Save the resized image as '2massj_small.png'\n",
    "cv.imwrite('2massj_small.png', img)\n",
    "\n",
    "# Display the smaller image in the Jupyter notebook\n",
    "Image(filename='2massj_small.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"#880000\">Full Color and Overlay With MeerKAT Image\n",
    "\n",
    "The above can be packaged up in a Python script with whatever minimum input and defaults you desired.  Repeat the processing for three different wavelengths and you can combine them (and optionally overlays of various sorts) into a full-color image. See our <a style=\"text-decoration: none;\" href=\"mViewer.ipynb\"> mViewer</a> sky visualization notebook for details.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## <font color=\"#880000\">H-Band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean out any old copy of the work tree for 2MASS H, then remake it\n",
    "# and the set of the subdirectories we will need.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "workdir = \"2MASSH\"\n",
    "\n",
    "os.chdir(home)\n",
    "\n",
    "try:\n",
    "    shutil.rmtree(workdir)\n",
    "except:\n",
    "    print(\"Can't delete work tree; probably doesn't exist yet.\\n\", flush=True)\n",
    "\n",
    "os.makedirs(workdir)  \n",
    "\n",
    "os.chdir(workdir)\n",
    "\n",
    "os.makedirs(\"raw\")\n",
    "os.makedirs(\"projected\")\n",
    "os.makedirs(\"diffs\")\n",
    "os.makedirs(\"corrected\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = \"2MASS H\"\n",
    "\n",
    "os.chdir(home + \"/2MASSH\")\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "# Use the Montage utility 'mArchiveDownload' to download the images from the specified dataset\n",
    "# covering the region defined by 'location' and 'size'. The images will be saved in the 'raw' directory.\n",
    "mArchiveDownload(dataset, location, size, \"raw\")\n",
    "rtn = mImgtbl(\"raw\", \"rimages.tbl\")\n",
    "\n",
    "print(\"Raw images:          \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSH\")\n",
    "\n",
    "# Reproject the raw images to match the frame defined by the output FITS header 'region.hdr'\n",
    "mProjExec(\"raw\", \"rimages.tbl\", \"../region.hdr\", projdir=\"projected\", quickMode=True)\n",
    "\n",
    "# Create a metadata table for the reprojected images, storing the result in 'pimages.tbl'\n",
    "mImgtbl(\"projected\", \"pimages.tbl\")\n",
    "\n",
    "# Determine the overlaps between the reprojected images, which is essential for background correction\n",
    "# The overlaps are stored in 'diffs.tbl'\n",
    "mOverlaps(\"pimages.tbl\", \"diffs.tbl\")\n",
    "\n",
    "# Generate difference images for overlapping regions and fit the differences to model the background\n",
    "# The difference images are saved in the 'diffs' directory, and the fit parameters are saved in 'fits.tbl'\n",
    "mDiffFitExec(\"projected\", \"diffs.tbl\", \"../region.hdr\", \"diffs\", \"fits.tbl\")\n",
    "\n",
    "# Model the background corrections based on the fits and create a table of corrections ('corrections.tbl')\n",
    "mBgModel(\"pimages.tbl\", \"fits.tbl\", \"corrections.tbl\", \"\")\n",
    "\n",
    "# Apply the background corrections to the reprojected images, saving the corrected images in 'corrected'\n",
    "mBgExec(\"projected\", \"pimages.tbl\", \"corrections.tbl\", \"corrected\")\n",
    "\n",
    "# Create a metadata table for the corrected images, storing the result in 'cimages.tbl'\n",
    "mImgtbl(\"corrected\", \"cimages.tbl\")\n",
    "\n",
    "#Coadd the background-corrected images into a final mosaic\n",
    "rtn = mAdd(\"corrected\", \"cimages.tbl\", \"../region.hdr\", \"2massh.fits\")\n",
    "\n",
    "# Render the final FITS mosaic as a PNG image using 'mViewer'\n",
    "# The image is rendered in grayscale ('-gray') with a logarithmic stretch ('gaussian-log') and saved as '2massh.png'\n",
    "mViewer(\"-ct 1 -gray 2massh.fits -2s max gaussian-log -out 2massh.png\", \"\", mode=2)\n",
    "\n",
    "img = cv.imread(\"2massh.png\")\n",
    "img = cv.resize(img, shape)\n",
    "cv.imwrite('2massh_small.png', img)\n",
    "\n",
    "# Display the smaller image in the Jupyter notebook\n",
    "Image(filename='2massh_small.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## <font color=\"#880000\">K-Band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean out any old copy of the work tree for 2MASS H, then remake it\n",
    "# and the set of the subdirectories we will need.\n",
    "\n",
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "workdir = \"2MASSK\"\n",
    "\n",
    "os.chdir(home)\n",
    "\n",
    "try:\n",
    "    shutil.rmtree(workdir)\n",
    "except:\n",
    "    print(\"Can't delete work tree; probably doesn't exist yet.\\n\", flush=True)\n",
    "\n",
    "os.makedirs(workdir)  \n",
    "\n",
    "os.chdir(workdir)\n",
    "\n",
    "os.makedirs(\"raw\")\n",
    "os.makedirs(\"projected\")\n",
    "os.makedirs(\"diffs\")\n",
    "os.makedirs(\"corrected\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "dataset = \"2MASS K\"\n",
    "\n",
    "# Use the Montage utility 'mArchiveDownload' to download the images from the specified dataset\n",
    "# covering the region defined by 'location' and 'size'. The images will be saved in the 'raw' directory.\n",
    "mArchiveDownload(dataset, location, size, \"raw\")\n",
    "\n",
    "# Use the Montage utility 'mImgtbl' to scan the downloaded images in the 'raw' directory\n",
    "# and generate a metadata table ('rimages.tbl') that describes the coverage and properties of each image.\n",
    "rtn = mImgtbl(\"raw\", \"rimages.tbl\")\n",
    "\n",
    "print(\"Raw images:          \" + str(rtn), flush=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home + \"/2MASSK\")\n",
    "\n",
    "# Reproject the raw images to match the frame defined by the output FITS header 'region.hdr'.\n",
    "# The reprojected images will be saved in the 'projected' directory.\n",
    "mProjExec(\"raw\", \"rimages.tbl\", \"../region.hdr\", projdir=\"projected\", quickMode=True)\n",
    "\n",
    "# Create a metadata table for the reprojected images, storing the result in 'pimages.tbl'\n",
    "mImgtbl(\"projected\", \"pimages.tbl\")\n",
    "\n",
    "# Determine the overlaps between the reprojected images, which is essential for background correction.\n",
    "# The overlaps are stored in 'diffs.tbl'\n",
    "mOverlaps(\"pimages.tbl\", \"diffs.tbl\")\n",
    "\n",
    "# Generate difference images for overlapping regions and fit the differences to model the background.\n",
    "# The difference images are saved in the 'diffs' directory, and the fit parameters are saved in 'fits.tbl'\n",
    "mDiffFitExec(\"projected\", \"diffs.tbl\", \"../region.hdr\", \"diffs\", \"fits.tbl\")\n",
    "\n",
    "# Model the background corrections based on the fits and create a table of corrections ('corrections.tbl')\n",
    "mBgModel(\"pimages.tbl\", \"fits.tbl\", \"corrections.tbl\", \"\")\n",
    "\n",
    "# Apply the background corrections to the reprojected images, saving the corrected images in 'corrected'\n",
    "mBgExec(\"projected\", \"pimages.tbl\", \"corrections.tbl\", \"corrected\")\n",
    "\n",
    "# Create a metadata table for the corrected images, storing the result in 'cimages.tbl'\n",
    "mImgtbl(\"corrected\", \"cimages.tbl\")\n",
    "\n",
    "rtn = mAdd(\"corrected\", \"cimages.tbl\", \"../region.hdr\", \"2massk.fits\")\n",
    "\n",
    "\n",
    "# Render the final FITS mosaic as a PNG image using 'mViewer'.\n",
    "# The image is rendered in grayscale ('-gray') with a logarithmic stretch ('gaussian-log') and saved as '2massk.png'.\n",
    "mViewer(\"-ct 1 -gray 2massk.fits -2s max gaussian-log -out 2massk.png\", \"\", mode=2)\n",
    "\n",
    "img = cv.imread(\"2massk.png\")\n",
    "img = cv.resize(img, shape)\n",
    "cv.imwrite('2massk_small.png', img)\n",
    "\n",
    "# Display the smaller image in the Jupyter notebook\n",
    "Image(filename='2massk_small.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <font color=\"#880000\">Create 2MASS Color Mosaic\n",
    "\n",
    "Now that we have FITS images for all three (J, H and K) 2MASS wavelengths, we can create a color-composite 2MASS PNG for the region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home)\n",
    "\n",
    "# Use the Montage utility 'mViewer' to create a color composite image from the three 2MASS FITS files.\n",
    "rtn = mViewer(\"-blue  2MASSJ/2massj.fits -1s max gaussian-log \" \n",
    "              \"-green 2MASSH/2massh.fits -1s max gaussian-log \" \n",
    "              \"-red   2MASSK/2massk.fits -1s max gaussian-log -out 2MASS.png\", \n",
    "              \"\", mode=2)\n",
    "\n",
    "# Set contrast and brightness values for post-processing the image\n",
    "contrast   = 2.0\n",
    "brightness = -1.0\n",
    "\n",
    "os.chdir(home)\n",
    "\n",
    "# Read the generated color composite PNG image into memory using OpenCV\n",
    "tmass = cv.imread(\"2MASS.png\")\n",
    "\n",
    "# Adjust the contrast and brightness of the image using OpenCV's addWeighted function\n",
    "# 'contrast' and 'brightness' are applied to enhance the visibility of the image details\n",
    "tmassStretch = cv.addWeighted(tmass, contrast, tmass, 0, brightness)\n",
    "\n",
    "cv.imwrite('2massStretch.png', tmassStretch)\n",
    "\n",
    "img = cv.resize(tmassStretch, shape)\n",
    "cv.imwrite('2MASS_small.png', img)\n",
    "\n",
    "# Display the smaller, processed image in the Jupyter notebook\n",
    "Image('2MASS_small.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and a grayscale image with exactly the same coverage from MeerKAT:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "os.chdir(home)\n",
    "\n",
    "# Use the Montage utility 'mViewer' to create a grayscale PNG rendering of the MeerKAT FITS file.\n",
    "rtn = mViewer(\"-gray MeerKAT_Galactic_Centre_1284MHz-StokesI.fits \"\n",
    "              \"-1s max gaussian-log -out MeerKAT.png\",\n",
    "              \"\", mode=2)\n",
    "\n",
    "# Set contrast and brightness values for post-processing the image.\n",
    "contrast   = 3.0\n",
    "brightness = -128.0\n",
    "\n",
    "meerkat = cv.imread('MeerKAT.png')\n",
    "\n",
    "# Adjust the contrast and brightness of the image using OpenCV's addWeighted function.\n",
    "meerkatStretch = cv.addWeighted(meerkat, contrast, meerkat, 0, brightness)\n",
    "cv.imwrite('meerkatStretch.png', meerkatStretch)\n",
    "\n",
    "\n",
    "# Resize the adjusted image to the dimensions specified in 'shape' for easier display.\n",
    "img = cv.resize(meerkatStretch, shape)\n",
    "cv.imwrite('MeerKAT_small.png', img)\n",
    "\n",
    "# Display the smaller, processed image in the Jupyter notebook.\n",
    "Image('MeerKAT_small.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## <font color=\"#880000\"> Overlay 2MASS and MeerKAT Images and Blend Together\n",
    "\n",
    "Finally, we use OpenCV to blend the images.  We are repeating some of the steps in the last two cells so that this cell can stand alone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "contrast   = 2.0\n",
    "brightness = -1.0\n",
    "\n",
    "os.chdir(home)\n",
    "\n",
    "tmassPath = os.path.join(home, '2MASS.png')\n",
    "tmass = cv.imread(tmassPath)\n",
    "\n",
    "tmassStretch = cv.addWeighted(tmass, contrast, tmass, 0, brightness)\n",
    "\n",
    "cv.imwrite('2massStretch.png', tmassStretch)\n",
    "\n",
    "# -------------------\n",
    "   \n",
    "contrast   = 3.0\n",
    "brightness = -128.0\n",
    "\n",
    "meerkatPath = os.path.join(home, 'MeerKAT.png')\n",
    "meerkat = cv.imread(meerkatPath)\n",
    "\n",
    "meerkatStretch = cv.addWeighted(meerkat, contrast, meerkat, 0, brightness)\n",
    "\n",
    "cv.imwrite('meerkatStretch.png', meerkatStretch)\n",
    "\n",
    "# -------------------\n",
    "# Set the blending parameters for combining the 2MASS and MeerKAT images\n",
    "alpha = 0.6\n",
    "beta  = 1.0 - alpha\n",
    "gamma = 48.0\n",
    "\n",
    "imgBlend = cv.addWeighted(tmassStretch, alpha, meerkatStretch, beta, gamma)\n",
    "\n",
    "cv.imwrite('MeerKAT_2MASS.png', imgBlend)\n",
    "\n",
    "img = cv.imread(\"MeerKAT_2MASS.png\")\n",
    "img = cv.resize(img, shape)\n",
    "cv.imwrite('MeerKAT_2MASS_small.png', img)\n",
    "\n",
    "Image(filename=\"MeerKAT_2MASS_small.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PNGs that we display here are essentially thumbnails; they have been shrunken by a factor of 10 in both dimensions.  See the PNG files on disk for full resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.datetime.now()\n",
    "print(now.strftime(\"%Y-%m-%d %H:%M:%S\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"#880000\">Acknowledgements\n",
    "\n",
    "Montage is funded by the National Science Foundation under Grant Numbers ACI-1440620,1642453 and 1835379, and was previously funded by the National Aeronautics and Space Administration's Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"#880000\">Need help?\n",
    "\n",
    "Contact the Montage help desk at http://vaoweb3.ipac.caltech.edu/cgi-bin/Helpdesk/nph-genTicketForm?projname=Montage&projmail=montage@ipac.caltech.edu"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}