Masking_with_fiona.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import fiona\n",
    "import rasterio\n",
    "import rasterio.mask\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read in a polygon saved from Global Mapper as a .shp file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fn = 'test_shape.shp'  # broadly delimits study area\n",
    "fn = 'random_poly.shp' # smaller area\n",
    "with fiona.open(fn, \"r\") as shapefile:\n",
    "    features = [feature[\"geometry\"] for feature in shapefile]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'type': 'Polygon', 'coordinates': [[(376394.1489175722, 4625369.23452655), (376401.85702172574, 4625386.317351971), (376419.3565014257, 4625393.400474707), (376444.3557581398, 4625382.775790604), (376463.3135278147, 4625358.401515307), (376443.9391038612, 4625340.4853813285), (376422.8980627935, 4625331.52731434), (376406.44021879, 4625340.27705419), (376394.1489175722, 4625369.23452655)]]}]\n"
     ]
    }
   ],
   "source": [
    "print(features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read in the same polygon saved from Global Mapper in geoJSON format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Polygon\n",
      "[[[376394.148918, 4625369.234527], [376401.857022, 4625386.317352], [376419.356501, 4625393.400475], [376444.355758, 4625382.775791], [376463.313528, 4625358.401515], [376443.939104, 4625340.485381], [376422.898063, 4625331.527314], [376406.440219, 4625340.277054], [376394.148918, 4625369.234527]]]\n",
      "{'type': 'Polygon', 'coordinates': [[[376394.148918, 4625369.234527], [376401.857022, 4625386.317352], [376419.356501, 4625393.400475], [376444.355758, 4625382.775791], [376463.313528, 4625358.401515], [376443.939104, 4625340.485381], [376422.898063, 4625331.527314], [376406.440219, 4625340.277054], [376394.148918, 4625369.234527]]]}\n",
      " \n",
      "{'type': 'Polygon', 'coordinates': [[[376394.148918, 4625369.234527], [376401.857022, 4625386.317352], [376419.356501, 4625393.400475], [376444.355758, 4625382.775791], [376463.313528, 4625358.401515], [376443.939104, 4625340.485381], [376422.898063, 4625331.527314], [376406.440219, 4625340.277054], [376394.148918, 4625369.234527]]]}\n"
     ]
    }
   ],
   "source": [
    "with open('random_poly.geojson') as f:\n",
    "    data = json.load(f)\n",
    "\n",
    "for feature in data['features']:\n",
    "    print(feature['geometry']['type'])\n",
    "    print(feature['geometry']['coordinates'])\n",
    "    print(feature['geometry'])\n",
    "\n",
    "print(\" \")\n",
    "jsonfeatures = data['features'][0]['geometry']\n",
    "print(jsonfeatures)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note the only difference is that the geojson version has no brackets around it.\n",
    "So, in the rasterio.mask.mask command below, both `features` and `[jsonfeatures]` work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 6800, 8000)\n"
     ]
    }
   ],
   "source": [
    "# read in a raster and mask it\n",
    "# the nodata regions that are cropped out can be declared nans or assigned a value. Global Mapper uses -32767.\n",
    "with rasterio.open(\"zmean.tif\") as src:\n",
    "    out_image, out_transform = rasterio.mask.mask(src, [jsonfeatures],\n",
    "                                                        crop=False, nodata=float('nan'))\n",
    "    out_meta = src.meta.copy()\n",
    "    \n",
    "print(np.shape(out_image))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAFcCAYAAADGebf+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHBlJREFUeJzt3XuUZWV55/Hvc+pUd9EN3S22mHAT\nkBDH8QZB7WhMNI4XMnhZ6gzqmCWQpYFEg7eVAZ0VTbxhjAySOGgngsyAxFsmZCZEZDLiZTleuKiI\nRCXBC4oCdqCv1V3n7Gf+OKe6q5redXbVOVW1u/b3s9ZZdS57v/0eqha//bz73e+OzESSJDVHa7k7\nIEmSlpbhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDtJe7\nA3MKXHtYkpooieXuwkpm5S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY/pIk\nNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM\n4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEv\nSVLDGP6SJDWM4S9JUsMY/pIkVRQR7Yg4KyJiufsyDMNfkqSKrrjiiql169Zddt111xXL3ZdhGP6S\nJFUQEe2LL76Y6667jne84x0czNW/4S9JUgVXXHHF1HOe8xw2bdrECSecwMFc/UdmLncfygU17pwk\nadEktaqqI6J98sknT33mM59h48aN3HHHHZx99tl84QtfaGWtg/TArPwlSRpguurfuHEjACeeeOJB\nXf1b+UuS6qdGlf/+Vf+0g7n6t/KXJGkO+1f90w7m6t/KX5JUPzWp/Muq/mkHa/Vv5S9JUomyqn/a\nfKv/iHh9RNwWEd+KiKsjYmKkHa7Iyl+SVD81qPwHVf3Tqlb/EXEU8EXg0Zm5KyI+DlybmR8ZeecH\nsPKvmWd/503L3QVJUs/L56r6p01X/8CzKrTZBg6JiDawBvjJ0L1cACv/mjj7zjczWRRsLwomi4LJ\nTD7/qPctd7ckaXnUo/L/8r333vvkQeEP8L3vfY+TTjrp05l52oA2zwPeCewCPpOZ/2k0vZ0fw3+Z\n/cmP38ZkJtu7XbbPDP8Zz2/+txcvdzclaWnVI/xvysxTqmzb6XQYHx+/E9gy4+3Nmbl5RnsPAT4F\nnAHcD3wC+GRmXjnCblfSXup/UPDRn1+4N9xXFy0oChgbA6CTCa0W7RlLRp9y2+vYlcmObpcfPv4v\nlqvbktQ4Rbc7n+3uycxNc2z274A7M/NegIj4G+ApgOG/0n131yUcOrlvqsXeX0BR0Gm1WEtCt2AS\nmGjNmJJRFHRbLY79xmuYzOSeJ3xgKbstSY1UFNUu4a+43Q+BTRGxht6w/zOBGxfcuSEY/kuoO/VB\nNtLmURMT3LF7N9vm3LiA1uz5mJ3pJ0XB4bf8HruKgl2/8sFF6q0kaZ6V/5wy8ysR8UngZnr/S78F\n2Dz3XovD8F9C2+/fzfoNq6Hd+89+x+7dez/rzPjZyYQxaBdzT3loR7Du5nPZesqli9RjSWq2oqgY\n/hW3y8y3Am8doksjYfgvofvu3gmw9wCgk8n39+w54LbtIoAu01djtiP2Vf7TioJOBIfcdI4jAJK0\nCIpuxWH/itvVheG/RK7/6Js49K7VFN2kKJL1h0/A+DgAd01NQacX7dO/kH2/mN4BwPai4NDWgZdl\ncARAkhbHqCv/ujD8l8iWn+2iKJJud99Q/mHrV9NetWrv6/Ful/unn8fMq1z2jQBM9a8GmKU/AuAB\ngCSN1ijP+deJ4b8E3nHmyxi7LTny+HUccXSS3dw7M/TQ9as5bvVqAH7a3357BHS7ey//69l3AEC3\n++ADAGCyzms2SNJBqDM1VW27zoNOzNaa4b8EJtasZc/kJD+64wGKbvILj+i9X3STI45ODt2w7wBg\npvaDLh158AFAZ+ZIQFFwxNd/38sAJWlErPy1YN1ulzWHHcbObdv48Z1bmdzZ4Zhf3gBAa6w3vD99\nANAZWL3vOwAY7x8ctIti1noBHgBI0miM+Dr/2jD8l0D2/yimDwB+/rPerP88cf3eP5iiSNYdPrF3\nDYB2BHQ6s1b6m34+zuw/spm/xLEI2h4ASNJIWPlrwYpul87UFO3xcdYcdhiTO3dy39076EwVHHvS\nhlnbrjt8ghNnnAK4f8Yf1HT4b+//bEewfca+7em5Av3TAUd+4zX8xOWAJWnBnO2vBXvzZVfyrrNf\nAUB7fJyJNWsAeGDLJD/87v0UxfpZ2685dBWPWjcxaxEg2DcHYGx6NKA/KbAdMXuSYP/n7oNsGEqS\n6sbKX0PpTE3RGhujPT5Oa2ysfwqgdwDAHb3Jf52pfWHdGgtOXPvgSYCw3y9t/5n/0wcERcF4BCfd\n+gd897GXLM6XkqQVznP+Gsof/Y+P8a6zX7H3D6Q9voqJNWvYua3LnskOW7dMAtBqRW8hoG6yjonS\nA4BZ9jsAGJv+IywefH8ASVJ1Vv4aWlF0KbpdVq2eAHoHAGPj4+zcvoPO1Ha2P7BnVvUPcGh3NSeu\nW007gvbUFO0ZCwFNa0ewvdulf+EA7Qi6/UsAv/qo9y3BN5Oklcnw19D+y0eu5sJXvZI9uydptVp0\npqYYGxvjsA0PZee2bb1TAMDUVHfvaoBF/+Y+x61bNaut/e8I2AY6nWTtWKt3Q6B2m48/8sIl+FaS\ntHI57K+R2DM5SdHtMjY+zviqVb05AK2xvZcBTh8AdKcKim6SRfbvB1Bw/PreKYCxWUv/7rOx3e5d\nHdAquOx4g1+ShmXlr5FYu67Nzu17aPWX7u1dAtj7bGLtWiZ37Nh7AABb6fQPAqY94rBVs679n2n6\nfgAXHfX2Reu/JDWJl/ppJN74gat477kvZ3LnLgBWTUz0rgRotWi1WrNGAGaG/rSimxy1YTWddu+z\n6cWAeh8WvOGoty3VV5GkFc+1/TUyGzYewuTODpM7p9j2rztZfcgaotVi9cQE7fFVew8Adm7fs3ef\nbtE7/9+Z6q0R8IgNvUmA09X+eLfLix725mX5PpK0Ujnsr5F51ds/zNXvO5fJHR22rt7NtgcmacU4\nu+lNGpleCXD7A/fPOgDI/iWA0FsO+BfXr6bd/w0+e/2bluGbSNLKVn3C38F1V1XDf5ls2HgInfUF\nE2vbjK1qsW3LJN0Zo0u9A4B1TO7Y0Z8jELPO/0+vBfCwwyc4cu1rlulbSNLKVrny95y/qjjtlRfx\nD1e8gcMfvoZWq8VYK9i6ZZKp3V2yKCi6XSbWrKXVarFz2zbuv28XnakunamidxBQJMc/5vzl/hqS\ntKI57K+RO+2VFwFw3ZVv3Pve1i2TTO3pDfV3pvbsnQMwuXMnW7fspOgmv/vOy5alv5LUNF7nr0Xz\nnFf0VuH72H/9fWDfAcDuyd7SvO3xVZz/l5cuW/8kqams/LXoznj9BwC49Pyz2Lplku7UFK/70OXL\n3CtJai6v89eSOfdCA1+S6qDoVhz2r7hdXRj+kiSVsPKXJKlhqp7LT8/5S5K0MjjbX5KkhnG2vyRJ\nDVP1xj5dw1+SpJXByl+SpIZxtr8kSQ1j5S9JUsM421+SpIax8pckqWEMf0mSGsZhf0mSGsbKX5Kk\nhlmpl/q1lrsDkiTVVdEtKj+qiIgNEfHJiPiniLg9In51kb/CAVn5S5JUYhEq//cDn87Ml0TEKmDN\nArs2FMNfkqQS1c/5D678I2Id8OvAmQCZuQfYs/DeLZzhL0lSico39qk22/8E4F7g8oh4PHATcF5m\n7lhwBxfIc/6SJJUout3KD+CIiLhxxuPV+zXXBk4BLs3Mk4EdwPlL/JX2dkSSJB3APK/zvyczN82x\n2V3AXZn5lf7rT2L4S5JUL5XP+Vc4SMjMn0bEjyLilzPzO8AzgW8P18OFMfwlSSpRfbZ/5RX+Xgtc\n1Z/p/y/AWQvr2XAMf0mSSlS9fr9q+Gfm14FTh+jSSBj+kiSVWITKvxYMf0mSSozynH+dGP6SJJXw\nrn6SJDWMlb8kSQ1j+EuS1DAO+0uS1DCdqWr33am4tn9tGP6SJJXIqtfvG/6SJK0MrbFq97/LWOSO\njJjhL0lSidZYtVTPOLjS3/CXJKlEq1Ux/DH8JUlaEaz8JUlqmLpW/hExAZwOPA04EtgFfAv4+8y8\nbdD+hr8kSSXqWPlHxNuA5wE3AF8B7gEmgJOAC/sHBm/MzG+WtWH4S5JUoqaV/9cy820ln10UEUcA\nx87VgOEvSVKJOlb+mfn3Az6/h95oQKlqFzBKktRArVZUfiyViDg1Ij4bEVdGxDERcX1EPBARX4uI\nk6u0YeUvSVKJOlb+wH8D3gpsAL4EvD4znxURz+x/9quDGrDylySpRGssKj+W0Hhm/kNmXg1kZn6S\n3pN/pDfxbyArf0mSSrTHq9XIsbRL+09GxLOB9UBGxAsz828j4jeASvcgNvwlSSrRalVc2z8XuSOz\nnQP8KVAAzwHOjYiPAD8GXl2lAcNfkqQSlc/5L+Glfpn5DXqhP+28/qMyw1+SpBKVr/Nfwtv6RcTh\nwGvoVfqXARcATwFuB96Vmf86qA0n/EmSVKKmE/6uBNYCpwKfBX4ReA+9JX4/UqUBK39JkkrUcdgf\nODIzfysiArgrM5/ef/8LEfH1Kg0Y/pIklajjsD/QioiHAIcBh0bEcZn5/Yh4KLCqSgOGvyRJJWpa\n+b8b+Kf+87OBv4qIBB4N/HGVBgx/SZJK1LHyz8yrI+LjQGRmJyKuAZ4A/Dgz767ShuEvSVKJmlb+\nAAGcFhGP6r++Hah0vh8Mf0mSStWx8o+II+nN8r8buIXegcDpwPsi4hmZ+ZNBbRj+kiSVqGnl/y7g\n0sy8eOabEfEH9OYDvHJQA4a/JEkl6lj5A5sy88wH9yEviYjvVGnA8JckqUR71Vil7aK7pOG/a47P\ndlZpwPCXJKlETSv/9RHxogO8H8C6Kg0Y/pIklajpOf/PAc8r+ezzVRow/CVJKlHTyv/8zPzZMA0Y\n/pIklahp5f+NiLgVuBr4VGY+MN8GvKufJEklWq2o9IiKIwQjchTwZ8DTgO9GxN9GxBkRcUjVBgx/\nSZJK1PGWvpnZzczrMvMs4BjgcuCFwJ0RcVWVNgx/SZJKVK38q84NAIiIsYi4JSL+97D9y8w9wLfp\nLe+7ld7NfQbynL8kSSWqVvSt+Z3zP49eWFe6LO9AIuJY4AzgZcBa4K+BF2Tm7VX2N/wlSSox6vCP\niKOBfw+8E3jDQvoUEV+id97/E8CrM/PG+bZh+EuSVKLqcH6rd6nfERExM4g3Z+bm/Ta9GPhD4LAh\nunUB8PnMzP0/iIgnZubXBjVg+EuSVGKelf89mbmpbJuIOL2/zU0R8fSF9ikzP7dfu48GXkrvFMAD\nwKmD2jD8JUkq0R6vtrZ/xVH/pwLPj4jfAiaAdRFxZWa+Yr79iohH0Av7lwEd4BHAqZn5/Sr7O9tf\nkqQSo5ztn5kXZObRmXkcvUr9/y4w+L8EXAuMAy/JzF8BtlUNfrDylySpVOVh/6Vd3vde4Gjg4cDD\ngO8BDzr/PxfDX5KkEosV/pl5A3DDvDvU2/cFEbEeeDHwxxFxIrAhIp6UmV+t0obhL0lSicqz/Zd2\neV/66/lfBlwWEQ+nd83/xRFxTGYeM2h/w1+SpBI1HfafpX+Hv0uAS/oTAQcy/CVJKlHHyj8iLqf8\nHH8CvzOoDcNfkqQSNa38D3RPgGOB1wGVrk00/CVJKlHHyj8zPzX9PCJOAN4M/DpwIfDhKm0Y/pIk\nlahp5U9E/BvgLcDJwHuBczKzU3V/w1+SpBJ1rPwj4hP0lvD9M+D1QJfeaoEAZOaWQW0Y/pIklahp\n5f9EehP73gS8sf/edAcSOGFQA4a/JEklWmPVVsFvPfgGe4umvzzwUAx/SZJKtMcr3gInli78R8Hw\nlySpRB3P+Y+C4S9JUomanvMfmuEvSVKJulf+EfEQ4Bhm5Hlm3jxoP8NfkqQSda78I+LtwJnAP7Nv\nud8EfnPQvoa/JEklal75/0fgkZm5Z747Gv6SJJWoc+UPfAvYANwz3x0Nf0mSStQ8/N8N3BIR3wJ2\nT7+Zmc8ftKPhL0lSiZoP+18BvAe4FSjms6PhL0lSiZpX/vdl5iUL2dHwlySpREa1UK+63YjdFBHv\nBv6O2cP+XuonSdJCdSqu2V91uxE7uf9z04z3vNRPkqRhTBbVTqV3Km43Spn5jIXuW/GOBZIkNU8n\ns/JjqUTEKyKiNL8j4pER8WtztWHlL0lSic6ItxuRh9K7xO8m4CbgXmACOBH4DeA+4Py5GjD8JUkq\nUcdz/pn5/oj4C3rn9p8KPA7YBdwO/HZm/nBQG4a/JEkl6hj+AJnZBa7vP+bN8JckqURdw39Yhr8k\nSSUMf0mSGqamE/6GZvhLklSizpV/RKwGXgwcx4w8z8w/GbSv4S9JUok6hz9wDfAAvcv9dg/YdhbD\nX5KkEjUP/6Mz87kL2dHwlySpRM3D/0sR8djMvHW+Oxr+kiSVqPmEv18DzoyIO+kN+weQmfm4QTsa\n/pIklah6Y5/uMtzYBzhtoTsa/pIklag6nN9dhmH/zPxBRDweeFr/rS9k5jeq7Otd/SRJKlHHu/pN\ni4jzgKuAI/qPKyPitVX2tfKXJKlE1XP53UXtRanfAZ6cmTsAIuI9wP8D/nzQjlb+kiSVGGXlHxHH\nRMRnI+L2iLitX7kPI5h93NHtvzeQlb8kSSVGfM6/A7wxM2+OiMOAmyLi+sz89gK7dznwlYj4n/3X\nLwQ+XGVHw1+SpBKjDP/MvBu4u/98W0TcDhwFLCj8M/OiiLiB3iV/AZyVmbdU2dfwlySpxGLN9o+I\n44CTga/Mt08RsS4zt0bE4cD3+4/pzw7PzC2D2jD8JUkqMc8Jf0dExI0z3t6cmZv33zYiDgU+Bbwu\nM7cuoFsfBU6nt6b/zKOO6L8+YVADhr8kSSXmWfnfk5mb5touIsbpBf9Vmfk3C+lTZp7e/3n8QvYH\nZ/tLklRqxLP9g96EvNsz86Jh+xYR/1jlvQOx8pckqcSIz/k/Ffht4NaI+Hr/vTdn5rXz6VNETABr\ngI0R8RD2Xd63DjiyShuGvyRJJUY82/+LVLwOf4DfBV5HL+hvmtHmVuADVRow/CVJKjFZw7X9M/P9\nwPsj4rWZOXA1vwMx/CVJKlG18i+Wfml/MvPPI+IxwKOBiRnv//dB+xr+kiSVqB7+y3Jjn7cCT6cX\n/tfSu8XvF4GB4e9sf0mSStT5rn7AS4BnAj/NzLOAxwOrq+xo5S9JUomqi/wUi9qLUrsys4iITkSs\nA+6hwgI/YPhLklSqzsP+wI0RsQH4S3qz/rcDX62yo+EvSVKJOod/Zv5e/+kHI+LTwLrM/GaVfQ1/\nSZJK1Dn8I+Ia4GPANZn5/fns64Q/SZJK1HzC30X0buf77Yj4RES8pL/630BW/pIklajzhL/M/Bzw\nuYgYA34TeBVwGb1lfudk+EuSVKLOw/4AEXEI8DzgDOAU4Ioq+xn+kiSVqHP4R8THgCcDn6a3pv8N\nmVlpEMLwlySpxGRRbUC/qLjdiF0OvDwzu/Pd0Ql/kiSV6MzjsQw+D1wQEZsBIuKXIuL0Kjsa/pIk\nlaj5bP/LgT3AU/qv7wLeUWVHw1+SpBI1D/9HZuafAlMAmbkLiCo7es5fkqQSVUM9lyf89/Rn+ydA\nRDwS2F1lR8NfkqQSNQ//t9Kb6X9MRFwFPBU4s8qOhr8kSSWqTuRbjujPzOsj4mZgE73h/vMy874q\n+xr+kiSVqGPlHxGn7PfW3f2fx0bEsZl586A2DH9JkkrUMfyB983xWdJb6ndOhr8kSSXqGP6Z+Yxh\n2/BSP0mSStTxUr+I+MMZz//Dfp+9q0obhr8kSSVqusLfS2c8v2C/z55bpQGH/SVJKlHHYX9mL+Sz\n/6I+LvIjSdIwqt7YJ5f2xj5Z8vxArw/I8JckqURNK//HR8RWelX+If3n9F9PVGnA8JckqUQdF/nJ\nzLFh2zD8JUkqUdPKf2iGvyRJJQx/SZIaxvCXJKlhMscrbriks/2HZvhLklSmWFVxO8NfkqSVwfCX\nJKlhDH9Jkhqmcvh3F7cfI2b4S5JUxvCXJKlhqoZ/Gv6SJK0MnUOrbdetFv4R8Vzg/cAY8FeZeeFC\nuzYMw1+SpDIjHPaPiDHgA8CzgLuAr0XE32Xmt4fo4YIY/pIklelWDP9qlf+TgDsy818AIuKvgRcA\nhr8kSbWRIz3nfxTwoxmv7wKePO8+jYDhL0lSmfkN+x8RETfOeHdzZm6e8ToOsOey3BTA8Jckqcz8\nwv+ezNw0x1Z3AcfMeH008JOFdm0Yhr8kSWVGe53/14BfiojjgR8DLwVevtCuDcPwlySpTOXw7wzc\nJDM7EfEa4Dp6l/pdlpm3DdO9hTL8JUkqUzn8W5U2y8xrgWsX3qHRMPwlSSoz4vCvC8NfkqQyla/z\nN/wlSVoZKl/nb/hLkrQyVF7bf/CEvzox/CVJKuOwvyRJDeOEP0mSGsZz/pIkNUx3rOJ2y7JE/4IZ\n/pIklWgVRaXtsiiW5w49C2T4S5JUIrqV1uyHbtfwlyRpJZhP+B9MDH9JkkrMZ9j/YIp/w1+SpBJW\n/pIkNYzhL0lSw8xn2P9gYvhLklTCyl+SpIYZm5qqtF0W3thHkqQVoXLlX1j5S5K0Ihj+kiQ1jBP+\nJElqGCt/SZIaxvCXJKlhHPaXJKlhrPwlSWoYw1+SpIZx2F+SpIax8pckqWEMf0mSGsa1/SVJapjq\nlb/n/CVJWhGc8CdJUsN4zl+SpIYx/CVJahiH/SVJahgrf0mSGsbwlySpYRz2lySpYaz8JUlqGMNf\nkqSGqTzsn8MN+0fEe4HnAXuAfwbOysz7h2p0Dq3FaliSpINddLuVH0O6HnhMZj4O+C5wwdCdn4OV\nvyRJJSrf2CeHC//M/MyMl18GXjJUgwMY/pIkHdgYnQ4RMXDDLPaM8t89G/jYKBvcX2TmYrY/nKDG\nnZMkLZpkcOIusoi49Bo45/kVtv0QcA78L+DIGW9vzszNM9r7P8AvHGD3t2TmNf1t3gKcCrwoFzGg\nDX9JUv3UI/yP3QQ/+BLM2Zk9wBOBb8JDhpmkFxGvBM4BnpmZOxfaThVO+JMk6QAy84dPoFfOz+Vy\n4MW97YcJ/ucC/xl4/mIHP1j5S5LqqAaVPwyu/kdY9d8BrAZ+3n/ry5l5zkLbG/jvGf6SpNqpSfgD\nnBuRpwEHOvf/IeBnwB9l1qa/VRj+kqT6qVH4l1X/o6r6l4Pn/CVJmkPZuf9RnOtfLlb+kqT6qVHl\nDw+u/g/mqh+s/CVJGmj/6v9grvrByl+SVEc1q/xhX/V/A/AkDt6qH6z8JUmqZLr6P4ODu+oHK39J\nUh3VsPKHXvX/WPjBrQdx1Q+GvySpjmoa/iuFw/6SJDWM4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEv\nSVLDGP6SJDWM4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY/pIkNYzhL0lS\nwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsMY\n/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDVMe7k7MKcklrsLkiStNFb+kiQ1jOEvSVLDGP6SJDWM\n4S9JUsMY/pIkNYzhL0lSwxj+kiQ1jOEvSVLDGP6SJDWM4S9JUsP8f7KLHtx7YEisAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2d8570da198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Change palette to highlight missing values\n",
    "import copy\n",
    "\n",
    "# Change palette to highlight missing values\n",
    "tcm = copy.copy(plt.cm.terrain) # work on a copy to prevent munting the original\n",
    "tcm.norm_clip=False # need this to make next line work\n",
    "tcm.set_bad(color='magenta',alpha=None)\n",
    "tcm.set_under('red',alpha=None)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,6))\n",
    "plt.axis('off');\n",
    "im=plt.imshow(out_image[0][:][:],cmap=tcm,vmin=-2,vmax=8)\n",
    "\n",
    "cbar = fig.colorbar(im, extend='both', shrink=0.8, ax=ax)\n",
    "cbar.set_label('Elevation (m NAVD88)')\n",
    "for ticklabel in ax.xaxis.get_ticklabels():\n",
    "    ticklabel.set_visible(False)"
   ]
  },
  {
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}