opencv_plus_ry.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\envs\\wradlib0_11_3\\lib\\site-packages\\h5py\\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "%matplotlib\n",
    "import wradlib\n",
    "import cv2\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "import matplotlib.patches as mpatches\n",
    "from matplotlib.collections import PatchCollection\n",
    "from scipy.ndimage import zoom\n",
    "import datetime\n",
    "import warnings\n",
    "warnings.simplefilter('once', DeprecationWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read sample data\n",
    "\n",
    "Data is from the German Weather Service: the so called RY product represents rainfall intensity composite for the whole of Germany in 5 minute intervals. \n",
    "\n",
    "Spatial resolution: `1 x 1 km`; spatial extent: `900 x 900 km`.\n",
    "\n",
    "**Information required from user**\n",
    "\n",
    "- specify the directory `datadir` where you store the RY data (unpack the ry archives there).\n",
    "- select a specific interval by commenting/uncommenting the `dtimes` lines.\n",
    "- decide whether you need to reduce the resolution (downsize the image by a `downsizeby`) in order to avoid memory problems (this becomes relevant once you solve the 2D-adveciton equation...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set data directory\n",
    "datadir = \"data/ry\"\n",
    "\n",
    "# Original grid dimensions\n",
    "nx = 900\n",
    "ny = 900\n",
    "\n",
    "# pixel size (in meters)\n",
    "dx = 1000.\n",
    "dy = 1000.\n",
    "\n",
    "# Downsize by factor \"downsizeby\"\n",
    "#    downsizeby = 1 will leave the dimensions unchanged,\n",
    "#    but for a 900x900 km grid, downsizing might be \n",
    "#    required in order to avoid MemoryError\n",
    "downsizeby = 1\n",
    "\n",
    "# interval between observations (in seconds)\n",
    "interval = 300\n",
    "\n",
    "# Set time window\n",
    "##dtimes = wradlib.util.from_to(\"2008-06-02 17:00:00\", \"2008-06-02 19:00:00\", interval)\n",
    "##dtimes = wradlib.util.from_to(\"2015-04-26 17:00:00\", \"2015-04-26 19:00:00\", interval)\n",
    "##dtimes = wradlib.util.from_to(\"2015-03-29 17:00:00\", \"2015-03-29 19:00:00\", interval)\n",
    "#dtimes = wradlib.util.from_to(\"2016-05-29 16:00:00\", \"2016-05-29 19:00:00\", interval)\n",
    "dtimes = wradlib.util.from_to(\"2016-05-23 04:00:00\", \"2016-05-23 08:00:00\", interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\src\\git\\heistermann\\wradlib\\wradlib\\trafo.py:127: RuntimeWarning: divide by zero encountered in log10\n",
      "  return 10. * np.log10(x)\n"
     ]
    }
   ],
   "source": [
    "# Compute grid dimensions and grid coordinates after resampling\n",
    "dx2, dy2 = dx*downsizeby, dy*downsizeby\n",
    "nx2, ny2 = int(nx/downsizeby), int(ny/downsizeby)\n",
    "\n",
    "X2, Y2 = np.meshgrid( np.arange(0,nx2*dx2, dx2), np.arange(0,ny2*dy2, dy2) )\n",
    "\n",
    "# Define container\n",
    "frames = np.zeros( (len(dtimes), nx2, ny2 ) )\n",
    "\n",
    "# Read the data, convert back to dBZ, and downsize\n",
    "#   (maybe also try with keeping mm/h instead of converting to dBZ?)\n",
    "for i, dtime in enumerate(dtimes):\n",
    "    fname = dtime.strftime( os.path.join(datadir, \"raa01-ry_10000-%y%m%d%H%M-dwd---bin\") )\n",
    "    frames[i] = zoom( wradlib.io.read_RADOLAN_composite(fname, missing=0)[0], 1./downsizeby, order=1)\n",
    "    frames[i] = wradlib.trafo.decibel( wradlib.zr.r2z(frames[i]) )\n",
    "    frames[i][frames[i]<0] = 0 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use OpenCV's Optical Flow to detect and track features\n",
    "\n",
    "This example uses the Lucas-Kanade Optical Flow implementation in OpenCV (see [here](http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_video/py_lucas_kanade/py_lucas_kanade.html)). We take the first frame, detect some Shi-Tomasi corner points in it, then we iteratively track those points over the subsequent images.\n",
    "\n",
    "The parameter dictionaries are certainly something to experiment with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FEATURE DETECTION: Parameters for ShiTomasi corner detection\n",
    "feature_params = dict( maxCorners = 200,\n",
    "                       qualityLevel = 0.2,\n",
    "                       minDistance = 7,\n",
    "                       blockSize = 21 )\n",
    "\n",
    "# FEATURE TRACKING: Parameters for Lucas Kanade (lk) Optical Flow technique\n",
    "lk_params = dict( winSize  = (20,20),\n",
    "                  maxLevel = 2,\n",
    "                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0))\n",
    "\n",
    "# Over which time steps (of the data we've read in) do you want to track\n",
    "trackstart = 0\n",
    "trackend = len(frames)-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "49"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(frames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Our approach requires 8 bit integers - so we need to normalize our radar data accordingly\n",
    "#   (there might be a more elegant solution...)\n",
    "minval = 0\n",
    "maxval = 59 # dBZ in this case\n",
    "iframes = frames.copy()\n",
    "iframes[iframes<minval] = minval\n",
    "iframes[iframes>maxval] = maxval\n",
    "iframes = ((iframes / maxval)*255).astype(np.uint8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 168 good features to track.\n"
     ]
    }
   ],
   "source": [
    "# Find good features to track...\n",
    "old = cv2.goodFeaturesToTrack(iframes[trackstart], mask = None, **feature_params)\n",
    "print(\"Found %d good features to track.\" % len(old) )\n",
    "\n",
    "# Set containers to collect results (time steps in rows, detected corners in columns)\n",
    "#   Tracking status\n",
    "sts = np.zeros((trackend,len(old)), dtype=np.bool)\n",
    "#   corner x coords\n",
    "x = np.zeros((trackend,len(old))) * np.nan\n",
    "#   corner y coords\n",
    "y = np.zeros((trackend,len(old))) * np.nan\n",
    "#   tracking error\n",
    "errs = np.zeros((trackend,len(old))) * np.nan\n",
    "#   Assign persistent corner IDs\n",
    "ids = np.arange(len(old))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Track good features\n",
    "for i in range(trackstart, trackend):\n",
    "    # track current corners in next image\n",
    "    new, st, err = cv2.calcOpticalFlowPyrLK(prevImg=iframes[i], nextImg=iframes[i+1], prevPts=old, nextPts=None, **lk_params)\n",
    "    success = st.ravel()==1\n",
    "    ids = ids[success]\n",
    "    sts[i, ids] = True\n",
    "    x[i, ids] = old[success,0,0]\n",
    "    y[i, ids] = old[success,0,1]\n",
    "    errs[i, ids] = err.ravel()[success]\n",
    "    # new corners will be old in the next loop\n",
    "    old = new[success]\n",
    "\n",
    "# Incremental euclidic distance from starting point\n",
    "trackdist = np.diff( np.sqrt( (x-x[0].reshape((1,-1)))**2 + (y-y[0].reshape((1,-1)))**2 ), axis=0 )\n",
    "trackdist = np.vstack( (np.zeros((1,trackdist.shape[1])), trackdist))\n",
    "\n",
    "# Plot feature persistence\n",
    "fig = plt.figure(figsize=(12,6))\n",
    "ax = fig.add_subplot(211)\n",
    "cb = plt.imshow(errs, interpolation=\"none\", cmap=\"summer\", vmax = 15)\n",
    "plt.xlabel(\"Feature ID\")\n",
    "plt.ylabel(\"Tracking time step\")\n",
    "plt.colorbar(cb, shrink=0.5)\n",
    "plt.title(\"Tracking error\")\n",
    "\n",
    "# Plot consistence of movement\n",
    "ax = fig.add_subplot(212)\n",
    "cb = plt.imshow(trackdist, interpolation=\"none\", cmap=\"bwr\", vmin=-5, vmax=5)\n",
    "plt.xlabel(\"Feature ID\")\n",
    "plt.ylabel(\"Tracking time step\")\n",
    "plt.colorbar(cb, shrink=0.75)\n",
    "plt.title(\"Incremental euclidian distance from starting point\")\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 143 good tracks and 25 bad tracks.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\envs\\wradlib0_11_3\\lib\\site-packages\\ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in less\n",
      "  if __name__ == '__main__':\n",
      "C:\\Anaconda3\\envs\\wradlib0_11_3\\lib\\site-packages\\ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in greater\n",
      "  if sys.path[0] == '':\n"
     ]
    }
   ],
   "source": [
    "# Find good tracks (but what is a \"good\" track...?)\n",
    "#   Certainly a lot of subjective criteria to play with...\n",
    "goodtrack = np.zeros(x.shape[1], dtype=np.bool)\n",
    "for i in range(len(goodtrack)):\n",
    "    # persistence of the track\n",
    "    if len(np.where(sts[:,i])[0]) < 2:\n",
    "        continue\n",
    "    # consistency of movement\n",
    "    if len(np.where(trackdist[:,i]<0)[0]) > 0:\n",
    "        continue\n",
    "    # tracking error\n",
    "    if len(np.where(errs[:,i]>15)[0]) > 5:\n",
    "        continue\n",
    "    goodtrack[i] = True\n",
    "print(\"Found %d good tracks and %d bad tracks.\" % \\\n",
    "      (len(np.where(goodtrack)[0]), len(goodtrack)-len(np.where(goodtrack)[0])) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize tracks: green=good track, red=bad track\n",
    "goodcolor = \"limegreen\"\n",
    "badcolor = \"red\"\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "# average reflectivity over entire tracking period as background image\n",
    "ax.imshow(np.mean(frames[trackstart:trackend], axis=0), origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "bad_line = plt.Line2D([], [], color=badcolor, label='Bad track')\n",
    "good_line = plt.Line2D([], [], color=goodcolor, label='Good track')\n",
    "plt.legend(handles=[bad_line, good_line], loc=\"upper left\")\n",
    "for i, isgood in enumerate(goodtrack):\n",
    "    ix = sts[:,i]\n",
    "    color = badcolor\n",
    "    if isgood:\n",
    "        color = goodcolor\n",
    "    ax.plot(x[ix,i], y[ix,i],marker=\"None\", color=color, markersize=14, linestyle=\"-\")\n",
    "    #ax.arrow(x[ix,i][-2], y[ix,i][-2],\n",
    "    #         np.diff(x[ix,i][-2:])[0], np.diff(y[ix,i][-2:])[0], \n",
    "    #         head_width=2, head_length=2, fc=color, ec=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(35, 168)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([898.])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[ix,i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Animate features\n",
    "#plt.rcParams['animation.ffmpeg_path'] = r\"E:\\install\\ffmpeg\\bin\\ffmpeg.exe\"\n",
    "\n",
    "# Prepare canvas\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "ax1 = plt.subplot(111,aspect=\"equal\")\n",
    "im1 = ax1.imshow(iframes[trackstart], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=\"limegreen\")\n",
    "ax1.plot(x[0,~goodtrack], y[0,~goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=\"red\")\n",
    "ax1.grid(color=\"white\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[trackstart].isoformat(), color=\"white\", fontsize=14)\n",
    "\n",
    "def animate(j):\n",
    "    im1.set_array(iframes[trackstart+j])\n",
    "    for line in plt.gca().get_lines():\n",
    "        if not line.get_linestyle()==\"None\":\n",
    "            line.remove()\n",
    "    for i, isgood in enumerate(goodtrack):\n",
    "        ix = np.where(sts[:j,i])[0]\n",
    "        color = \"red\"\n",
    "        if isgood:\n",
    "            color = \"limegreen\"\n",
    "        ax1.plot(x[ix,i], y[ix,i], marker=\"None\", color=color, markersize=14, linestyle=\"-\")\n",
    "    tstamp1.set_text(dtimes[trackstart+j].isoformat())\n",
    "    return im1\n",
    "\n",
    "# ATTENTION: THIS IS SLOW - Rendering each frame of the animation might take more time than the interval between the frames\n",
    "#    This can cause the temporal sequence to be confused in the matplotlib interactive mode.\n",
    "#    The animation thus looks better if saved as movie, or you have to increase the interval argument\n",
    "#    Animation not shown in notebook if you use %pylab inline\n",
    "ani = animation.FuncAnimation(fig, animate, frames=np.arange(trackstart, trackend-1), interval=400, blit=False)\n",
    "ani.save(\"featurescv.gif\", writer=\"imagemagick\", dpi=150)\n",
    "#ani.save(\"features.mp4\", writer=animation.FFMpegWriter(bitrate=2000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update tracked corners for each time step of the considered tracking period\n",
    "\n",
    "Until now, we only tracked those corners which we detected in the initial time step. We now want to add new tracks with each addtional time step, and follow these as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "List of # corners in each time step:\n",
      " [168, 161, 131, 109, 119, 177, 140, 200, 166, 190, 167, 198, 183, 200, 200, 200, 111, 200, 193, 200, 200, 200, 200, 200, 200, 200, 185, 200, 200, 165, 112, 162, 126, 149, 200, 200, 178, 199, 139, 123, 115, 153, 154, 200, 200, 136, 170, 200]\n"
     ]
    }
   ],
   "source": [
    "init_crns = [cv2.goodFeaturesToTrack(iframes[i], mask = None, **feature_params) for i in range(trackstart, trackend)]\n",
    "print(\"List of # corners in each time step:\\n\", [len(crn) for crn in init_crns ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this function wraps up everything which we already did above for a single set of corners\n",
    "def tracker(old, frameset, lk_params):\n",
    "    # Set containers to collect results (time steps in rows, corners in columns)\n",
    "    #   Tracking status\n",
    "    sts = np.zeros((trackend,len(old)), dtype=np.bool)\n",
    "    #   corner x coords\n",
    "    x = np.zeros((trackend,len(old))) * np.nan\n",
    "    #   corner y coords\n",
    "    y = np.zeros((trackend,len(old))) * np.nan\n",
    "    #   tracking error\n",
    "    errs = np.zeros((trackend,len(old))) * np.nan\n",
    "    #   Assign persistent corner IDs\n",
    "    ids = np.arange(len(old))\n",
    "    # Track good features\n",
    "    for i in range(len(frameset)-1):\n",
    "        # track current corners in next image\n",
    "        new, st, err = cv2.calcOpticalFlowPyrLK(prevImg=frameset[i], nextImg=frameset[i+1],\n",
    "                                                prevPts=old, nextPts=None, **lk_params)\n",
    "        success = st.ravel()==1\n",
    "        ids = ids[success]\n",
    "        sts[i, ids] = True\n",
    "        x[i, ids] = new[success,0,0]\n",
    "        y[i, ids] = new[success,0,1]\n",
    "        errs[i, ids] = err.ravel()[success]\n",
    "        # new corners will be old in the next loop\n",
    "        old = new[success]\n",
    "\n",
    "    # Incremental euclidic distance from starting point\n",
    "    trackdist = np.diff( np.sqrt( (x-x[0].reshape((1,-1)))**2 + (y-y[0].reshape((1,-1)))**2 ), axis=0 )\n",
    "    trackdist = np.vstack( (np.zeros((1,trackdist.shape[1])), trackdist))\n",
    "\n",
    "    # Find good tracks (but what is a \"good\" track...?)\n",
    "    goodtrack = np.zeros(x.shape[1], dtype=np.bool)\n",
    "    for i in range(len(goodtrack)):\n",
    "        # persistence of the track\n",
    "        if len(np.where(sts[:,i])[0]) < 2:\n",
    "            continue\n",
    "        # consistency of movement\n",
    "        if len(np.where(trackdist[:,i]<0)[0]) > 0:\n",
    "            continue\n",
    "        # tracking error\n",
    "        if len(np.where(errs[:,i]>15)[0]) > 5:\n",
    "            continue\n",
    "        goodtrack[i] = True\n",
    "    \n",
    "    return sts, x, y, errs, goodtrack    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\envs\\wradlib0_11_3\\lib\\site-packages\\ipykernel_launcher.py:39: RuntimeWarning: invalid value encountered in less\n",
      "C:\\Anaconda3\\envs\\wradlib0_11_3\\lib\\site-packages\\ipykernel_launcher.py:42: RuntimeWarning: invalid value encountered in greater\n"
     ]
    }
   ],
   "source": [
    "sts_ls, x_ls, y_ls, errs_ls, goodtrack_ls = [], [], [], [], []\n",
    "for i, crns in enumerate(init_crns):\n",
    "    sts, x, y, errs, goodtrack = tracker(crns, iframes[i:], lk_params)\n",
    "    sts_ls.append(sts)\n",
    "    x_ls.append(x)\n",
    "    y_ls.append(y)\n",
    "    errs_ls.append(errs)\n",
    "    goodtrack_ls.append(goodtrack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize tracks:\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "# average reflectivity as background image\n",
    "ax.imshow(np.mean(frames[trackstart:trackend], axis=0), origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect track properties (not in inline mode!)]\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "colors = [ plt.cm.spring(i) for i in np.linspace(0,254, len(goodtrack_ls)).astype(\"i4\") ]\n",
    "for j, goodtrack in enumerate(goodtrack_ls[:-2]):\n",
    "    sts, x, y = sts_ls[j], x_ls[j], y_ls[j]\n",
    "    for i, isgood in enumerate(goodtrack):\n",
    "        ix = sts[:,i]\n",
    "        # HERE WE DO NOT PLOT THE BAD TRACKS\n",
    "        color = \"none\"\n",
    "        if isgood:\n",
    "            color = colors[j]\n",
    "        ax.plot(x[ix,i], y[ix,i],marker=\"None\", color=color, linestyle=\"-\", alpha=0.4)\n",
    "        #ax.arrow(x[ix,i][-2], y[ix,i][-2],\n",
    "        #         np.diff(x[ix,i][-2:])[0], np.diff(y[ix,i][-2:])[0], \n",
    "        #         head_width=2, head_length=2, fc=color, ec=color, alpha=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ATTENTION: THIS ANIMATION TAKES A LONG WHILE (SEVERAL MINUTES) AND MIGHT STILL BE BUGGY\n",
    "\n",
    "# Prepare canvas\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "ax1 = plt.subplot(111,aspect=\"equal\")\n",
    "im1 = ax1.imshow(iframes[trackstart], origin=\"lower\", cmap=\"gray\", interpolation=\"none\", vmin=10,  vmax=60)\n",
    "plt.xlabel(\"Easting (km)\")\n",
    "plt.ylabel(\"Northing (km)\")\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "#ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=colors[0])\n",
    "ax1.grid(color=\"white\")\n",
    "plt.xlim(150,450)\n",
    "plt.ylim(550,900)\n",
    "tstamp1 = ax1.text(160, 560, dtimes[trackstart].isoformat(), color=\"white\", fontsize=14)\n",
    "colors = [ plt.cm.rainbow(i) for i in np.linspace(0,254, len(goodtrack_ls)).astype(\"i4\") ]\n",
    "def animate(j):\n",
    "    im1.set_array(iframes[trackstart+j])\n",
    "    tstamp1.set_text(dtimes[0+j].isoformat())\n",
    "    for line in plt.gca().get_lines():\n",
    "        line.remove()\n",
    "        #if not line.get_linestyle()==\"None\":\n",
    "        #    line.remove()\n",
    "    for k, goodtrack in enumerate(goodtrack_ls[:j]):\n",
    "        sts, x, y = sts_ls[k], x_ls[k], y_ls[k]\n",
    "        for i, isgood in enumerate(goodtrack):\n",
    "            ix = np.where(sts[:j,i])[0]\n",
    "            # HERE WE DO NOT PLOT THE BAD TRACKS\n",
    "            color = \"none\"\n",
    "            if isgood:\n",
    "                color = colors[k]\n",
    "            #ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle=\"None\", marker=\"o\", mfc=\"None\", mec=color, alpha=0.4)\n",
    "            ax1.plot(x[ix,i], y[ix,i],marker=\"None\", color=color, linestyle=\"-\", alpha=0.3)\n",
    "# ATTENTION: THIS IS SLOW - Rendering each frame of the animation might take more time than the interval between the frames\n",
    "#    This can cause the temporal sequence to be confused.\n",
    "#    The animation thus looks better if saved as movie, or you have to increase the interval argument\n",
    "ani = animation.FuncAnimation(fig, animate, frames=np.arange(trackstart, trackend), interval=400, blit=False)\n",
    "ani.save(\"featurescv.gif\", writer=\"imagemagick\", dpi=150)\n",
    "#ani.save(\"features2.avi\", dpi=500, bitrate=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment with SIFT/SURF feature detection and description\n",
    "\n",
    "See [SIFT and SURF](http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_table_of_contents_feature2d/py_table_of_contents_feature2d.html) for feature detection. Right now, this does not seem to add value as compared to the Optical Flow approach above. Features seem to be much less persistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SURF\n",
    "surf = cv2.xfeatures2d.SURF_create(3000)\n",
    "\n",
    "kplist = []\n",
    "deslist= []\n",
    "\n",
    "for i in range(trackstart, trackend):\n",
    "    kp, des = surf.detectAndCompute(iframes[i],None)\n",
    "    kplist.append(kp)\n",
    "    deslist.append(des)\n",
    "    print(\"Found %d keypoints in step %d.\" % (len(kp), i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "# average reflectivity as background image\n",
    "ax.imshow(frames[0], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "plt.xlabel(\"Easting (# pixels)\")\n",
    "plt.ylabel(\"Northing (# pixels)\")\n",
    "plt.title(\"[Zoom in to inspect feature properties (not in inline mode)]\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[0].isoformat(), color=\"white\", fontsize=14)\n",
    "plt.grid(color=\"white\")\n",
    "plt.xlim(0,nx/downsizeby)\n",
    "plt.ylim(0,nx/downsizeby)\n",
    "patches = []\n",
    "for kp_ in kplist[0]:\n",
    "    if kp_.size > 5:\n",
    "        circle = mpatches.Circle(kp_.pt, kp_.size, fill=False, edgecolor=\"red\")\n",
    "    #ax.add_patch(circle)\n",
    "    patches.append(circle)\n",
    "collection = PatchCollection(patches, facecolor=\"none\", edgecolor=\"red\")\n",
    "ax.add_collection(collection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make list of patch collections for all timesteps\n",
    "def collect(kp):\n",
    "    patches = []\n",
    "    for kp_ in kp:\n",
    "        if (kp_.size > 10) and (kp_.size < 50):\n",
    "            circle = mpatches.Circle(kp_.pt, kp_.size, fill=False, edgecolor=\"red\")\n",
    "            patches.append(circle)\n",
    "    return(PatchCollection(patches, facecolor=\"none\", edgecolor=\"red\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Animate features\n",
    "_plot_style = dict(markersize=12, markeredgewidth=2,\n",
    "                       markerfacecolor='none', markeredgecolor='r',\n",
    "                       marker='o', linestyle='none')\n",
    "_pcm_style = dict(cmap=plt.cm.spectral, vmin=0., vmax=30.)\n",
    "\n",
    "# Prepare canvas\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax1 = plt.subplot(111,aspect=\"equal\")\n",
    "im1 = ax1.imshow(iframes[0], origin=\"lower\", cmap=\"gray\", interpolation=\"none\")\n",
    "ax1.add_collection(collect(kplist[0]))\n",
    "ax1.grid(color=\"white\")\n",
    "tstamp1 = ax1.text(25, 850, dtimes[0].isoformat(), color=\"white\", fontsize=14)\n",
    "\n",
    "\n",
    "def animate(i):\n",
    "    im1.set_array(iframes[trackstart+i])\n",
    "    ax1.collections = []\n",
    "    ax1.add_collection(collect(kplist[trackstart+i]))\n",
    "    tstamp1.set_text(dtimes[trackstart+i].isoformat())\n",
    "    return im1\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate, frames=np.arange(trackstart, trackend-1), interval=200, blit=False)\n",
    "#ani.save(\"features_surf.avi\", dpi=400, bitrate=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Match features (brute force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According [Bowler et al. (2004)](http://www.sciencedirect.com/science/article/pii/S0022169403004591), maximum advection velocity of rainfall objects is about 130 km/h which is roughly 10 km (pixels) in 5 minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "maxveloc = 10.\n",
    "# Detect initial feature set \n",
    "detector = cv2.xfeatures2d.SURF_create(3000)\n",
    "kp1, des1 = detector.detectAndCompute(iframes[trackstart],None)\n",
    "\n",
    "# create BFMatcher object\n",
    "bf = cv2.BFMatcher()\n",
    "\n",
    "kp1_ls = []\n",
    "kp2_ls = []\n",
    "\n",
    "for i in range(trackstart+1, trackend):\n",
    "    kp2, des2 = detector.detectAndCompute(iframes[i],None)\n",
    "    matches = bf.knnMatch(des1, des2, k=1)\n",
    "    # Select matches to keep\n",
    "    kp1_, des1_, kp2_, des2_  = [], [], [], []\n",
    "    for match in matches:\n",
    "        match=match[0]\n",
    "        xy = np.vstack( (kp1[match.queryIdx].pt, kp2[match.trainIdx].pt) )\n",
    "        eucdist = np.sqrt( (xy[0,0] - xy[1,0])**2 + (xy[0,1] - xy[1,1])**2 )\n",
    "        if eucdist < maxveloc:\n",
    "            kp1_.append( kp1[match.queryIdx] )\n",
    "            des1_.append( np.array( des1[match.queryIdx] ) )\n",
    "            kp2_.append( kp2[match.trainIdx] )\n",
    "            des2_.append( np.array( des2[match.trainIdx] ) )\n",
    "    kp1_ls.append(kp1_)\n",
    "    kp2_ls.append(kp2_)\n",
    "    # Update initial feature set\n",
    "    kp1, des1 = kp2_, np.array( des2_ )"
   ]
  },
  {
   "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.6.4"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "119px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}