01 - Introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geologic Image Processing in Python\n",
    "===============================\n",
    "\n",
    "As a geoscientist, some of the most useful and frequently-used computational tools fall under the broad category of image processing.  It's more than working with photographs or satellite imagery, though.  All \"image processing\" means in this context is working with data that's on a regular grid. For example, a digital elevation model is every bit as much an image as a core photograph is. \n",
    "\n",
    "Outline for Today\n",
    "--------------------------\n",
    "\n",
    "  * Overview / Introduction\n",
    "  * Seamount Detection Example\n",
    "    - Thresholding\n",
    "    - Filtering\n",
    "    - Segmentation\n",
    "    - Simplification\n",
    "  * Slope and Hillshade of Topographic Data\n",
    "    - Gradients\n",
    "    - Hillshade\n",
    "  * Toe of Slope Detection\n",
    "    - Laplacian\n",
    "    - Skeltonization\n",
    "  * Lineament Analysis from Aerial Photography\n",
    "    - Edge detection\n",
    "    - Hough Transform\n",
    "    - Structure Tensor\n",
    "  * Grain Detection in Thin Sections\n",
    "    - Color-based segmentation\n",
    "    - SLIC\n",
    "    - Image moments\n",
    "   \n",
    "\n",
    "#### Goals\n",
    "\n",
    "This tutorial will introduce you to some core image processing methods by solving a handful of realistic tasks related to geology and geophysics.  The goal is to gain familiarity with key \"building blocks\" and terminology so that you can understand how to use common Python libraries such as `scipy.ndimage` and `sklearn` in your day-to-day work.  For many of you, these may seem like simple tasks and things that are trivial to accomplish in ArcGIS, ImageJ, or Photoshop.  However, the terminology is a bit different when working with image processing and computer vision libraries.  Many operations are called very different things, or are broken into smaller pieces.  Therefore, it's important to understand how to string the fundamental operations that are usually exposed in programming libraries into the higher-level operations you're used to thinking about.  \n",
    "\n",
    "There are a lot of great tutorials out there already for the libraries we'll be working with.  However, there are not many geoscience-focused examples freely available.  It's much easier to see how methods can be applied to your domain when there are examples of familiar tasks related to what you're doing.  \n",
    "\n",
    "Often, particularly with image processing, what we want to do as scientisits is significantly different than the tasks that image processing tutorials are aimed at.  The methods are the same, but thousands of examples of manipulating cat photos doesn't always help form a connection to an analysis you're stuck on.  Hopefully this tutorial can provide a bit of a bridge between the two worlds.\n",
    "\n",
    "This is not meant to be a complete introduction to image processing, or even a complete introduction to common geoscience image processing problems. Those of you already familiar with image processing methods will notice that we completely gloss over some very important points and do not fully explain many underlying methods.  The goal here, however, is to give you a quick overview of what's possible by combinging a few well-known methods.  Hopefully after this tutorial you feel comfortable enough to start experimenting and learning more on your own.\n",
    "\n",
    "#### Libraries\n",
    "\n",
    "We'll focus on using a combination of `numpy`, `scipy.ndimage`, and `sklearn`. `sklearn` is a leading image processing library that has many very nice features and exposes many advanced methods.  `scipy.ndimage` is a bit more low-level, but has the advantage of both working in 3D (or N-D) and focusing on efficient implementations of common operations.   We'll also use libraries like `rasterio` for reading and writing geospatial data, but we won't dive deeply into the details of working with geospatial data.\n",
    "\n",
    "#### Notes\n",
    "\n",
    "Because the Transform2020 tutorial is remote, it's difficult to provide the type of hands-on help we normally would.  Therefore, this set of notebooks is meant to be a \"cookbook\" demonstrating common tasks and illustrating underlying principles through specific examples.  We won't go over all of the details, but hopefully you can come back to these examples later and adjust them to your needs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Additional Resources\n",
    "\n",
    "This tutorial assumes basic Python knowledge and at least some familiarity with `numpy`.  However, if you're new to all of this, that's perfectly okay!  In that case, focus on how the operations are chained together more than the details of the code. If you're looking for tutorials on learning Python for scientific purposes, here are some resources you may find useful:\n",
    "\n",
    "Scipy-lectures is a great introduction to scientific computing in python: https://scipy-lectures.org/index.html  I can't recommend this strongly enough, honestly.  If you've picked up basic python syntax and looking for a place to start actually applying python, this is a great practical guide to the tools you'll need to know most.  It has sections on most of the libraries we'll use today, and I'd recommend browsing through it even if you're an expert.\n",
    "\n",
    "The (new) official `scipy.ndimage` tutorial is also a good overview of a library we'll be using extensively here: https://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html  However, it's most useful if you're already familiar with the basics of both image processing and numpy/scipy operations.  You may find it a bit dense otherwise, but it's full of excellent examples of all of the key functionality in `scipy.ndimage`.\n",
    "\n",
    "There are many `skimage` tutorials out there, but the gallery is a great place to start: https://scikit-image.org/docs/dev/auto_examples/index.html  Image processing operations are often visual, so it's not uncommon to suspect something is easily accomplished but not know the name of the operation.  The `skimage` gallery is something of a visual index to many of the operations in the library."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diving in: Identifying Seamounts\n",
    "------------------------------------------------\n",
    "\n",
    "Let's get started with some concrete examples.  If you'd like to see where we're going, you can jump straight to the complete example and run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "%load examples/seamount_detection.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to try to detect, count, and calculate areas of seamounts based on bathymetry data.  Along the way, we'll cover the following image processing concepts:\n",
    "\n",
    "  * Array representation\n",
    "  * Thresholding\n",
    "  * Filtering\n",
    "  * Segmentation\n",
    "  \n",
    "Let's start by loading our data from a geotiff and taking a look at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import rasterio as rio\n",
    "\n",
    "from context import data\n",
    "\n",
    "# Let's load data from a geotiff using rasterio...\n",
    "with rio.open(data.gebco.seamounts, 'r') as src:\n",
    "    bathy = src.read(1)\n",
    "    \n",
    "print(bathy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So `bathy` is a 2D array with integer values.  The units are meters relative to sea level.  Note that more or less everything is negative: This is GEBCO bathymetry data from a of the Western Pacific near the Marianas Trench.\n",
    "\n",
    "Let's take a look at what this data looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "im = ax.imshow(bathy, cmap='gray')\n",
    "fig.colorbar(im, orientation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a single array of data that we're displaying as grayscale.  Let's go ahead and add some color to that. We'll discuss color in images in more detail later, but this is a good chance to briefly introduce using colormaps in matplotlib:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "im = ax.imshow(bathy, cmap='Blues_r', vmax=0)\n",
    "im.cmap.set_over('green') # Just display any land as green...\n",
    "fig.colorbar(im, orientation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so we're looking at a large number of seamounts rising up above the abyssal plain.  It's really obvious where they are visually, but it would be nice to be able to quickly identify them programatically.  For example, we might want to look at their distribution by area or volume, or to just get a count without manually counting all of them.\n",
    "\n",
    "### Thresholding\n",
    "\n",
    "The simplest approach we could take would be to threshold the bathymetry data.  The abyssal plain is usually around 4km depth due to the relatively constant thickness and density of oceanic crust.  Therefore, we could try thresholding out anything above 3500 meters.  Because this is a numpy array, the operation is quite simple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_threshold = bathy > -3500\n",
    "print(simple_threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have an array of True/False values. This is a boolean array, and some of the operations we'll work with today only operate on these sort of boolean True/False arrays. \n",
    "\n",
    "Often, in image processing, you'll convert the True/False representation into a 1/0 representation.  Behind the scenes, the `True` values above can be efficiently converted into `1` and the `False` values into `0`. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \"view\" only changes the way we're interepting the underlying data.\n",
    "# If you're not familiar with \"view\" vs \"astype\", use \"astype\". All I'm\n",
    "# showing here is that it's seamless to go from True/False --> 1/0\n",
    "print(simple_threshold.view(np.uint8)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, enough beating around the bush. Let's take a look at what we've accomplished:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True, constrained_layout=True)\n",
    "axes[0].imshow(bathy, cmap='Blues_r', vmax=0)\n",
    "axes[1].imshow(simple_threshold.view(np.uint8))\n",
    "\n",
    "for ax in axes.flat:\n",
    "    ax.set(xticks=[], yticks=[])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The yellow regions are `True` in the boolean array. Note that we capture many seamounts out in the abyssal plain, but classify the entire volcanic arc and forearc in the west as a single large seamount.  We also miss a lot of smaller seamounts that \n",
    "\n",
    "Let's take a second to make a fancier display so we can explore what we capture and what we don't.  I'm going to use a quick utility included with this tutorial to allow toggling of different overlays on the plot.  We'll re-use this throughout the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from context import utils\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(bathy, cmap='Blues_r', vmax=0).cmap.set_over('green') \n",
    "\n",
    "# We'll mask any False values so that they're transparent\n",
    "im = ax.imshow(np.ma.masked_where(~simple_threshold, simple_threshold),\n",
    "               vmin=0, vmax=1, label='>3500 mbsl')\n",
    "\n",
    "ax.set(xticks=[], yticks=[])\n",
    "\n",
    "utils.Toggler(im).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, we're doing an okay job of detecting large seamounts, but an awful job of detecting smaller ones and arc volcanoes.  This is because we're using a fixed elevation threshold to determine whether or not something is a seamount.\n",
    "\n",
    "Visually, we'd determine whether or not a pixel is part of a seamount based on the area around it. We're looking for features that rise up from the surrounding topography.  However, that \"base level\" of topography varies throughout our study area.  Therefore we need a way of finding the \"background\" elevation.  Remember that -- we'll come back to it soon.\n",
    "\n",
    "### Filters\n",
    "\n",
    "Filters (and convolution, which is verly closely related) are an ubiquitous concept in image processing.  Filtering an image is a type of \"moving window\" operation.  For each pixel in the image, we take a region around it and apply some operation based on that region to define a new pixel value.  Most commonly used filters involve multiplying each pixel in the region by a weight and then summing (i.e. a convoluion). The array of weights is usually [called a kernel](https://en.wikipedia.org/wiki/Kernel_(image_processing)). This allows blurring, sharpening, edge detection, and many other useful operations.  \n",
    "\n",
    "Other filers aren't defined by weights, but by more flexible operations.  A simple example of this is a median filter, where the value of the pixel is the median of the pixels in some window surrounding it.  To calculate a median, we need to sort all the pixels we're using and find the one in the middle -- it can't be defined by a weighted average.  As a result, it can become slow if a large window is used.  (Note: in practice, there are some shortcuts to the \"full\" sort, but either way a median filter is slower than filters defined by weights.)\n",
    "\n",
    "Let's use one of the simplest possible filters: [a uniform filter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.uniform_filter.html).  A uniform filter is an average of all pixel values in a square region.  It's simple, but it's fast.  In practice, it blurs the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.ndimage\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "im = ax.imshow(bathy, cmap='Blues_r', vmax=0)\n",
    "\n",
    "ax.set(xticks=[], yticks=[])\n",
    "\n",
    "# Values are width of the square window in _pixels_\n",
    "def update(value):\n",
    "    blurred = scipy.ndimage.uniform_filter(bathy, value)\n",
    "    im.set_data(blurred)\n",
    "\n",
    "utils.Slider(ax, 0, 150, update, start=50).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so we mentioned this was a moving window. We need to think about what happens at the edges.  By default, most functions in `scipy.ndimage` use \"reflect\" boundary conditions.  That's perfect for this use case, but it's good to have a look at the other options.  The `mode` kwarg controls how boundaries are handled: (this is the same for any function where boundaries matter in `scipy.ndimage`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy.ndimage.uniform_filter?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using filters in seamount detection\n",
    "\n",
    "Going back to our seamount detection problem, let's use the blurred / uniform-filtered bathymetry to define the \"background\" elevation.  Seamounts or other peak-like features will be significantly higher than the background elevation and trenches or other trough-like features will be significantly below it.  Therefore, we can identify seamounts by comparing the uniform-filtered bathymetry data to the original bathymetry data.  Anything that's more than some amount higher than the background elevation, we'll consider a seamount:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# We'll use a simple filter to define the local background elevation and\n",
    "# assume anything more than 500m above the background is a seamount\n",
    "blurred = scipy.ndimage.uniform_filter(bathy, 150)\n",
    "better_threshold = bathy > (blurred + 500)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(bathy, cmap='Blues_r', vmax=0).cmap.set_over('green') \n",
    "\n",
    "# Compare this to our original 3500m threshold\n",
    "im1 = ax.imshow(np.ma.masked_where(~simple_threshold, simple_threshold),\n",
    "                vmin=0, vmax=1, label='>3500m')\n",
    "im2 = ax.imshow(np.ma.masked_where(~better_threshold, better_threshold),\n",
    "                vmin=0, vmax=1, label='Above Background')\n",
    "\n",
    "ax.set(xticks=[], yticks=[])\n",
    "fig.tight_layout()\n",
    "\n",
    "utils.Toggler(im1, im2).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleanup of Thresholded Regions\n",
    "\n",
    "As you can see, we've done a better job identifying smaller seamounts, and we're no longer identifying the entire volcanic arc as one gigantic seamount.\n",
    "\n",
    "However, we're also identifying some very tiny features that we'd rather leave out, and the boundaries of each feature are very rough.   It would be nice to \"clean up\" our classification so that we have smoother boundaries, holes are filled in, and very small regions are excluded.\n",
    "\n",
    "To do this, we'll operate on our thresholded boolean array, rather than operating on the bathymetry data directly.  \n",
    "\n",
    "First, let's fill any \"holes\" in our classified areas.  Anything that's fully surrounded by something we're considering a seamount is clearly a seamount as well.  There aren't many holes in our classification, but if you look around, you can find several small ones.  They may not matter much for this case, but it's an extremely common post-processing step in classifications, and one that's very frequently handy to know how to do.\n",
    "\n",
    "To fill holes, we'll rely on [mathematical morphology](https://en.wikipedia.org/wiki/Mathematical_morphology) (which is a cornerstone of image processing that comes from geology, by the way).  Most mathematical morphology operators work on boolean arrays similar to what we have.  It provides ways of identifying connected and non-connected regions, buffering, eroding, and similar key operations on classified images.  It's straight-forward to identify holes that are surrounded by True values in a boolean array, and `scipy.ndimage` even gives us a one-step function to do it:  [`scipy.ndimage.binary_fill_holes`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.binary_fill_holes.html#scipy.ndimage.binary_fill_holes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filling holes in our classified regions is straight-forward\n",
    "filled = scipy.ndimage.binary_fill_holes(better_threshold)\n",
    "\n",
    "# And let's compare the results:\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(bathy, cmap='Blues_r', vmax=0).cmap.set_over('green') \n",
    "\n",
    "im1 = ax.imshow(np.ma.masked_where(~better_threshold, better_threshold),\n",
    "                vmin=0, vmax=1, label='Original')\n",
    "im2 = ax.imshow(np.ma.masked_where(~filled, filled),\n",
    "                vmin=0, vmax=1, label='Filled')\n",
    "\n",
    "ax.set(xticks=[], yticks=[])\n",
    "fig.tight_layout()\n",
    "\n",
    "utils.Toggler(im1, im2).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's do something a bit more interesting.  Let's clean up our classification with a single filter.  This will smooth boundaries and eliminate very small features we're not interested in.\n",
    "\n",
    "Remember the median filter we very briefly talked about earlier?  You can also apply a median filter to boolean True/False data.  In that case, the pixel will be whatever the majority of pixels in the region around it are.  The large a region (a.k.a. kernel) we choose, the smoother/simpler the result will be.  \n",
    "\n",
    "Let's have a look at that and experiment with different filter sizes. We'll use [`scipy.ndimage.median_filter`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.median_filter.html#scipy.ndimage.median_filter) for this.\n",
    "\n",
    "Quick note: If we had more than one class (e.g. say this was a land use classification with water, urban, forest, and agriculture classes), the equivalent \"cleanup\" operator would be a [majority filter](https://scikit-image.org/docs/stable/api/skimage.filters.rank.html#majority).  For a boolean True/False array, a median filter is equivalent to a majority filter.\n",
    "\n",
    "This will take a bit to run. A median filter is a much slower operation than the other filters we've applied so far..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "widths = [3, 5, 11, 21]\n",
    "versions = ([filled] +\n",
    "            [scipy.ndimage.median_filter(filled, x) for x in widths])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(bathy, cmap='Blues_r', vmax=0)\n",
    "\n",
    "ims = []\n",
    "for width, version in zip([0] + widths, versions):\n",
    "    transparent = np.ma.masked_equal(version.view(np.uint8), 0)\n",
    "    im = ax.imshow(transparent, vmin=0, vmax=1, label=f\"{width}x{width}\")\n",
    "    ims.append(im)\n",
    "\n",
    "ax.set(xticks=[], yticks=[])\n",
    "fig.tight_layout()\n",
    "\n",
    "utils.Toggler(*ims).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's choose a final width and we'll fill holes again, as the median filter can introduce _new_ holes into the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleaned = scipy.ndimage.median_filter(better_threshold, 13)\n",
    "cleaned = scipy.ndimage.binary_fill_holes(cleaned)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so we've simplified things reasonably well, but we still have a single True/False array.  \n",
    "\n",
    "### Identifying individual features\n",
    "\n",
    "Originally we said we were going to count the seamounts and calculate areas.  How do we go from a bunch of True/False regions to a count?  The answer is to use mathematical morphology again.  We can identify and label each distinct group of `True` pixels that touch each other but do not touch any other `True` pixels.  In `scipy.ndimage`, this is the [`label`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.label.html#scipy.ndimage.label) function.  It assigns a unique value to each group of pixels that touch each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels, count = scipy.ndimage.label(cleaned)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(bathy, cmap='Blues_r', vmax=0)\n",
    "ax.imshow(np.ma.masked_equal(labels, 0), cmap='tab20b')\n",
    "\n",
    "ax.set(xticks=[], yticks=[],\n",
    "       title=f'There are {count} different seamounts!')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, great! `labels` is now an array where each separate region has a unqiue value. The colormap doesn't quite show it, but if you hover over each seamount, you'll see that the value is different.\n",
    "\n",
    "Remember we mentioned calculating area?  Let's go ahead and get a pixel count for each seamount.  We could do something like:\n",
    "\n",
    "```\n",
    "seamount_area = []\n",
    "for i in range(1, count + 1):\n",
    "    seamount_area.append((labels == i).sum())\n",
    "```\n",
    "\n",
    "However, that's a bit inefficient in this case, and `scipy.ndimage` has a builtin function to do similar \"zonal statistics\", so let's go ahead and use [`scipy.ndimage.sum`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.sum.html) to get a pixel count for each individual seamount and then plot a histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pixel_counts = scipy.ndimage.sum(cleaned, labels, np.arange(count)+1)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.hist(pixel_counts, bins=40)\n",
    "ax.set(ylabel='Number of seamounts', xlabel='Area in pixels')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a histogram, but what the heck does \"area in pixels\" mean?  Units matter!\n",
    "\n",
    "### Calculating Area (for a geographic WGS84 raster)\n",
    "\n",
    "Let's convert that to something more sensible.  We're working with data projected in geographic WGS84 in this case.  We can find the cellsize of the pixels, but it's in degrees.  Because it's in degrees, the area of the pixel varies by latitude.  \n",
    "\n",
    "Normally, we'd just multiply by a constant factor to convert pixel counts to area, but because we're working in degrees, the area varies depending on location.\n",
    "\n",
    "We can very closely approximate the true area with a simple conversion.  A degree of latitude is a constant ~111.3 km. A degree of longitude is ~111.3 km at the equator. At other latitudes, multiply by the cosine of the latitude to get the width of a degree of longitude.  If you need more precision, you can reproject into a proper projection, but this is precise enough for most purposes.\n",
    "\n",
    "Let's use `rasterio` to get the latitude and longitude of each pixel and use that to calculate the area of each pixel, then we'll sum over each seamount's region to get their total area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate area for each pixel (memory-inefficient, but fine for now)\n",
    "i, j = np.mgrid[:bathy.shape[0], :bathy.shape[1]]\n",
    "with rio.open(data.gebco.seamounts, 'r') as src:\n",
    "    cellsize = src.transform.a # Only because pixels are square an north-south\n",
    "    lon, lat = src.xy(j, i)\n",
    "area = (cellsize * 111.3)**2 * np.cos(np.radians(lat))\n",
    "\n",
    "areas = scipy.ndimage.sum(area, labels, np.arange(count)+1)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.hist(areas, bins=40)\n",
    "ax.set(ylabel='Number of seamounts', xlabel='Area in $km^2$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick aside... Note the shape of the histogram.  It's very similar to a [log-normal distribution](https://en.wikipedia.org/wiki/Log-normal_distribution)  (i.e. see the long tail on the right hand side).  This is very common in any natural process that involves areas or volumes.  As an arm-wavy explanation: If you have many independent normally distributed variables (say, width and height) and you multiply them together, the resulting distribution will be approximately log-normal.  (We don't actually expect seamount size to be as simple as width x height in this case, but still, expect log-normal distributions instead of normal anytime you start working with areas or volumes.)  Be careful about the way you interpet statistics like standard deviation or mode anytime area or volume pops up.\n",
    "\n",
    "### But I Need to Use This Data Elsewhere...\n",
    "\n",
    "I don't want to spend too long talking about geopatial data formats and I/O.  However, it's the natural next step to ask.  Therefore, I've included an example that both saves what we've done as a geotiff and vectorizes our regions and saves them as a shapefile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load examples/seamount_detection_saving.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Review\n",
    "\n",
    "Okay, let's quickly review the steps we took and have a look at a complete, stand-alone example.  To detect seamounts, we did the following:\n",
    "\n",
    "  1. Load bathymetry data into an array\n",
    "  2. Tried thresholding based on a constant elevation (didn't work well)\n",
    "  3. Estimated background elevations with a local average (uniform filter)\n",
    "  4. Detected seamounts as being regions more than 500m above the background (worked well)\n",
    "  5. Cleaned up our detection to remove small regions and have smoother boundaries (median filter)\n",
    "  6. Filled in any holes in our detected seamounts\n",
    "  7. Separated each connected region of pixels into a unique seamount (i.e. count how many)\n",
    "  8. Calculated the area distribution of seamounts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load examples/seamount_detection.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take-home Question:\n",
    "-------------------------------\n",
    "\n",
    "If you're comfortable with what we've done so far, here's a short challenge for you to try later:\n",
    "\n",
    "Some of these seamounts the surface and form islands.  We usually don't consider a feature to be a seamount if a part of it reaches the surface.  **Can you exclude all seamounts that reach the surface?**\n",
    "\n",
    "(Hint: This mostly needs numpy boolean indexing.  There's also another zonal statistics function in [`scipy.ndimage`](https://docs.scipy.org/doc/scipy/reference/ndimage.html#measurements) that you might find useful, but there's definitely more than one way to accomlish this.)\n",
    "\n",
    "If you're very familiar with numpy alread and that seems too obvious, try some of the [mathematical morpohology operations](https://docs.scipy.org/doc/scipy/reference/ndimage.html#morphology) we didn't talk about. How would you join together seamounts that are within some distance of X pixels of each other so that they're considered a single connected region? (Hint: The most efficient way will change the shape of the regions slightly. There's more than one way to approach it, though.)\n",
    "\n",
    "There are example answers in the `examples` folder for both of these if you get stuck. However, there's very much more than one right way to do it, so I'd encourage experimentation!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Break!\n",
    "----------\n",
    "\n",
    "Let's take five minutes. My voice needs a break!\n",
    "\n",
    "Next Section\n",
    "------------------\n",
    "\n",
    "We'll move on to talking about gradients in images next. We'll start out with simple slope and hillshade calculations, but then use the same methods to identify separate grains in thin section images.\n",
    "\n",
    "[02 - Gradients and Edges](./02%20-%20Gradients%20and%20Edges.ipynb)"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}