Masking_array_with_shapely.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masking a raster array with shapely\n",
    "\n",
    "Although you can use `fiona` to mask a raster as you read it in with `rasterio`, I wanted to use the mask on arrays that have already been loaded. This works with a regular grid and a polygon that is completely enclosed in the raster.\n",
    "\n",
    "TODO - Test with polygons that extend beyond the raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon, Point\n",
    "import shapely.vectorized\n",
    "import json\n",
    "# import fiona\n",
    "import numpy as np\n",
    "import rasterio # rasterio and OGR don't mix...see rasterio docs\n",
    "import rasterio.mask\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use `rasterio` to read in a geotiff\n",
    "\n",
    "Use `src.transform` to determine the coordinates of the corner points.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "RasterioIOError",
     "evalue": "E:\\2017_Karen_Sandwich_maps\\2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif: No such file or directory",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mCPLE_OpenFailedError\u001b[0m                      Traceback (most recent call last)",
      "\u001b[1;32mrasterio\\_base.pyx\u001b[0m in \u001b[0;36mrasterio._base.DatasetBase.__init__\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mrasterio\\_shim.pyx\u001b[0m in \u001b[0;36mrasterio._shim.open_dataset\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mrasterio\\_err.pyx\u001b[0m in \u001b[0;36mrasterio._err.exc_wrap_pointer\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mCPLE_OpenFailedError\u001b[0m: E:\\2017_Karen_Sandwich_maps\\2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif: No such file or directory",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[1;31mRasterioIOError\u001b[0m                           Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-3cdbb46752eb>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mfnames\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;34m\"2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;31m# read band 1 (elevation...only band in these files) into specifed window\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[1;32mwith\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfdir\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mfnames\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0msrc\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m     \u001b[0mraster\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msrc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\envs\\IOOS\\lib\\site-packages\\rasterio\\env.py\u001b[0m in \u001b[0;36mwrapper\u001b[1;34m(*args, **kwds)\u001b[0m\n\u001b[0;32m    419\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    420\u001b[0m         \u001b[1;32mwith\u001b[0m \u001b[0menv_ctor\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msession\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msession\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 421\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    422\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    423\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\envs\\IOOS\\lib\\site-packages\\rasterio\\__init__.py\u001b[0m in \u001b[0;36mopen\u001b[1;34m(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)\u001b[0m\n\u001b[0;32m    214\u001b[0m         \u001b[1;31m# None.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    215\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mmode\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'r'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 216\u001b[1;33m             \u001b[0ms\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mDatasetReader\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdriver\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdriver\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    217\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0mmode\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'r+'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    218\u001b[0m             \u001b[0ms\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget_writer_for_path\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdriver\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdriver\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mrasterio\\_base.pyx\u001b[0m in \u001b[0;36mrasterio._base.DatasetBase.__init__\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mRasterioIOError\u001b[0m: E:\\2017_Karen_Sandwich_maps\\2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif: No such file or directory"
     ]
    }
   ],
   "source": [
    "# read in a .tif file...this will result in a 32-bit array\n",
    "fdir = \"E:\\\\2017_Karen_Sandwich_maps\\\\\"\n",
    "fnames = (\"2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif\")\n",
    "# read band 1 (elevation...only band in these files) into specifed window\n",
    "with rasterio.open(fdir+fnames) as src:\n",
    "    raster = src.read(1,)\n",
    "\n",
    "    # upper left corner\n",
    "    x1,y1=src.transform*(0,0)\n",
    "    # lower right corner\n",
    "    x2,y2 = src.transform*(src.width, src.height)\n",
    "    nx, ny = (src.width, src.height)\n",
    "\n",
    "print(x1,y1)\n",
    "print(x2,y2)\n",
    "\n",
    "# Set missing values to NaN\n",
    "raster[raster<=-99]=np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10226, 11987)\n",
      "(10226, 11987)\n"
     ]
    }
   ],
   "source": [
    "print(np.shape(raster))\n",
    "# next step could also be done with np.mgrid, but not quite as intuitively\n",
    "X,Y = np.meshgrid(np.linspace(x1,x2,nx),\\\n",
    "                  np.linspace(y1,y2,ny))\n",
    "\n",
    "print(np.shape(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read a polygon that was saved in geojson format and build a shapely polygon\n",
    "\n",
    "This seems harder than it should be. This version assumes that there is only one `features` in the .geojson file. If you have multiple `features`, you have to do something like this:\n",
    "\n",
    "```\n",
    "for feature in data['features']:\n",
    "    print(feature['geometry']['type'])\n",
    "    #print(feature['geometry']['coordinates'])\n",
    "    #print(feature['geometry'])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11, 2)\n",
      "area of polygon:  5986.593737889501\n"
     ]
    }
   ],
   "source": [
    "with open('owash_poly.geojson') as f:\n",
    "    data = json.load(f)\n",
    "\n",
    "# convert to array\n",
    "a=np.squeeze(data['features'][0]['geometry']['coordinates'])\n",
    "print(np.shape(a))\n",
    "# shapely.geometry.Polygon requires a sequence of coordinate tuples\n",
    "#(not an array)\n",
    "poly = Polygon(tuple(map(tuple,a)))\n",
    "# print(poly.wkt)\n",
    "print(\"area of polygon: \",poly.area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use shapely to make a boolean mask\n",
    "\n",
    "This takes a while...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = shapely.vectorized.contains(poly,X,Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Various options for applying the mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-49-afbb248f0ec7>, line 6)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-49-afbb248f0ec7>\"\u001b[1;36m, line \u001b[1;32m6\u001b[0m\n\u001b[1;33m    mm[mm~=1]=np.nan\u001b[0m\n\u001b[1;37m         ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "print(np.shape(m))\n",
    "print(np.min(m),np.max(m))\n",
    "# this stays 32-bit\n",
    "mras1 = m*raster\n",
    "mm = (1*m)\n",
    "mm[mm~=1]=np.nan\n",
    "mras2 = (1*m)*raster\n",
    "print(np.nansum(mras1),np.nansum(mras2))\n",
    "print(np.nanmax(mras1),np.nanmax(mras2))\n",
    "print(mras1.dtype,mras2.dtype,raster.dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nan"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nan*False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}