notebooks/05_composites.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RGBs and Other Composites\n",
    "\n",
    "Meteorological satellite instrument bands provide a look at the Earth at various wavelengths. Each band represents the electro-magnetic radiation observed at one wavelength or range of wavelengths. While this is very useful, flexible, and requires a very limited amount of processing to visualize, it does require that the viewer know what to expect from each band or wavelength. In addition to viewing individual bands, another common option for viewing satellite data is to combine multiple bands in to one image such as RGB (Red, Green, Blue) images. By assigning certain bands to particular channels in the output RGB image we can provide an easy way to identify features in the data. Instead of having to know what a particular wavelength is sensitive to, we can now focus on what colors in particular RGB recipe mean.\n",
    "\n",
    "Using our ABI CONUS data from before, let's look at some of the built-in RGB composites provided by Satpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ../init_notebook.py\n",
    "from satpy import Scene\n",
    "from glob import glob\n",
    "\n",
    "filenames = glob('../data/abi_l1b/20180511_texas_fire_abi_l1b_conus/*.nc')\n",
    "scn = Scene(reader='abi_l1b', filenames=filenames)\n",
    "scn.available_composite_names()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above list shows what composites Satpy knows how to make and that it has the inputs for. A lot of the recipes in Satpy come from standard recipes from science agencies like EUMETSAT, NOAA, or NASA. Let's first look at the `airmass` RGB. The `airmass` RGB is made of the following bands:\n",
    "\n",
    "```\n",
    "R: C08 - C10\n",
    "G: C12 - C13\n",
    "B: C08\n",
    "```\n",
    "\n",
    "The red channel is the difference between the C08 (6.185µm) and C10 (7.34µm) bands, the green channel is the difference between the C12 (9.61µm) and C13 (10.35µm) bands, and the blue channel is the C08 (6.185µm) band.\n",
    "\n",
    "The EUMETSAT EUMeTrain manual summarizes the use of the `airmass` RGB as:\n",
    "\n",
    "> As both water vapor channels (WV 6.2 µm and WV 7.3 µm) are included in this composite, the main applications are the detection of dynamic processes, such as rapid cyclogenesis, jet streams and PV anomalies.\n",
    "\n",
    "Read more on the [EUMeTrain website](http://www.eumetrain.org/data/3/306/navmenu.php?tab=9&page=4.0.0).\n",
    "\n",
    "In Satpy, we can load preconfigured composites like any normal file-provided product. Satpy will look at the recipe, determine what bands need to loaded, what modifications/corrections need to be applied to satisfy the recipe, and then give us the `DataArray` representing that RGB back."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scn.load(['airmass'])\n",
    "scn['airmass']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key part of the above xarray output is that this is the first 3 dimensional DataArray we've worked with in Satpy with the `(bands: 3, y: 1500, x: 2500)` dimensions. The `bands` dimension tells us what each \"band\" represents. In this case Red, Green, and Blue channels of an RGB image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scn['airmass'].coords['bands']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's save this product to a PNG image to get an idea of what an `airmass` RGB looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.diagnostics import ProgressBar\n",
    "with ProgressBar():\n",
    "    scn.save_datasets(writer='simple_image')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see in the image how the various bands combine in different ways and certain colors can represent specific states of the Earth.\n",
    "\n",
    "Earlier in this tutorial we talked about \"enhancing\" data; the process of scaling data from observed scientific data values to image values (ex. 0-255 integers). This process was automatically done for us by the `save_datasets`. If we want to display our data in a matplotlib figure we need to do this manually. Satpy provides a `get_enhanced_image` utility function to simplify this.\n",
    "\n",
    "We can tell xarray what dimension of our `DataArray` represents the \"RGB\" axis. We do this by passing the `rgb` keyword argument with the name of the dimension. In Satpy this dimension should always be called `'bands'`. Without this matplotlib would not be able to plot our RGB properly.\n",
    "\n",
    "Lastly, we'll tell the `imshow` method that we *know* that our RGB will be between 0 and 1 with the `vmin` and `vmax` keyword arguments so there is no need to compute the min/max of the data.\n",
    "\n",
    "Let's plot our `airmass` RGB composite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "from satpy.writers import get_enhanced_image\n",
    "\n",
    "plt.figure()\n",
    "img = get_enhanced_image(scn['airmass'])\n",
    "# get DataArray out of `XRImage` object\n",
    "img_data = img.data\n",
    "img_data.plot.imshow(rgb='bands', vmin=0, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compositing Different Resolutions\n",
    "\n",
    "When generating the `airmass` RGB we were fortunate in that all of the input data was the same resolution. Let's look at the `natural_color` RGB now, where the input resolutions differ. The `natural_color` RGB is useful as described by EUMeTrain for seeing snow and ice:\n",
    "\n",
    "> Snow on the ground as well as ice over mountains, frozen lakes and sea ice appear cyan in the Natural Colour RGB images.\n",
    "\n",
    "See the [EUMeTrain website](http://www.eumetrain.org/resources/MSG_natcolour_cyan_snow.html) for more information.\n",
    "\n",
    "For the ABI instrument this composite is a combination of the following channels:\n",
    "\n",
    "```\n",
    "R: C05 (1000m)\n",
    "G: C03 (1000m)\n",
    "B: C02 (500m)\n",
    "```\n",
    "\n",
    "This means that we are trying to combine arrays of different shapes. The only way we know to make data the same resolution and make it easier to compare or combine, is to resample. Let's load this composite and see how Satpy handles generating a composite like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scn.load(['natural_color'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'natural_color' in scn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scn.missing_datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We didn't do anything wrong, but Satpy doesn't have a `natural_color` composite in the `Scene`. Satpy knows that the reason it wasn't able to generate this composite was because the data was not all the same resolution. The `Scene` knows that we need to resample before we can be able to generate this composite. Let's do that now by using the `native` resampler. We'll also use the `min_area` method to aggregate the data to the lowest resolution in the Scene; 2km in this case (from the `airmass` dependencies):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_scn = scn.resample(scn.min_area(), resampler='native')\n",
    "new_scn['natural_color']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Once resampling had completed, the `Scene` attempted to regenerate any composites that were previous requested. Once successfully generated the `natural_color` RGB shows up like any other product. Let's follow the same procedure we used for plotting the `airmass` RGB, but for the `natural_color` RGB. Note we are using the `new_scn` because the `natural_color` RGB only exists in this resampled `Scene`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with ProgressBar():\n",
    "    plt.figure()\n",
    "    img = get_enhanced_image(new_scn['natural_color'])\n",
    "    # get DataArray out of `XRImage` object\n",
    "    img_data = img.data\n",
    "    img_data.plot.imshow(vmin=0, vmax=1, rgb='bands')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom RGBs\n",
    "\n",
    "Now we've learned how to look at the available builtin composites, load them, and generate them. Let's make our own custom RGB. Let's recreate our Scene from before and load three channels that we'd like to put in an RGB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = glob('../data/abi_l1b/20180511_texas_fire_abi_l1b_conus/*.nc')\n",
    "scn = Scene(reader='abi_l1b', filenames=filenames)\n",
    "scn.available_dataset_names()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_band = 'EDITME'\n",
    "g_band = 'EDITME'\n",
    "b_band = 'EDITME'\n",
    "scn.load([r_band, g_band, b_band])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `Scene` object now has three channels loaded. We'll now use the special `GenericCompositor` from Satpy's builtin composites to join these arrays together for us and make an RGB. Let's see what happens when we attempt this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from satpy.composites import GenericCompositor\n",
    "comp = GenericCompositor('my_rgb')\n",
    "my_rgb = comp((scn[r_band], scn[g_band], scn[b_band]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the channels you choose are different resolutions, you should have gotten an `IncompatibleAreas` exception from the above operation. Just like with `natural_color` before, this is Satpy's compositors telling us that our inputs are at different resolutions and they need to be resampled.\n",
    "\n",
    "Let's resample using the `native` resampler to fix the resolutions and try regenerating our composite with the newly resampled `Scene`.\n",
    "\n",
    "Do the below step even if you didn't receive an error. It isn't necessary, but will make the remaining instructions easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_scn = scn.resample(resampler='native')\n",
    "my_rgb = comp((new_scn[r_band], new_scn[g_band], new_scn[b_band]))\n",
    "my_rgb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We successfully generated an RGB `DataArray` object. Let's plot it and see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "img = get_enhanced_image(my_rgb)\n",
    "# get DataArray out of `XRImage` object\n",
    "img_data = img.data\n",
    "img_data.plot.imshow(vmin=0, vmax=1, rgb='bands')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above image may look a little strange (really bright colors, really dull colors). This is because Satpy, unless configured otherwise, will scale RGBs between their minimum and maximum values, but this is usually not desired. We told xarray and matplotlib that our RGB's data was between 0 and 1, but depending on the channels you chose it was actually much higher than that (ex. 0-120% for reflectance, 200-250K for brightness temperatures).\n",
    "\n",
    "### Exercise:\n",
    "\n",
    "**Time: 10-15 minutes**\n",
    "\n",
    "To find the ranges that make our data look the best, let's linearly scale our data below. We'll scale each band separately where the first number represents our minimum, the second our maximum, and the third is our minimum again. The end result should be an RGB `DataArray` where \"main\" part of our data is between 0 and 1. Keep in mind what type of data we are scaling; reflectances or brightness temperatures.\n",
    "\n",
    "Remember the limits that produce the best image for you. We'll use them in the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_norm = (new_scn[r_band] - 220) / (260 - 220)\n",
    "g_norm = (new_scn[g_band] - 220) / (260 - 220)\n",
    "b_norm = (new_scn[b_band] - 0) / (120 - 0)\n",
    "my_rgb = comp((r_norm, g_norm, b_norm))\n",
    "\n",
    "plt.figure()\n",
    "my_rgb.plot.imshow(vmin=0, vmax=1, rgb='bands')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configure Composite Recipes\n",
    "\n",
    "Now that we've found an RGB recipe that we like, we may want to reuse it later. We could write a python script or module that generates this, but that requires that we know what bands go in to the RGB when loading from the Scene and can make using the composite more difficult.\n",
    "\n",
    "Satpy was originally built to help automatic and repeated operational processing. To help with this Satpy uses a series of YAML configuration files that are builtin to the package, but can also be customized by users. This tutorial comes with skeleton YAML files to be filled in for this exercise that we will fill in now.\n",
    "\n",
    "1. Add the composite recipe\n",
    "\n",
    "   To add the band combinations to a YAML configuration, first navigate to the `notebooks/composites/` directory [here](/tree/notebooks/composites) and open the `abi.yaml` file in a text editor. Or you can use Jupyter Notebook's editor by going to [abi.yaml](/edit/notebooks/composites/abi.yaml). We want to add the below text to this file:\n",
    "   \n",
    "   ```yaml\n",
    "   composites:\n",
    "     my_rgb:\n",
    "       compositor: !!python/name:satpy.composites.GenericCompositor\n",
    "       prerequisites:\n",
    "       - name: YOUR_R_CHANNEL\n",
    "       - name: YOUR_G_CHANNEL\n",
    "       - name: YOUR_B_CHANNEL\n",
    "    ```\n",
    "    \n",
    "    Since this `composites` directory is in our current directory Satpy will find it and include our composite as something we can load. We've told it our composite's name is `my_rgb`, that we want to use the `GenericCompositor` class, and that it depends on the three channels we've specified.\n",
    "    \n",
    "2. Add the enhancement recipe\n",
    "\n",
    "   To add our color limits we will edit a separate file in `notebooks/enhancements/abi.yaml`. Click [here](/edit/notebooks/enhancements/abi.yaml) to edit the file using Jupyter Notebook's editor. We want our enhancement to look like this:\n",
    "   \n",
    "   ```yaml\n",
    "   enhancements:\n",
    "     my_rgb:\n",
    "       name: my_rgb\n",
    "       operations:\n",
    "       - name: color_limits\n",
    "         method: !!python/name:satpy.enhancements.stretch\n",
    "         kwargs: {stretch: 'crude', min_stretch: [RMIN, GMIN, BMIN], max_stretch: [RMAX, GMAX, BMAX]}\n",
    "   ```\n",
    "   \n",
    "   This will use the builtin `stretch` method to do a simple linear scaling between our specified minimum and maximum values. You should replace the `XMIN` and `XMAX` values with the values you came up with in the above exercise.\n",
    "   \n",
    "\n",
    "Now that we've configured our composite for creation and scaling, we can load it like any other composite in Satpy. Note how we don't have to explicitly preload our input channels, our recipe tells the `Scene` what it needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = glob('../data/abi_l1b/20180511_texas_fire_abi_l1b_conus/*.nc')\n",
    "scn = Scene(reader='abi_l1b', filenames=filenames)\n",
    "scn.available_composite_names()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scn.load(['my_rgb'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like before, if our composite needs multiple resolutions then we need to resample to make all inputs the same size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_scn = scn.resample(resampler='native')\n",
    "new_scn['my_rgb']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We successfully generated our RGB, let's plot it. We'll use the `get_enhanced_image` function like we did at the beginning of this lesson to apply our configured enhancement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "# get DataArray out of `XRImage` object with .data\n",
    "img_data = get_enhanced_image(new_scn['my_rgb']).data\n",
    "img_data.plot.imshow(vmin=0, vmax=1, rgb='bands')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we've loaded our composite from a pre-configured YAML file and displayed it just like the `airmass` and `natural_color` RGB we started with. We could have used other methods we learned about earlier, like `crop` or `save_datasets`, to do other operations. We'll keep using this RGB in the next couple lessons and continue to show how reusing a configured recipe can be make our work easier."
   ]
  },
  {
   "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.6.7"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}