examples/Getting_Started_Extracting_Vectors_From_Images.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our previous tutorial on `Analysis with Satellite Images`, we looked at using band indices and a bit of band algebra to identify features. The next step for many people is to extract those features as vectors for further analysis.\n",
    "\n",
    "## Raster vs. Vector\n",
    "In the world of spatial data, one way of dividing up the world is raster data. Raster data represents the world as arrays with single values per array coordinate. Examples would be satellite images, digital elevation models, and other data layers that are arranged on a regular grid . Vector data, on the other hand, represent the world through geometric objects: points, line, and polygons. Examples of vector data would be building footprints, road networks, and points for store locations.\n",
    "\n",
    "Each data type has strengths and weaknesses, making them appropriate for different kinds of analysis. This  <a href=\"https://gis.stackexchange.com/questions/7077/what-are-raster-and-vector-data-in-gis-and-when-to-use\" target=\"_blank\" rel=\"noopener\">stack exchange</a> discussion is a good introduction to the topic. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Previous Analysis\n",
    "In the previous tutorial we used basic segmentation to extract water features. Let's start by redoing the analysis to get the results into this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gbdxtools import CatalogImage\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_id = '2012052117443370000011625025'\n",
    "bbox = '-104.61636543273927, 38.23971962412238, -104.58580970764162, 38.26183422148863'\n",
    "\n",
    "#No Pan sharpening or Acomp\n",
    "pueblo_image = CatalogImage(catalog_id, band_type=\"MS\", bbox=map(float, bbox.split(\",\")))\n",
    "\n",
    "#Pull out the separate bands so we can calculate a water index\n",
    "p_green = pueblo_image[1,:,:].astype(np.float32)\n",
    "p_nir = pueblo_image[3,:,:].astype(np.float32)\n",
    "\n",
    "#Calculate the water index\n",
    "nwdi = ((p_green - p_nir)/(p_green + p_nir))\n",
    "\n",
    "#Create a band (raster) that has True for all the water locations and False everywhere else\n",
    "nwdipoint3 = nwdi.copy()\n",
    "nwdipoint3 = nwdipoint3 >= 0.3\n",
    "\n",
    "#plot the color image\n",
    "fig = plt.figure(figsize=(20,20))\n",
    "a = fig.add_subplot(1,2,1)\n",
    "a.axis('off')\n",
    "a.set_title('R, G, B')\n",
    "imgplot = plt.imshow(pueblo_image.rgb(bands=[2,1,0]))\n",
    "\n",
    "a = fig.add_subplot(1,2,2)\n",
    "a.axis('off')\n",
    "a.set_title('Water')\n",
    "imgplot = plt.imshow(nwdipoint3, cmap=\"Greys_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to take the water areas in the black and white image and make them into polygons with data indicating they are water."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting Rasters To Vectors\n",
    "\n",
    "The next step is to turn our raster of water/not-water into vectors to view them in other ways, such as in GIS software. One of the best exchange formats for simple vector data is  <a href=\"https://tools.ietf.org/html/rfc7946\" target=\"_blank\" rel=\"noopener\">geojson</a>. Most GIS software and web mapping frameworks understand GeoJSON.\n",
    "\n",
    "Projections and datums define how we take a spherical 3-D object, the earth, and turn it into a flat 2-D object, like a map or a satellite image. If you ever plot your vectors on a map and they are in the wrong place, it's probably because your data is in two different projections or datums. \n",
    "\n",
    "All the satellite images you get in notebooks come unprojected with a datum of WGS84 and units of decimal degrees. You may also see this referred to as ESPG:4326.  \n",
    "\n",
    "#### Take-aways on Projections and Datums\n",
    "\n",
    "1. All DG satellite images will be in unprojected WGS84 Decimal Degrees (EPSG:4326).\n",
    "2. Therefore, all the vectors you extract from the images will be in EPSG:4326.\n",
    "3. The GeoJSON specification requires that coordinates are in EPSG:4326.\n",
    "4. It is likely that your other spatial data will not be in this coordinate system. You will probably need to reproject either the imported GeoJSON or your other map data. Here is documentation on how to reproject in  <a href=\"https://docs.qgis.org/2.8/en/docs/training_manual/vector_analysis/reproject_transform.html\" target=\"_blank\" rel=\"noopener\">QGIS</a> and  <a href=\"http://pro.arcgis.com/en/pro-app/tool-reference/data-management/an-overview-of-projections-and-transformations-toolset.htm\" target=\"_blank\" rel=\"noopener\">ArcGIS</a>\n",
    "5. It's important to understand the  <a href=\"http://www.geo.hunter.cuny.edu/~jochen/gtech201/lectures/lec6concepts/map%20coordinate%20systems/how%20to%20choose%20a%20projection.htm\" target=\"_blank\" rel=\"noopener\">properties</a> of the coordinate system you use - it is impossible to take a 3-D object and make it 2-D without  <a href=\"http://www.geography.hunter.cuny.edu/~jochen/GTECH361/lectures/lecture04/concepts/Map%20coordinate%20systems/Map%20projections%20and%20distortion.htm\" target=\"_blank\" rel=\"noopener\">introducing distortion</a>. \n",
    "\n",
    "\n",
    "### Rasterio\n",
    "\n",
    "<a href=\"https://rasterio.readthedocs.io/en/latest/index.html\" target=\"_blank\" rel=\"noopener\">Rasterio</a> is a Python library geared towards working with raster data using Python native idioms. It also knows how to work with GeoJSON.\n",
    "\n",
    "### Installing New Python libraries\n",
    "To get even more functionality with GeoJSON, we are going to install \n",
    " <a href=\"https://github.com/frewsxcv/python-geojson\" target=\"_blank\" rel=\"noopener\">python-geojson</a>. This package is not included in GBDX Notebooks by default. In the world of Python, the most common package installer is  <a href=\"https://pypi.org/project/pip\" target=\"_blank\" rel=\"noopener\">pip</a>. To install a package we normally would use the pip command on the command line. But in a notebook we can use an exclamation point (*!*) to run any command in a code window like it was on the command line. This will install the package into your Python environment so it is available for use in your notebooks. NOTE: if you reset your workspace, all your hand-installed libraries will be removed and need to be reinstalled.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the ! to execute an operating system command in the environment under the notebook\n",
    "print(\"Installing Deps...\")\n",
    "!pip install geojson rasterio -qq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import rasterio.features\n",
    "import rasterio.warp\n",
    "from geojson import Feature, Polygon, FeatureCollection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#look at the projection of the original image\n",
    "print(\"Coordinate system is \" + pueblo_image.metadata['georef']['spatialReferenceSystemCode'])\n",
    "\n",
    "\n",
    "# We need to cast our binary mask to an array and make sure it is an integer array. Trues = 1 and false = 0\n",
    "int_nwdi = np.asarray(nwdipoint3,np.uint8)\n",
    "\n",
    "features = []\n",
    "\n",
    "for geom, val in rasterio.features.shapes(int_nwdi, transform = pueblo_image.affine):\n",
    "    if val == 1.0:\n",
    "        #make a Polygon Feature for water\n",
    "        new_feature = Feature(geometry = geom, properties={\"LandCover\": \"water\"})\n",
    "    else:\n",
    "        new_feature = Feature(geometry = geom, properties={\"LandCover\": \"not_water\"})\n",
    "    features.append(new_feature)\n",
    "\n",
    "# Now put our features in a featureCollection\n",
    "feature_collection = FeatureCollection(features)\n",
    "\n",
    "feature_collection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export GeoJSON to File\n",
    "\n",
    "Now that we have our GeoJSON we want to get it out of the notebook and into a location where our GIS software can access it. We are going to write the JSON to the filesystem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "# open a file for writing\n",
    "with open('pueblo_water.json', 'w') as output:\n",
    "    json.dump(feature_collection, output)\n",
    "print(\"Done\")              "
   ]
  },
  {
   "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}