diff --git a/requirements.txt b/requirements.txt
index 9da18a57..39681e9b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,5 +12,8 @@ reproject
 dust_extinction
 gala
 synphot
-pyvo==1.2.0
+pyvo
 pvextractor
+photutils
+requests
+ccdproc
diff --git a/tutorials/image-registration/image-registration.ipynb b/tutorials/image-registration/image-registration.ipynb
new file mode 100644
index 00000000..f21de6c2
--- /dev/null
+++ b/tutorials/image-registration/image-registration.ipynb
@@ -0,0 +1,763 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "5524f2e1",
+   "metadata": {},
+   "source": [
+    "# Image Registration and Combination\n",
+    "\n",
+    "## Author\n",
+    "[C. E. Brasseur](https://ceb8.github.io/)\n",
+    "\n",
+    "## Learning Goals\n",
+    "* Correcting image World Coordinate System (WCS) mapping\n",
+    "* Aligning one image to another\n",
+    "* Co-adding images\n",
+    "\n",
+    "\n",
+    "## Keywords\n",
+    "FITS, WCS, DaoFind, ccdproc, source extraction, wcs fitting, reproject, image alignment, coordinate crossmatch\n",
+    "\n",
+    "\n",
+    "## Summary\n",
+    "Often we want to combine images, either by co-adding to improve signal-to-noise, or subtracting to bring out variable stars, or combining different wavebands to colorize an image.  In this tutorial I will lead you through this entire process, starting with two images from the Las Cumbres Observatory.\n",
+    "\n",
+    "First we must correct the pixel-to-world coordinate mapping, for which we will crossmatch sources in the image against the Gaia catalog. These the images then must be projected onto the same alignement, which we do using `reproject`. Finally we use `ccdproc` to combine the images.\n",
+    "\n",
+    "I assume the learner is comfortable using FITS files and WCS objects (see the [FITS-images](https://learn.astropy.org/tutorials/FITS-images.html) and [FITS-cubes](https://learn.astropy.org/tutorials/FITS-cubes.html) tutorials).\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "af185b76",
+   "metadata": {},
+   "source": [
+    "## Table of Contents\n",
+    "\n",
+    "- [1. Imports](#1.-Imports)\n",
+    "- [2. Download the data](#2.-Download-the-data)\n",
+    "- [3. Correct the WCS info](#3.-Correct-the-WCS-info)\n",
+    "- [4. Image Registration](#4.-Image-Registration)\n",
+    "- [5. Coadd the images](#5.-Coadd-the-images)\n",
+    "- [6. Save stacked image](#6.-Save-stacked-image)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c5d78420",
+   "metadata": {},
+   "source": [
+    "## 1. Imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1271c318",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# For downloading observations\n",
+    "import requests\n",
+    "import json\n",
+    "\n",
+    "# For Querying Gaia\n",
+    "from astroquery.mast import Catalogs\n",
+    "\n",
+    "# For source finding\n",
+    "from photutils.detection import DAOStarFinder\n",
+    "\n",
+    "# For WCS fitting\n",
+    "from astropy.wcs.utils import fit_wcs_from_points\n",
+    "\n",
+    "# For image alignment\n",
+    "from reproject import reproject_interp\n",
+    "\n",
+    "# For image co-adding\n",
+    "from ccdproc import CCDData, Combiner\n",
+    "\n",
+    "# For interaction with FITS image files\n",
+    "from astropy.io import fits\n",
+    "from astropy.wcs import WCS\n",
+    "\n",
+    "# Other useful astropy imports\n",
+    "from astropy.stats import sigma_clipped_stats\n",
+    "from astropy.table import Table\n",
+    "from astropy.coordinates import SkyCoord\n",
+    "import astropy.units as u\n",
+    "\n",
+    "# For plotting\n",
+    "import matplotlib\n",
+    "%matplotlib inline\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib.colors import LogNorm"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4069768e",
+   "metadata": {},
+   "source": [
+    "## 2. Download the data\n",
+    "\n",
+    "For this tutorial we will use two images taken by the Las Cumbres Observatory in 2015. We will download two images of the globular cluster NGC 1866 in the V-band. These images can be searched for at the [LCO Archive](https://archive.lco.global/), however for this tutorial we will assume we already know the files we want."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2c2092d7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_lco_observation(frame):\n",
+    "    \"\"\"\n",
+    "    Given a frame value download the associated LCO obvservation fits file.\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    frame_url = f\"https://archive-api.lco.global/frames/{frame}\"\n",
+    "    result = requests.get(frame_url)\n",
+    "    \n",
+    "    file_info = json.loads(result.content)\n",
+    "    data = requests.get(file_info[\"url\"])\n",
+    "    \n",
+    "    with open(file_info[\"filename\"], 'wb') as FLE:\n",
+    "        FLE.write(data.content)\n",
+    "        \n",
+    "    return file_info[\"filename\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10ff9cd3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "file_1 = get_lco_observation(915222)\n",
+    "file_2 = get_lco_observation(934536)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "364ae1ee",
+   "metadata": {},
+   "source": [
+    "### Plotting the data\n",
+    "\n",
+    "Here we open both files and plot them individually and together. We can see from the composite image that the two images are offset from each other, so we will need to align them befor we can co-add them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f72ad44d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hdu_1 = fits.open(file_1)\n",
+    "img_1 = hdu_1[1].data\n",
+    "hdr_1 = hdu_1[1].header\n",
+    "hdu_1.close()\n",
+    "\n",
+    "hdu_2 = fits.open(file_2)\n",
+    "img_2 = hdu_2[1].data\n",
+    "hdr_2 = hdu_2[1].header\n",
+    "hdu_2.close()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bdd1c518",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axs = plt.subplots(1, 3, figsize=(12,4))\n",
+    "\n",
+    "for ax in axs:\n",
+    "    ax.set_axis_off()\n",
+    "\n",
+    "axs[0].imshow(img_1, cmap='Oranges_r', norm=LogNorm(vmin=10,vmax=60))\n",
+    "axs[2].imshow(img_1, cmap='Oranges_r', norm=LogNorm(vmin=10,vmax=60), alpha=0.5)\n",
+    "\n",
+    "axs[1].imshow(img_2, cmap='Blues_r', norm=LogNorm(vmin=10,vmax=60))\n",
+    "axs[2].imshow(img_2, cmap='Blues_r', norm=LogNorm(vmin=10,vmax=60), alpha=0.5)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0f977c42",
+   "metadata": {},
+   "source": [
+    "## 3. Correct the WCS info\n",
+    "\n",
+    "Before we can align the images we need to make sure they both have correct WCS information.The observation files include WCS information, however it is not accurate enough for our purposes.\n",
+    "\n",
+    "We will first focus on the process for a single file (`hdr_1` and `img_1`).\n",
+    "\n",
+    "To do this we need to identify a handful of bright sources in our image (x,y pixel coordinates) and match them to canonical world coordinates (ra,dec) for those sources. We can then use the list of x,y and canonical ra,dec pairs to build the corrected WCS object.\n",
+    "\n",
+    "**Note:** Loading the WCS information from the observation header gives a warning. This is not a problem, it is astropy telling us it has \"fixed\" the header WCS keywords to make them conform to the official standards."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1f0375b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init_wcs = WCS(hdr_1)\n",
+    "print(init_wcs)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c1074d6f",
+   "metadata": {},
+   "source": [
+    "### a. Finding source locations in the image.\n",
+    "\n",
+    "We will use `DAOStarFinder` (from the [detection module](https://photutils.readthedocs.io/en/stable/detection.html) in photutils) to identify sources within our image. We estimate the background level (median) and background noise (std)using [sigma-clipped statistics](https://docs.astropy.org/en/stable/api/astropy.stats.sigma_clipped_stats.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "68bcd96d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mean, median, std = sigma_clipped_stats(img_1, sigma=3.0)    \n",
+    "print(f\"Image mean: {mean:.1f}\\nImage median: {median:.1f}\\nImage std: {std:.1f}\") "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "71f47efc",
+   "metadata": {},
+   "source": [
+    "We now create a `DAOStarFinder` looking for stars with a FWHM of ~6 px and peaks at least 120 standard deviations above the background. We set the threshold so high because we want to use only a handful of the brightest stars for our reference objects."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a58b89b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "daofind = DAOStarFinder(fwhm=6, threshold=120*std) # We only want the brightest sources\n",
+    "sources = daofind(img_1 - median)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "25763a68",
+   "metadata": {},
+   "source": [
+    "### b. Gaia sources\n",
+    "\n",
+    "We will use the [Gaia Catalog](https://gea.esac.esa.int/archive/) through [astroquery](https://astroquery.readthedocs.io/en/latest/mast/mast.html#catalog-queries) to get the canonical sky coordinates for our reference stars."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "88d96bfc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Pull out the center coordinate of the image\n",
+    "center_ra,center_dec = init_wcs.wcs.crval\n",
+    "coord = SkyCoord(ra=center_ra, dec=center_dec, unit=(u.degree, u.degree), frame='icrs')\n",
+    "\n",
+    "# Performing the query\n",
+    "gaia_catalog = Catalogs.query_region(coordinates=coord, radius=0.095*u.deg, catalog=\"Gaia\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e32268ac",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Filtering down to brightest sources not in the dense center\n",
+    "\n",
+    "gaia_bright = gaia_catalog[gaia_catalog[\"phot_g_mean_mag\"] < 15]\n",
+    "gaia_bright = gaia_bright[gaia_bright['distance'] > 3]\n",
+    "\n",
+    "gaia_bright['coord'] = SkyCoord(gaia_bright['ra'], gaia_bright['dec'], unit=u.deg)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "28a1db92",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Using the observation WCS to get the image coordinates for the Gaia sources\n",
+    "\n",
+    "pxs = init_wcs.world_to_pixel(gaia_bright['coord'])\n",
+    "gaia_bright[\"x\"] = pxs[0]\n",
+    "gaia_bright[\"y\"] = pxs[1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9105c5dd",
+   "metadata": {},
+   "source": [
+    "### c. Plotting image and Gaia sources together\n",
+    "\n",
+    "Before we embark on correcting the image WCS object we plot both the image sources from `DAOStarFinder`, and the canonical sources from Gaia. We can see that there are a few image sources (blue) without Gaia counterparts (red), and that the Gaia sources are offset from the image sources. The offset is due to errors in the image WCS, which is what we are aiming to correct."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "75270527",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(7,7))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(img_1, cmap='gray', norm=LogNorm(vmin=10,vmax=60))\n",
+    "\n",
+    "ax.scatter(sources['xcentroid'], sources['ycentroid'], ec='#44aeed', fc=\"none\", s=150, lw=2, label=\"Image source\")\n",
+    "ax.scatter(gaia_bright[\"x\"], gaia_bright[\"y\"], ec='#e71f71', fc=\"none\", s=50, lw=2, label=\"Gaia source\")\n",
+    "\n",
+    "ax.legend(fontsize=13)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "be100d36",
+   "metadata": {},
+   "source": [
+    "### d. Catalog crossmatch\n",
+    "\n",
+    "Now that we have our reference sources from the image and Gaia, we need to pair the two lists properly. To do this we will crossmatch the list of bright Gaia sources with the list of image sources using [`match_to_catalog_sky`](https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord.match_to_catalog_sky). \n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d82b1607",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Getting sky coordinates for image sources\n",
+    "sources[\"coord\"] = init_wcs.pixel_to_world(sources['xcentroid'], sources['ycentroid'])\n",
+    "\n",
+    "# Performing the cross match\n",
+    "idx, d2d, d3d = gaia_bright['coord'].match_to_catalog_sky(sources[\"coord\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9ac1f459",
+   "metadata": {},
+   "source": [
+    "`match_to_catalog_sky` returns the indexes of the input catalog (`sources`) that correspondes to matched source in base catalog (`gaia_bright`), and the 2- and 3-D distances between the sources. We use that indesing to update the pixel positions of the `gaia_bright` sources, and add a column with the 2D distances, which are the calculated offsets between the gaia source pixel positions and the true image pixel positions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "63552490",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Update the x/y columns in the gaia catalog\n",
+    "gaia_bright[\"x\"] = sources['xcentroid'][idx]\n",
+    "gaia_bright[\"y\"] = sources['ycentroid'][idx]\n",
+    "\n",
+    "# Add the distance column (we have no distances so no need for the 3d distances)\n",
+    "gaia_bright[\"d2d\"] = d2d.deg\n",
+    "\n",
+    "print(\"Nearest Source Distances:\")\n",
+    "for row in gaia_bright:\n",
+    "    print(f\"{row['d2d']:.3f} deg\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eb04777c",
+   "metadata": {},
+   "source": [
+    "We print out the distances to make sure we're indeed connecting the correct catalog entries. All of them look good (i.e. small and similar to each other), so we can move on to use [`fit_wcs_from_points`](https://docs.astropy.org/en/stable/api/astropy.wcs.utils.fit_wcs_from_points.html) to build the corrected WCS."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e12414b4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "corrected_wcs_1 = fit_wcs_from_points([gaia_bright[\"x\"],gaia_bright[\"y\"]], gaia_bright[\"coord\"])\n",
+    "print(corrected_wcs_1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1eb477b7",
+   "metadata": {},
+   "source": [
+    "Now we can update the Gaia x/y coordinates, using the corrected WCS, and plot the sources again. This time they line up correctly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "daa56e21",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pxs = corrected_wcs_1.world_to_pixel(gaia_bright['coord'])\n",
+    "gaia_bright[\"x\"] = pxs[0]\n",
+    "gaia_bright[\"y\"] = pxs[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4317cd77",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(7,7))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(img_1, cmap='gray', norm=LogNorm(vmin=10,vmax=60))\n",
+    "\n",
+    "ax.scatter(sources['xcentroid'], sources['ycentroid'], ec='#44aeed', fc=\"none\", s=150, lw=2, label=\"Image source\")\n",
+    "ax.scatter(gaia_bright[\"x\"], gaia_bright[\"y\"], ec='#e71f71', fc=\"none\", s=50, lw=2, label=\"Gaia source\")\n",
+    "\n",
+    "ax.legend(fontsize=13)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fe42bccc",
+   "metadata": {},
+   "source": [
+    "### e. Doing the same for the other image\n",
+    "\n",
+    "Next we need to go through all the same steps for the second image."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "34eafced",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Finding source locations on the image\n",
+    "\n",
+    "init_wcs = WCS(hdr_2)\n",
+    "\n",
+    "mean, median, std = sigma_clipped_stats(img_2, sigma=3.0)    \n",
+    "\n",
+    "daofind = DAOStarFinder(fwhm=6.2, threshold=120*std) # We only want the brightest sources\n",
+    "sources = daofind(img_2 - median)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2f88c903",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Using the observation WCS to get the image coordinates for the Gaia sources\n",
+    "\n",
+    "pxs = init_wcs.world_to_pixel(gaia_bright['coord'])\n",
+    "gaia_bright[\"x\"] = pxs[0]\n",
+    "gaia_bright[\"y\"] = pxs[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fe7e0d56",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(5,5))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(img_2, cmap='gray', norm=LogNorm(vmin=10,vmax=60))\n",
+    "\n",
+    "ax.scatter(sources['xcentroid'], sources['ycentroid'], ec='#44aeed', fc=\"none\", s=150, lw=2, label=\"Image source\")\n",
+    "ax.scatter(gaia_bright[\"x\"], gaia_bright[\"y\"], ec='#e71f71', fc=\"none\", s=50, lw=2, label=\"Gaia source\")\n",
+    "\n",
+    "ax.legend(fontsize=13)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "60832d29",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Catalog crossmatch\n",
+    "sources[\"coord\"] = init_wcs.pixel_to_world(sources['xcentroid'], sources['ycentroid'])\n",
+    "\n",
+    "idx, d2d, d3d = gaia_bright['coord'].match_to_catalog_sky(sources[\"coord\"])\n",
+    "\n",
+    "gaia_bright[\"x\"] = sources['xcentroid'][idx]\n",
+    "gaia_bright[\"y\"] = sources['ycentroid'][idx]\n",
+    "gaia_bright[\"d2d\"] = d2d.deg"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e0981f5f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Correcting the WCS\n",
+    "corrected_wcs_2 = fit_wcs_from_points([gaia_bright[\"x\"],gaia_bright[\"y\"]], gaia_bright[\"coord\"])\n",
+    "\n",
+    "pxs = corrected_wcs_2.world_to_pixel(gaia_bright['coord'])\n",
+    "gaia_bright[\"x\"] = pxs[0]\n",
+    "gaia_bright[\"y\"] = pxs[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7ee2e3ff",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(5,5))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(img_2, cmap='gray', norm=LogNorm(vmin=10,vmax=60))\n",
+    "\n",
+    "ax.scatter(sources['xcentroid'], sources['ycentroid'], ec='#44aeed', fc=\"none\", s=150, lw=2, label=\"Image source\")\n",
+    "ax.scatter(gaia_bright[\"x\"], gaia_bright[\"y\"], ec='#e71f71', fc=\"none\", s=50, lw=2, label=\"Gaia source\")\n",
+    "\n",
+    "ax.legend(fontsize=13)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3c86735e",
+   "metadata": {},
+   "source": [
+    "## 4. Image Registration\n",
+    "\n",
+    "At this point we have two images, each with a WCS object that we know is correct and consistent. The next step is to use the [`reproject`](https://reproject.readthedocs.io/en/stable/) package to register the images together. \n",
+    "\n",
+    "\n",
+    "We will align `img_2` onto `img_1`. This means that we will interpolate `img_2` onto the pixel grid of `img_1`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c6f234c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "img_2_aligned, footprint = reproject_interp((img_2, corrected_wcs_2), corrected_wcs_1, shape_out=img_2.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9a54608c",
+   "metadata": {},
+   "source": [
+    "Plotting the images together, we can see they are no longer offset. \n",
+    "\n",
+    "Note also the orange border on the left and top. This shows us how `img_2` had to be shifted to align with `img_1`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dcc2954e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(7,7))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(img_1, cmap='Oranges_r', norm=LogNorm(vmin=10,vmax=60), alpha=0.5)\n",
+    "ax.imshow(img_2_aligned, cmap='Blues_r', norm=LogNorm(vmin=10,vmax=60), alpha=0.5)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5da6207a",
+   "metadata": {},
+   "source": [
+    "## 5. Coadd the images\n",
+    "\n",
+    "A common reason to align images is to allow the user to coadd the images to increase signal-to-noise and decrease the background noise. Here we will use the [image combination](https://ccdproc.readthedocs.io/en/latest/image_combination.html) functionality from [ccdproc](https://ccdproc.readthedocs.io/en/latest/index.html), which offers a variety of useful options. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "af0f6b77",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load each image into a CCDData frame\n",
+    "\n",
+    "ccd_frame_1 = CCDData(img_1, unit=u.dimensionless_unscaled)\n",
+    "ccd_frame_2 = CCDData(img_2_aligned, unit=u.dimensionless_unscaled)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "85d7684b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Build the image combiner\n",
+    "\n",
+    "combiner = Combiner([ccd_frame_1, ccd_frame_2])\n",
+    "combiner.data_arr.mask[np.isnan(combiner.data_arr)]=True"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d93d9346",
+   "metadata": {},
+   "source": [
+    "We produce a combined image by averaging the input images."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7881465b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coadd_img = combiner.average_combine()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "08973bae",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(7,7))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(coadd_img, cmap='gray', norm=LogNorm(vmin=10,vmax=60))\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "62aa8d62",
+   "metadata": {},
+   "source": [
+    "We can take the mediane of the input images instead."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e058c10c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coadd_img = combiner.median_combine()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d33c301d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(7,7))\n",
+    "ax.set_axis_off()\n",
+    "\n",
+    "ax.imshow(coadd_img, cmap='gray', norm=LogNorm(vmin=10,vmax=60))\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eb57b8ac",
+   "metadata": {},
+   "source": [
+    "## 6. Save stacked image\n",
+    "\n",
+    "Now that we've processed our image, we can save the result to a new FITS file. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8788878b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hdu = fits.PrimaryHDU(coadd_img, header=corrected_wcs_1.to_header())\n",
+    "hdu.writeto(\"combined_img.fits\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a065772c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython"
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tutorials/image-registration/requirements.txt b/tutorials/image-registration/requirements.txt
new file mode 100644
index 00000000..a64bad67
--- /dev/null
+++ b/tutorials/image-registration/requirements.txt
@@ -0,0 +1,9 @@
+numpy
+requests
+json
+matplotlib
+astropy
+astroquery
+reproject
+ccdproc
+photutils