add rasterio filtering and simplifying example

Author: jreiberkyle

jupyter-notebooks/analytics-snippets/roads_as_vector.ipynb

@@ -6,7 +6,12 @@
    "source": [
     "# Roads as Vectors\n",
     "\n",
-    "This notebook demonstrates converting the roads raster that is the output of the Analaytics feed into a vector dataset."
+    "This notebook demonstrates converting the roads raster that is the output of the Analaytics feed into a vector dataset.\n",
+    "\n",
+    "It demonstrates the following techniques for converting to vector:\n",
+    "1. GDAL CLI\n",
+    "2. Rasterio (no processing)\n",
+    "3. Rasterio (with filtering and simplification)"
    ]
   },
   {
@@ -23,7 +28,9 @@
     "import rasterio\n",
     "from rasterio import features as rfeatures\n",
     "from rasterio.enums import Resampling\n",
-    "from rasterio.plot import show"
+    "from rasterio.plot import show\n",
+    "import shapely\n",
+    "from shapely.geometry import shape as sshape"
    ]
   },
   {
@@ -42,7 +49,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Identify road feed feature for download\n",
+    "## Obtain Analytics Raster\n",
+    "\n",
+    "### Identify road feed feature for download\n",
     "\n",
     "We want to download the most recent feature from the feed for road detection in Kirazli, Turkey."
    ]
@@ -100,7 +109,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Download Quad Raster"
+    "### Download Quad Raster"
    ]
   },
   {
@@ -212,7 +221,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Visualize Roads Image\n",
+    "### Visualize Roads Image\n",
     "\n",
     "The output of the analytics road detection is a boolean image where road pixels are given a value of True and non-road pixels are given a value of False."
    ]
@@ -237,7 +246,7 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.axes._subplots.AxesSubplot at 0x7fd86fae4690>"
+       "<matplotlib.axes._subplots.AxesSubplot at 0x7f64d938b550>"
       ]
      },
      "execution_count": 10,
@@ -276,7 +285,45 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Convert to Vector using Rasterio (no filtering)\n",
+    "## Convert Roads to Vector Features\n",
+    "\n",
+    "### GDAL Command-Line Interface (CLI)\n",
+    "\n",
+    "GDAL provides a python script that can be run via the CLI. It is quite easy to run and fast, though it doesn't allow for some of the convenient pixel-space filtering and processing that rasterio provides and we will use later on."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdal_output_filename = os.path.join('data', 'test_gdal.shp')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Creating output data/test_gdal.shp of format ESRI Shapefile.\n",
+      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
+     ]
+    }
+   ],
+   "source": [
+    "!gdal_polygonize.py $filename $gdal_output_filename"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Rasterio - no filtering\n",
     "\n",
     "In this section we use rasterio to convert the binary roads raster into a vector dataset. The vectors are written to disk as a shapefile. The shapefile can be imported into geospatial programs such as QGIS or ArcGIS for visualization and further processing.\n",
     "\n",
@@ -285,7 +332,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
@@ -329,57 +376,67 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Convert to Vectors with GDAL CLI\n",
+    "### Rasterio - Filtering and Simplifying\n",
     "\n",
-    "GDAL provides a python script that can be run via the CLI. It is quite easy to run and fast, though it doesn't allow for some of the convenient pixel-space filtering and processing that rasterio provides and we will use later on."
+    "In this section, we use `shapely` to filter the road vectors by size and simplify them so we don't have a million pixel edges."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 14,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "'data/1176-1272_2019-07-01.tif'"
-      ]
-     },
-     "execution_count": 27,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "gdal_output_filename = os.path.join('data', 'test_gdal.shp')"
+    "def roads_as_vectors_with_filtering(filename, min_pixel_size=5): \n",
+    "    with rasterio.open(filename) as dataset:\n",
+    "        roads = dataset.read(1)\n",
+    "        road_mask = roads == 255 # mask non-road pixels\n",
+    "\n",
+    "        # we skip transform on vectorization so we can perform filtering in pixel space\n",
+    "        road_shapes = rfeatures.shapes(roads, mask=road_mask, connectivity=8)\n",
+    "        road_geometries = (r for r, _ in road_shapes)\n",
+    "        geo_shapes = (sshape(g) for g in road_geometries)\n",
+    "\n",
+    "        # filter to shapes bigger than min_pixel_size\n",
+    "        geo_shapes = (s for s in geo_shapes if s.area > min_pixel_size)\n",
+    "    \n",
+    "        # simplify so we don't have a million pixel edge points\n",
+    "        tolerance = 1 #1.5\n",
+    "        geo_shapes = (g.simplify(tolerance, preserve_topology=False)\n",
+    "                      for g in geo_shapes)\n",
+    "\n",
+    "        # apply image transform    \n",
+    "        # rasterio transform: (a, b, c, d, e, f, 0, 0, 1), c and f are offsets\n",
+    "        # shapely: a b d e c/xoff f/yoff\n",
+    "        d = dataset.transform\n",
+    "        shapely_transform = [d[0], d[1], d[3], d[4], d[2], d[5]]\n",
+    "        proj_shapes = (shapely.affinity.affine_transform(g, shapely_transform)\n",
+    "                       for g in geo_shapes)\n",
+    "        \n",
+    "        road_geometries = (shapely.geometry.mapping(s) for s in proj_shapes)\n",
+    "        \n",
+    "        crs = dataset.crs\n",
+    "    return (road_geometries, crs)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "\r\n",
-      "gdal_polygonize [-8] [-nomask] [-mask filename] raster_file [-b band|mask]\r\n",
-      "                [-q] [-f ogr_format] out_file [layer] [fieldname]\r\n",
-      "\r\n"
+      "wrote 838 geometries to data/test_filt.shp\n"
      ]
     }
    ],
    "source": [
-    "!gdal_polygonize.py $filename $gdal_output_filename"
+    "road_geometries_filt, crs = roads_as_vectors_with_filtering(filename)\n",
+    "output_filename = os.path.join('data', 'test_filt.shp')\n",
+    "save_as_shapefile(output_filename, road_geometries_filt, crs)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {