workshop_part_4_segmentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 4: Image segmentation\n",
    "\n",
    "\n",
    "In this example we will show segmentation of a Landsat scene.\n",
    "We show two segmentation modules: [i.segment](https://grass.osgeo.org/grass-stable/manuals/i.segment.html) and the addon [i.superpixels.slic](https://grass.osgeo.org/grass-stable/manuals/addons/i.superpixels.slic.html).\n",
    "Note that each segmentation algorithm is designed for different purpose, so we can't directly compare them.\n",
    "\n",
    "First, let's create a new mapset *segmentation* in nc_spm_08_grass7 location:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/nc_spm_08_grass7/segmentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize GRASS session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"nc_spm_08_grass7\", \"segmentation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install the addon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.extension\", extension=\"i.superpixels.slic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation\n",
    "Imagery modules typically work with *imagery groups*. We first list the landsat raster data and then create an imagery group:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps = gs.list_grouped(type=\"raster\", pattern=\"lsat*\")[\"PERMANENT\"]\n",
    "print(maps)\n",
    "gs.run_command(\"i.group\", group=\"landsat\", subgroup=\"landsat\", input=maps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we derive NDVI to see some of the effects of segmentation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"lsat7_2002_30\")\n",
    "gs.run_command(\"i.vi\", red=\"lsat7_2002_30\",  nir=\"lsat7_2002_40\", output=\"ndvi\", viname=\"ndvi\")\n",
    "\n",
    "ndvi_map = gj.Map()\n",
    "ndvi_map.d_rast(map=\"ndvi\")\n",
    "ndvi_map.d_legend(raster=\"ndvi\", at=(2, 50, 1, 5))\n",
    "ndvi_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Segmentation with i.superpixels.slic\n",
    "Superpixels can be defined as a group of pixels that share common characteristics and are useful in computer vision and image processing.\n",
    "Here we run [i.superpixels.slic](https://grass.osgeo.org/grass-stable/manuals/addons/i.superpixels.slic.html) and convert the resulting raster to vector for better visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"i.superpixels.slic\", input=\"landsat\", output=\"superpixels\", num_pixels=1000, compactness=0.5)\n",
    "gs.run_command(\"r.to.vect\", input=\"superpixels\", output=\"superpixels\", type=\"area\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can play with *compactness* and *num_pixels* parameters and see how the resulting segmentation changes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "superpixels_map = gj.Map()\n",
    "superpixels_map.d_rast(map=\"ndvi\")\n",
    "superpixels_map.d_vect(map=\"superpixels\", width=1, color=\"black\", fill_color=\"none\")\n",
    "superpixels_map.d_legend(raster=\"ndvi\", at=(2, 50, 1, 5))\n",
    "superpixels_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For fun, let's do zonal statistics on the results. We compute the median NDVI value within each segment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.stats.quantile\", base=\"superpixels\", cover=\"ndvi\", output=\"superpixels_ndvi\")\n",
    "\n",
    "superpixels_map = gj.Map()\n",
    "superpixels_map.d_rast(map=\"superpixels_ndvi\")\n",
    "superpixels_map.d_legend(raster=\"superpixels_ndvi\", at=(2, 50, 1, 5))\n",
    "superpixels_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Segmentation with i.segment\n",
    "\n",
    "Next, we do the same, but with i.segment to see the different behavior. Note that i.segment uses *region growing* algorithm by default, but *mean shift* is also available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"i.segment\", group=\"landsat\", output=\"segments\", threshold=0.5, minsize=50)\n",
    "gs.run_command(\"r.to.vect\", input=\"segments\", output=\"segments\", type=\"area\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "segments_map = gj.Map()\n",
    "segments_map.d_rast(map=\"ndvi\")\n",
    "segments_map.d_vect(map=\"segments\", width=1, color=\"black\", fill_color=\"none\")\n",
    "segments_map.d_legend(raster=\"ndvi\", at=(2, 50, 1, 5))\n",
    "segments_map.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.stats.quantile\", base=\"segments\", cover=\"ndvi\", output=\"segments_ndvi\")\n",
    "\n",
    "segments_map = gj.Map()\n",
    "segments_map.d_rast(map=\"segments_ndvi\")\n",
    "segments_map.d_legend(raster=\"superpixels_ndvi\", at=(2, 50, 1, 5))\n",
    "segments_map.show()"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}