Dual Pol Algorithm Working Copy Custom Cell Tracking.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "## You are using the Python ARM Radar Toolkit (Py-ART), an open source\n",
      "## library for working with weather radar data. Py-ART is partly\n",
      "## supported by the U.S. Department of Energy as part of the Atmospheric\n",
      "## Radiation Measurement (ARM) Climate Research Facility, an Office of\n",
      "## Science user facility.\n",
      "##\n",
      "## If you use this software to prepare a publication, please cite:\n",
      "##\n",
      "##     JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mwilson41\\AppData\\Local\\Continuum\\Anaconda3\\envs\\radar\\lib\\site-packages\\urllib3\\contrib\\pyopenssl.py:46: DeprecationWarning: OpenSSL.rand is deprecated - you should use os.urandom instead\n",
      "  import OpenSSL.SSL\n"
     ]
    }
   ],
   "source": [
    "### This script will plot some graphics showing the relative locations of the ZDR and KDP maxes in the storm.\n",
    "import matplotlib.pyplot as plt\n",
    "import pyart\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "from metpy.units import atleast_1d, check_units, concatenate, units\n",
    "from matplotlib.patches import PathPatch\n",
    "from matplotlib.path import Path\n",
    "from siphon.radarserver import RadarServer\n",
    "rs = RadarServer('http://thredds-aws.unidata.ucar.edu/thredds/radarServer/nexrad/level2/S3/')\n",
    "#rs = RadarServer('http://thredds.ucar.edu/thredds/radarServer/nexrad/level2/IDD/')\n",
    "from datetime import datetime, timedelta\n",
    "from siphon.cdmr import Dataset\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from cartopy.io.shapereader import Reader\n",
    "from cartopy.feature import ShapelyFeature\n",
    "from metpy.units import atleast_1d, check_units, concatenate, units\n",
    "from shapely.geometry import polygon as sp\n",
    "import pyproj \n",
    "import shapely.ops as ops\n",
    "from shapely.ops import transform\n",
    "from shapely.geometry.polygon import Polygon\n",
    "from functools import partial\n",
    "from shapely import geometry\n",
    "import netCDF4\n",
    "from scipy import ndimage as ndi\n",
    "from skimage.feature import peak_local_max\n",
    "from skimage import data, img_as_float\n",
    "from pyproj import Geod\n",
    "from metpy.calc import get_wind_dir\n",
    "import matplotlib.lines as mlines\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Comment out or don't run this cell if you don't want to plot tornado tracks\n",
    "tors = np.genfromtxt('BoxingDayTors.csv', skip_header = 0, usecols =(5,10,15,16,17,18), delimiter = ',')\n",
    "time_tor = tors[:,0]\n",
    "start_torlats = np.asarray(tors[:,2])\n",
    "start_torlons = np.asarray(tors[:,3])\n",
    "end_torlats = np.asarray(tors[:,4])\n",
    "end_torlons = np.asarray(tors[:,5])\n",
    "rating = tors[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "in loop\n",
      "23\n",
      "0.454479\n",
      "2018-02-21 15:06:09.171000\n",
      "its this line\n",
      "heres the problem\n",
      "almost gridding\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "got gradient\n",
      "Set up our projection\n",
      "(-94.475833333333341, 29.655555555555555, <cartopy.crs.PlateCarree object at 0x000001F4006EFD00>)\n",
      "plotting\n",
      "1\n",
      "-93.65111555556936\n",
      "past mask\n",
      "47.7814\n",
      "found a storm\n",
      "found a big storm\n",
      "plotted a big storm\n",
      "made it to beginning of loop\n",
      "2\n",
      "-93.75234677705434\n",
      "added polygon\n",
      "2\n",
      "-93.58193847904353\n",
      "added polygon\n",
      "1\n",
      "-93.16596655610545\n",
      "past mask\n",
      "46.3471\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-93.26064986941574\n",
      "past mask\n",
      "46.5229\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-92.90038061562052\n",
      "past mask\n",
      "45.9702\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-93.35334698361922\n",
      "past mask\n",
      "48.4764\n",
      "found a storm\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "(0, -93.752346777054342)\n",
      "[]\n",
      "(1, -93.581938479043529)\n",
      "[]\n",
      "(2, -93.165966556105445)\n",
      "[]\n",
      "(3, -93.260649869415744)\n",
      "[]\n",
      "(4, -92.90038061562052)\n",
      "[]\n",
      "(5, -93.353346983619218)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "first dataframe\n",
      "oops kdp\n",
      "means there's a kdp problem\n",
      "KDP-ZDR separation didt work\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "0.455788\n",
      "2018-02-21 15:08:10.032000\n",
      "its this line\n",
      "heres the problem\n",
      "almost gridding\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "got gradient\n",
      "Set up our projection\n",
      "(-94.475833333333341, 29.655555555555555, <cartopy.crs.PlateCarree object at 0x000001F445897B48>)\n",
      "plotting\n",
      "1\n",
      "-93.6442911374783\n",
      "past mask\n",
      "47.9102\n",
      "found a storm\n",
      "found a big storm\n",
      "plotted a big storm\n",
      "made it to beginning of loop\n",
      "2\n",
      "-93.75285415450608\n",
      "[  1.00822357  42.87584139  56.16980309  47.81469709  84.00754685\n",
      "  68.58600253]\n",
      "Poly lon -93.75285415450608\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 0\n",
      "added polygon\n",
      "2\n",
      "-93.58655574121767\n",
      "[ 26.32280014  17.62477421  46.76474644  33.95331976  65.51636722\n",
      "  41.45373101]\n",
      "Poly lon -93.58655574121767\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 1\n",
      "added polygon\n",
      "2\n",
      "-93.5739953712764\n",
      "[ 44.03276325   2.10657702  58.38908595  44.64009476  67.60957104\n",
      "  26.0518242 ]\n",
      "Poly lon -93.5739953712764\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-93.16347570379371\n",
      "past mask\n",
      "46.5253\n",
      "found a storm\n",
      "[ 56.32142846  55.93252446   2.3944759   12.23819006  32.14902391\n",
      "  59.35251641]\n",
      "Poly lon -93.16347570379371\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-93.26657265077093\n",
      "past mask\n",
      "46.5866\n",
      "found a storm\n",
      "[ 47.52959814  41.44986434  16.1909302    2.73466675  36.21366719\n",
      "  46.61987248]\n",
      "Poly lon -93.26657265077093\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 3\n",
      "added polygon\n",
      "1\n",
      "-92.89580458577885\n",
      "past mask\n",
      "46.0749\n",
      "found a storm\n",
      "[ 84.69014136  67.56463586  35.72617235  37.68368903   2.31701571\n",
      "  55.4160453 ]\n",
      "Poly lon -92.89580458577885\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 4\n",
      "added polygon\n",
      "1\n",
      "-93.34771125725695\n",
      "past mask\n",
      "48.32\n",
      "found a storm\n",
      "[ 68.89977097  28.93242807  62.41913188  50.30383002  56.80202437\n",
      "   1.18678985]\n",
      "Poly lon -93.34771125725695\n",
      "[-93.75234678 -93.58193848 -93.16596656 -93.26064987 -92.90038062\n",
      " -93.35334698]\n",
      "[0 1 2 3 4 5]\n",
      "storm id 5\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "(0, -93.752854154506082)\n",
      "[]\n",
      "(1, -93.586555741217666)\n",
      "[]\n",
      "(2, -93.573995371276396)\n",
      "[]\n",
      "(3, -93.163475703793708)\n",
      "[]\n",
      "(4, -93.266572650770925)\n",
      "[]\n",
      "(5, -92.895804585778848)\n",
      "[]\n",
      "(6, -93.347711257256947)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "oops kdp\n",
      "means there's a kdp problem\n",
      "KDP-ZDR separation didt work\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "0.454098\n",
      "2018-02-21 15:09:20.891000\n",
      "its this line\n",
      "heres the problem\n",
      "almost gridding\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "got gradient\n",
      "Set up our projection\n",
      "(-94.475833333333341, 29.655555555555555, <cartopy.crs.PlateCarree object at 0x000001F400D959E8>)\n",
      "plotting\n",
      "1\n",
      "-93.64031009668133\n",
      "past mask\n",
      "47.9223\n",
      "found a storm\n",
      "found a big storm\n",
      "plotted a big storm\n",
      "made it to beginning of loop\n",
      "2\n",
      "-93.7467423546865\n",
      "[  1.19361934  25.97689131  43.79752803  55.78567601  46.99901466\n",
      "  84.16404058  68.57374907]\n",
      "Poly lon -93.7467423546865\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 0\n",
      "added polygon\n",
      "2\n",
      "-93.6765966791881\n",
      "nope\n",
      "2\n",
      "-93.57718057629195\n",
      "[ 28.09509514   1.02493144  19.08785545  45.2502257   31.7720316\n",
      "  65.03351628  41.7380906 ]\n",
      "Poly lon -93.57718057629195\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 6\n",
      "added polygon\n",
      "2\n",
      "-93.57063561060812\n",
      "[ 46.82445329  21.49760974   1.88610367  57.93904854  43.43638772\n",
      "  67.57935267  25.76259623]\n",
      "Poly lon -93.57063561060812\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-93.16201495586067\n",
      "past mask\n",
      "46.7342\n",
      "found a storm\n",
      "[ 56.55278901  45.09605136  55.59343728   1.82095207  13.33362211\n",
      "  32.66512861  58.50975771]\n",
      "Poly lon -93.16201495586067\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-93.26404802714288\n",
      "past mask\n",
      "46.5566\n",
      "found a storm\n",
      "[ 48.54499541  32.04445623  41.02854545  15.90445672   2.04294145\n",
      "  36.46057317  45.61272176]\n",
      "Poly lon -93.26404802714288\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 3\n",
      "added polygon\n",
      "1\n",
      "-92.89325381652354\n",
      "past mask\n",
      "46.2307\n",
      "found a storm\n",
      "[ 85.59474008  66.20176381  67.22528674  34.99265976  37.9375816\n",
      "   1.33015917  55.02732786]\n",
      "Poly lon -92.89325381652354\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 4\n",
      "added polygon\n",
      "1\n",
      "-93.34408306235841\n",
      "past mask\n",
      "48.3751\n",
      "found a storm\n",
      "[ 70.42389182  43.29520217  27.73196741  60.64158408  48.07295199\n",
      "  55.7498162    0.66030079]\n",
      "Poly lon -93.34408306235841\n",
      "[-93.75285415 -93.58655574 -93.57399537 -93.1634757  -93.26657265\n",
      " -92.89580459 -93.34771126]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 5\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "(0, -93.746742354686504)\n",
      "[]\n",
      "(1, -93.577180576291951)\n",
      "[]\n",
      "(2, -93.570635610608122)\n",
      "[]\n",
      "(3, -93.16201495586067)\n",
      "[]\n",
      "(4, -93.264048027142877)\n",
      "[]\n",
      "(5, -92.89325381652354)\n",
      "[]\n",
      "(6, -93.344083062358408)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "oops kdp\n",
      "means there's a kdp problem\n",
      "KDP-ZDR separation didt work\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "0.452831\n",
      "2018-02-21 15:11:12.915000\n",
      "its this line\n",
      "heres the problem\n",
      "almost gridding\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "got gradient\n",
      "Set up our projection\n",
      "(-94.475833333333341, 29.655555555555555, <cartopy.crs.PlateCarree object at 0x000001F445846150>)\n",
      "plotting\n",
      "1\n",
      "-93.63460322186826\n",
      "past mask\n",
      "48.0472\n",
      "found a storm\n",
      "found a big storm\n",
      "plotted a big storm\n",
      "made it to beginning of loop\n",
      "2\n",
      "-93.72397316394488\n",
      "[  3.28686572  23.63636812  42.56590813  53.74752434  45.00684156\n",
      "  81.98822106  65.97932598]\n",
      "Poly lon -93.72397316394488\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 0\n",
      "added polygon\n",
      "2\n",
      "-93.57490522376386\n",
      "[ 28.35832947   1.64862531  19.31421572  45.08856909  31.6503849\n",
      "  65.05020079  40.89733415]\n",
      "Poly lon -93.57490522376386\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 6\n",
      "added polygon\n",
      "2\n",
      "-93.52782843564465\n",
      "[ 37.66544032  10.9100085   11.82827464  46.08124456  31.66126509\n",
      "  61.09349804  31.56814033]\n",
      "Poly lon -93.52782843564465\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 6\n",
      "added polygon\n",
      "2\n",
      "-93.56235076775661\n",
      "[ 48.39521109  23.63539462   2.76391712  58.12425805  43.53472798\n",
      "  67.47295328  24.25633699]\n",
      "Poly lon -93.56235076775661\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-93.15953139761099\n",
      "past mask\n",
      "46.5409\n",
      "found a storm\n",
      "[ 56.2548471   43.79141647  55.31832376   2.17016775  13.24449122\n",
      "  32.14068167  56.95785734]\n",
      "Poly lon -93.15953139761099\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-93.25922369914873\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "past mask\n",
      "46.4864\n",
      "found a storm\n",
      "[ 48.84874443  31.11279168  40.65563972  16.30224657   2.58303125\n",
      "  35.93654776  43.68665636]\n",
      "Poly lon -93.25922369914873\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 3\n",
      "added polygon\n",
      "1\n",
      "-92.88892549193648\n",
      "past mask\n",
      "46.3435\n",
      "found a storm\n",
      "[ 85.5792149   65.75781457  67.37380769  35.04316172  37.94751236\n",
      "   1.6456802   54.44729496]\n",
      "Poly lon -92.88892549193648\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 4\n",
      "added polygon\n",
      "1\n",
      "-93.33785387016242\n",
      "past mask\n",
      "48.5164\n",
      "found a storm\n",
      "[ 70.79225368  43.9684441   27.71500622  60.18912429  47.47530752\n",
      "  55.58775761   1.59250577]\n",
      "Poly lon -93.33785387016242\n",
      "[-93.74674235 -93.57718058 -93.57063561 -93.16201496 -93.26404803\n",
      " -92.89325382 -93.34408306]\n",
      "[0 6 1 2 3 4 5]\n",
      "storm id 5\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "(0, -93.723973163944876)\n",
      "[]\n",
      "(1, -93.574905223763864)\n",
      "[]\n",
      "(2, -93.527828435644651)\n",
      "[]\n",
      "(3, -93.562350767756612)\n",
      "[]\n",
      "(4, -93.159531397610991)\n",
      "[]\n",
      "(5, -93.259223699148734)\n",
      "[]\n",
      "(6, -92.888925491936476)\n",
      "[]\n",
      "(7, -93.337853870162419)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "oops kdp\n",
      "means there's a kdp problem\n",
      "KDP-ZDR separation didt work\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n"
     ]
    }
   ],
   "source": [
    "storm_relative_dir = 180\n",
    "srdir = 180\n",
    "query = rs.query()\n",
    "#Here, set the initial time of the archived radar loop you want.\n",
    "dt = datetime(2017,7, 11, 23, 0) # Our specified time\n",
    "station = 'KMVX'\n",
    "query.stations(station).time_range(dt, dt + timedelta(hours=1.7))\n",
    "cat = rs.get_catalog(query)\n",
    "cat.datasets\n",
    "f = 27\n",
    "n = 1\n",
    "storm_index = 0\n",
    "scan_index = 0\n",
    "g = Geod(ellps='sphere')\n",
    "for item in sorted(cat.datasets.items()):\n",
    "    # After looping over the list of sorted datasets, pull the actual Dataset object out\n",
    "    # of our list of items and access over CDMRemote\n",
    "    try:\n",
    "        ds = item[1]\n",
    "        radar1 = pyart.io.nexrad_cdm.read_nexrad_cdm(ds.access_urls['OPENDAP'])\n",
    "        #Now let's calculate and plot specific differential phase.\n",
    "        #Pre-SAILS # for below: 17\n",
    "        #SAILS #: 25\n",
    "        if radar1.nsweeps > 25:\n",
    "            continue\n",
    "        for i in range(radar1.nsweeps):\n",
    "            print('in loop')\n",
    "            print(radar1.nsweeps)\n",
    "            radar = radar1.extract_sweeps([i])\n",
    "            if ((np.mean(radar.elevation['data']) < .65) and (np.max(np.asarray(radar.fields['differential_reflectivity']['data'])) != np.min(np.asarray(radar.fields['differential_reflectivity']['data'])))):\n",
    "                n = n+1\n",
    "                print(np.mean(radar.elevation['data']))\n",
    "                time_start = netCDF4.num2date(radar.time['data'][0], radar.time['units'])\n",
    "                print(time_start)\n",
    "                kdp_dict = pyart.retrieve.kdp_proc.kdp_maesaka(radar)\n",
    "                print('its this line')\n",
    "                radar.add_field('KDP', kdp_dict[0])\n",
    "                print('heres the problem')\n",
    "                # mask out last 10 gates of each ray, this removes the \"ring\" around the radar.\n",
    "                radar.fields['differential_reflectivity']['data'][:, -10:] = np.ma.masked\n",
    "                ref_ungridded = radar.fields['reflectivity']['data']\n",
    "                refl_c = np.copy(ref_ungridded)\n",
    "                ref_c = ma.masked_where(refl_c < 20., refl_c)\n",
    "                #Get ungridded ZDR\n",
    "                zdr_ungridded = radar.fields['differential_reflectivity']['data']\n",
    "                zdrl_c = np.copy(zdr_ungridded)\n",
    "                zdr_c = ma.masked_where(refl_c < 20, zdrl_c)\n",
    "\n",
    "                ungrid_lons = radar.gate_longitude['data']\n",
    "                ungrid_lats = radar.gate_latitude['data']\n",
    "                # exclude masked gates from the gridding\n",
    "                gatefilter = pyart.filters.GateFilter(radar)\n",
    "                gatefilter.exclude_masked('differential_reflectivity')\n",
    "                print('almost gridding')\n",
    "                #Now let's grid the data\n",
    "                grid = pyart.map.grid_from_radars(\n",
    "                    (radar,), gatefilters=(gatefilter, ),\n",
    "                    grid_shape=(1, 500, 500),\n",
    "                    grid_limits=((200, 200), (-123000.0, 123000.0), (-123000.0, 123000.0)),\n",
    "                    fields=['differential_reflectivity','reflectivity','KDP','cross_correlation_ratio'])\n",
    "                #Get gate altitudes\n",
    "                gate_altitude = radar.gate_altitude['data'][:]\n",
    "                #Get the data from the grid\n",
    "                ZDR = grid.fields['differential_reflectivity']['data'][0]\n",
    "                REF = grid.fields['reflectivity']['data'][0]\n",
    "                KDP = grid.fields['KDP']['data'][0]\n",
    "                CC = grid.fields['cross_correlation_ratio']['data'][0]\n",
    "\n",
    "                #Mask everything below 20dbz\n",
    "                #import numpy.ma as ma\n",
    "                ZDRmasked1 = ma.masked_where(REF < 20, ZDR)\n",
    "                REFmasked = ma.masked_where(REF < 20, REF)\n",
    "                #Use a 45 dBZ mask for KDP to only get areas in the storm core\n",
    "                KDPmasked = ma.masked_where(REF < 50, KDP)\n",
    "                KDPmasked = ma.filled(KDPmasked, fill_value = -2)\n",
    "                #Filter out spots not in forward flank using gradient direction\n",
    "                print('made it to smoothing')\n",
    "                smoothed_ref1 = ndi.gaussian_filter(REFmasked, sigma = 2, order = 0)\n",
    "                REFgradient = np.asarray(np.gradient(smoothed_ref1))\n",
    "                REFgradient[0,:,:] = ma.masked_where(REF < 20, REFgradient[0,:,:])\n",
    "                REFgradient[1,:,:] = ma.masked_where(REF < 20, REFgradient[1,:,:])\n",
    "                print('made it through gradient')\n",
    "                grad_dir1 = get_wind_dir(REFgradient[1,:,:] * units('m/s'), REFgradient[0,:,:] * units('m/s'))\n",
    "                grad_dir = ma.masked_where(REF < 20, grad_dir1)\n",
    "                #srdir = storm_relative_dir\n",
    "                grad_ffd = np.abs(np.arctan2(np.sin(grad_dir * units('degrees')-srdir * units('degrees')), np.cos(grad_dir * units('degrees')-srdir * units('degrees'))))\n",
    "                grad_ffd = grad_ffd.to('degrees')\n",
    "                print('got gradient')\n",
    "                #Add a fill value for the ZDR mask so that contours will be closed\n",
    "                ZDRmasked2 = ma.masked_where(grad_ffd > 120 * units('degrees'), ZDRmasked1)\n",
    "                ZDRmasked = ma.masked_where(CC < .60, ZDRmasked2)\n",
    "                ZDRmasked = ma.filled(ZDRmasked, fill_value = -2)\n",
    "\n",
    "                rlons = grid.point_longitude['data']\n",
    "                rlats = grid.point_latitude['data']\n",
    "                rlons_2d = rlons[0,:,:]\n",
    "                rlats_2d = rlats[0,:,:]\n",
    "                cenlat = radar.latitude['data'][0]\n",
    "                cenlon = radar.longitude['data'][0]\n",
    "                #Let's set up the map projection!\n",
    "                print('Set up our projection')\n",
    "                crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)\n",
    "\n",
    "                # Set up our array of latitude and longitude values and transform to \n",
    "                # the desired projection.\n",
    "\n",
    "                tlatlons = crs.transform_points(ccrs.LambertConformal(central_longitude=265, central_latitude=25, standard_parallels=(25.,25.)),rlons[0,:,:],rlats[0,:,:])\n",
    "                tlons = tlatlons[:,:,0]\n",
    "                tlats = tlatlons[:,:,1]\n",
    "\n",
    "                # Limit the extent of the map area, must convert to proper coords.\n",
    "                LL = (cenlon-1.5,cenlat-1.5,ccrs.PlateCarree())\n",
    "                UR = (cenlon+1.5,cenlat+1.5,ccrs.PlateCarree())\n",
    "                print(LL)\n",
    "\n",
    "                # Get data to plot state and province boundaries\n",
    "                states_provinces = cfeature.NaturalEarthFeature(\n",
    "                        category='cultural',\n",
    "                        name='admin_1_states_provinces_lakes',\n",
    "                        scale='50m',\n",
    "                        facecolor='none')\n",
    "                fname = 'cb_2016_us_county_20m/cb_2016_us_county_20m.shp'\n",
    "                fname2 = 'cb_2016_us_state_20m/cb_2016_us_state_20m.shp'\n",
    "                counties = ShapelyFeature(Reader(fname).geometries(),ccrs.PlateCarree(), facecolor = 'none', edgecolor = 'black')\n",
    "                states = ShapelyFeature(Reader(fname2).geometries(),ccrs.PlateCarree(), facecolor = 'none', edgecolor = 'black')\n",
    "                fig=plt.figure(n,figsize=(30.,25.))\n",
    "                ax = plt.subplot(111,projection=ccrs.PlateCarree())\n",
    "                ax.coastlines('50m',edgecolor='black',linewidth=0.75)\n",
    "                #ax.add_feature(states_provinces,edgecolor='black',linewidth=0.5)\n",
    "                ax.add_feature(counties, edgecolor = 'black', linewidth = 0.5)\n",
    "                ax.add_feature(states, edgecolor = 'black', linewidth = 1.5)\n",
    "                ax.set_extent([LL[0],UR[0],LL[1],UR[1]])\n",
    "                REFlevels = np.arange(20,73,2)\n",
    "                print('plotting')\n",
    "                #ax.contourf(rlons[0,:,:],rlats[0,:,:],REFmasked,REFlevels,cmap = plt.cm.gist_ncar)\n",
    "                refp = ax.pcolormesh(ungrid_lons, ungrid_lats, ref_c, cmap=plt.cm.gist_ncar, vmin = 10, vmax = 73)\n",
    "                #zdrp = ax.pcolormesh(ungrid_lons, ungrid_lats, zdr_c, cmap=plt.cm.nipy_spectral, vmin = -2, vmax = 6)\n",
    "                #Plot local maxes in reflectivity\n",
    "                # image_max is the dilation of im with a 20*20 structuring element\n",
    "                # It is used within peak_local_max function\n",
    "                smoothed_ref = ndi.gaussian_filter(REFmasked, sigma = 3, order = 0)\n",
    "                #image_max = ndi.maximum_filter(smoothed_ref, size=15, mode='constant')\n",
    "\n",
    "                # Comparison between image_max and im to find the coordinates of local maxima\n",
    "                #coordinates = peak_local_max(smoothed_ref, min_distance=20)\n",
    "                #ref_maxes = REFmasked[coordinates[:,0], coordinates[:,1]]\n",
    "                #max_lons = rlons[0,coordinates[:,0], coordinates[:,1]]\n",
    "                #max_lats = rlats[0,coordinates[:,0], coordinates[:,1]]\n",
    "                #max_lons_c = max_lons[ref_maxes > 45]\n",
    "                #max_lats_c = max_lats[ref_maxes > 45]\n",
    "                \n",
    "                REFlev = [45]\n",
    "                REFlev1 = [50]\n",
    "                refc = ax.contour(rlons[0,:,:],rlats[0,:,:],smoothed_ref,REFlev, alpha=.01)\n",
    "\n",
    "                ref_areas = []\n",
    "                max_lons_c = []\n",
    "                max_lats_c = []\n",
    "                storm_ids = []\n",
    "                \n",
    "                proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "                           pyproj.Proj(init='epsg:3857'))\n",
    "\n",
    "                #Look for reflectivity centroids\n",
    "                for col in refc.collections:\n",
    "                    # Loop through all polygons that have the same intensity level\n",
    "                    for contour_path in col.get_paths(): \n",
    "                        # Create the polygon for this intensity level\n",
    "                        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "                            print(1)\n",
    "                            cpa = np.asarray(cp[:])\n",
    "                            x = cpa[:,0]\n",
    "                            y = cpa[:,1]\n",
    "                            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "                            if ncp == 0:\n",
    "                                poly = new_shape\n",
    "                            else:\n",
    "                                # Remove the holes if there are any\n",
    "                                poly = poly.difference(new_shape)\n",
    "\n",
    "                        # do something with polygon\n",
    "                        #print(poly.area) \n",
    "                        print(poly.centroid.x)\n",
    "                        s_new = transform(proj, poly)\n",
    "                        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "                        boundary = np.asarray(poly.boundary.xy)\n",
    "                        polypath = Path(boundary.transpose())\n",
    "                        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "                        maskr = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "                        meanr = np.mean(smoothed_ref[maskr])\n",
    "                        print('past mask')\n",
    "                        print(meanr)\n",
    "                        if projected_area > 10 * units('km^2') and meanr > REFlev[0]:\n",
    "                            print('found a storm')\n",
    "                            #For big blobs with embedded supercells, find the embedded storm cores\n",
    "                            #Normal 'big storm' cutoff 300 km^2\n",
    "                            if projected_area > 300 * units('km^2'):\n",
    "                                print('found a big storm')\n",
    "                                rlon_2 = rlons[0,:,:]\n",
    "                                rlat_2 = rlats[0,:,:]\n",
    "                                #smoothed_ref_m = ma.MaskedArray(smoothed_ref, mask=maskr)\n",
    "                                smoothed_ref_m = ma.masked_where(maskr==False, smoothed_ref)\n",
    "                                smoothed_ref_m = ma.filled(smoothed_ref_m, fill_value = -2)\n",
    "                                rlon2m = ma.MaskedArray(rlon_2, mask=maskr)\n",
    "                                rlat2m = ma.MaskedArray(rlat_2, mask=maskr)\n",
    "                                refc1 = ax.contour(rlon2m,rlat2m,smoothed_ref_m,REFlev1, linewidths = 3, linestyle = '--', alpha=.01)\n",
    "                                #refc1 = ax.contour(rlon_2[maskr],rlat_2[maskr],smoothed_ref[maskr],REFlev1, colors = 'g', linewidths = 3)\n",
    "                                print('plotted a big storm')\n",
    "                                #Look for reflectivity centroids\n",
    "                                for col1 in refc1.collections:\n",
    "                                    # Loop through all polygons that have the same intensity level\n",
    "                                    print('made it to beginning of loop')\n",
    "                                    for contour_path1 in col1.get_paths(): \n",
    "                                        # Create the polygon for this intensity level\n",
    "                                        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                                        for ncp1,cp1 in enumerate(contour_path1.to_polygons()):\n",
    "                                            print(2)\n",
    "                                            cpa1 = np.asarray(cp1[:])\n",
    "                                            x1 = cpa1[:,0]\n",
    "                                            y1 = cpa1[:,1]\n",
    "                                            new_shape1 = geometry.Polygon([(i[0], i[1]) for i in zip(x1,y1)])\n",
    "                                            if ncp1 == 0:\n",
    "                                                poly1 = new_shape1\n",
    "                                            else:\n",
    "                                                # Remove the holes if there are any\n",
    "                                                poly1 = poly1.difference(new_shape)\n",
    "\n",
    "                                        # do something with polygon\n",
    "                                        #print(poly.area) \n",
    "                                        print(poly1.centroid.x)\n",
    "                                        s_new1 = transform(proj, poly1)\n",
    "                                        projected_area1 = (transform(proj, poly1).area * units('m^2')).to('km^2')\n",
    "                                        if projected_area1 > 10 * units('km^2'):\n",
    "                                            ref_areas.append((projected_area1))\n",
    "                                            max_lons_c.append((poly1.centroid.x))\n",
    "                                            max_lats_c.append((poly1.centroid.y))\n",
    "                                            if scan_index == 0:\n",
    "                                                storm_ids.append((storm_index))\n",
    "                                                storm_index = storm_index + 1\n",
    "                                            else:\n",
    "                                                #dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                                max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "                                                max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "                                                storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "                                                dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                                for i in range(max_lons_p.shape[0]):\n",
    "                                                    distance_track = g.inv(poly1.centroid.x, poly1.centroid.y,\n",
    "                                                                           max_lons_p[i], max_lats_p[i])\n",
    "                                                    dist_track[i] = distance_track[2]/1000.\n",
    "                                                print(dist_track)\n",
    "                                                print('Poly lon', poly1.centroid.x)\n",
    "                                                print(max_lons_p)\n",
    "                                                print(storm_ids_p)\n",
    "                                                if np.min(dist_track) < 10.0:\n",
    "                                                    storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "                                                    print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                                else:\n",
    "                                                    storm_ids.append((storm_index))\n",
    "                                                    print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                                    storm_index = storm_index + 1\n",
    "                                            print('added polygon')\n",
    "                                        else:\n",
    "                                            print('nope')\n",
    "                            else:\n",
    "                                ref_areas.append((projected_area))\n",
    "                                max_lons_c.append((poly.centroid.x))\n",
    "                                max_lats_c.append((poly.centroid.y))\n",
    "                                if scan_index == 0:\n",
    "                                    storm_ids.append((storm_index))\n",
    "                                    storm_index = storm_index + 1\n",
    "                                else:\n",
    "                                    #dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                    max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "                                    max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "                                    storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "                                    dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                    for i in range(max_lons_p.shape[0]):\n",
    "                                        distance_track = g.inv(poly.centroid.x, poly.centroid.y,\n",
    "                                                               max_lons_p[i], max_lats_p[i])\n",
    "                                        dist_track[i] = distance_track[2]/1000.\n",
    "                                    print(dist_track)\n",
    "                                    print('Poly lon', poly.centroid.x)\n",
    "                                    print(max_lons_p)\n",
    "                                    print(storm_ids_p)\n",
    "                                    if np.min(dist_track) < 10.0:\n",
    "                                        storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "                                        print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                    else:\n",
    "                                        storm_ids.append((storm_index))\n",
    "                                        print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                        storm_index = storm_index + 1\n",
    "                                print('added polygon')\n",
    "                                            \n",
    "                    #print(s_new)\n",
    "                max_lons_c = np.asarray(max_lons_c)\n",
    "                max_lats_c = np.asarray(max_lats_c)\n",
    "                #ZDRlevels = np.arange(3, np.max(ZDRmasked)+((np.max(ZDRmasked))-3)/2, (np.max(ZDRmasked))-3)\n",
    "                zdrlev = [3.5]\n",
    "                kdplev = [1.5]\n",
    "                #ZDRlevels = np.arange(3,5.5,.5)\n",
    "                #ZDRlevels1 = np.arange(5,10,.5)\n",
    "                #KDPlevels = np.arange(.75, np.max(KDPmasked)+((np.max(KDPmasked))-1.5)/2, (np.max(KDPmasked))-1.5)\n",
    "                #KDPlevels = np.arange(.75,1.75,.25)\n",
    "                #KDPlevels1 = np.arange(1.5,10,.25)\n",
    "                zdrc = ax.contour(rlons[0,:,:],rlats[0,:,:],ZDRmasked,zdrlev,linewidths = 2, colors='k', alpha = .01)\n",
    "                #zrdc = ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked,ZDRlevels,linewide = .01, colors='pink', alpha = .8)\n",
    "                #ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked,ZDRlevels1,linewide = .01, colors='crimson', alpha = .8)\n",
    "\n",
    "                #kdpc = ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels,linewide = .01, colors ='green', alpha = .5)\n",
    "                kdpc = ax.contour(rlons[0,:,:],rlats[0,:,:],KDPmasked,kdplev,linewidths = 2, colors='green', alpha = .8)\n",
    "                #ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels1,linewide = .01, colors ='b', alpha = .5)\n",
    "                print('made it here')\n",
    "                plt.savefig('testfig.png')\n",
    "\n",
    "                #proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "                #           pyproj.Proj(init='epsg:3857'))\n",
    "                zdr_areas = []\n",
    "                zdr_centroid_lon = []\n",
    "                zdr_centroid_lat = []\n",
    "                zdr_mean = []\n",
    "                zdr_cc_mean = []\n",
    "                zdr_max = []\n",
    "                zdr_storm_lon = []\n",
    "                zdr_storm_lat = []\n",
    "                zdr_dist = []\n",
    "                zdr_forw = []\n",
    "                zdr_back = []\n",
    "                #print(\"here too\")\n",
    "                for col in zdrc.collections:\n",
    "                    # Loop through all polygons that have the same intensity level\n",
    "                    #print('hi')\n",
    "                    for contour_path in col.get_paths(): \n",
    "                        # Create the polygon for this intensity level\n",
    "                        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                        #print('hi')\n",
    "                        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "                            #print('hi')\n",
    "                            cpa = np.asarray(cp[:])\n",
    "                            x = cpa[:,0]\n",
    "                            y = cpa[:,1]\n",
    "                            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "                            if ncp == 0:\n",
    "                                poly = new_shape\n",
    "                                #print('hi')\n",
    "                            else:\n",
    "                                # Remove the holes if there are any\n",
    "                                poly = poly.difference(new_shape)\n",
    "                                #print('hi')\n",
    "\n",
    "                        # do something with polygon\n",
    "                        #print(poly.area) \n",
    "                        #print(poly.centroid.x)\n",
    "                        s_new = transform(proj, poly)\n",
    "                        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "                        boundary = np.asarray(poly.boundary.xy)\n",
    "                        polypath = Path(boundary.transpose())\n",
    "                        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "                        mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "                        mean = np.mean(ZDRmasked[mask])\n",
    "                        mean_cc = np.mean(CC[mask])\n",
    "                        if projected_area > 1 * units('km^2') and mean > zdrlev[0] and mean_cc > .88:\n",
    "                            g = Geod(ellps='sphere')\n",
    "                            dist = np.zeros((np.asarray(max_lons_c).shape[0]))\n",
    "                            forw = np.zeros((np.asarray(max_lons_c).shape[0]))\n",
    "                            back = np.zeros((np.asarray(max_lons_c).shape[0]))\n",
    "                            for i in range(dist.shape[0]):\n",
    "                                        distance_1 = g.inv(poly.centroid.x, poly.centroid.y,\n",
    "                                                               max_lons_c[i], max_lats_c[i])\n",
    "                                        #print(distance_1[2]/1000)\n",
    "                                        #print(distance_1)\n",
    "                                        back[i] = distance_1[1]\n",
    "                                        if distance_1[1] < 0:\n",
    "                                            back[i] = distance_1[1] + 360\n",
    "                                        forw[i] = np.abs(back[i] - storm_relative_dir)\n",
    "                                        dist[i] = distance_1[2]/1000.\n",
    "                            #print(dist.shape)\n",
    "                            if (forw[np.where(dist == np.min(dist))[0][0]] < 90 and np.min(dist) < 20.0) or (forw[np.where(dist == np.min(dist))[0][0]] < 140 and np.min(dist) < 6.0):\n",
    "                                zdr_storm_lon.append((max_lons_c[np.where(dist == np.min(dist))[0][0]]))\n",
    "                                zdr_storm_lat.append((max_lats_c[np.where(dist == np.min(dist))[0][0]]))\n",
    "                                zdr_dist.append(np.min(dist))\n",
    "                                zdr_forw.append(forw[np.where(dist == np.min(dist))[0][0]])\n",
    "                                zdr_back.append(back[np.where(dist == np.min(dist))[0][0]])\n",
    "                                zdr_areas.append((projected_area))\n",
    "                                zdr_centroid_lon.append((poly.centroid.x))\n",
    "                                zdr_centroid_lat.append((poly.centroid.y))\n",
    "                                zdr_mean.append((mean))\n",
    "                                zdr_cc_mean.append((mean_cc))\n",
    "                                zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "                                patch = PathPatch(polypath, facecolor='none', alpha=.8, edgecolor = 'blue', linewidth = 3)\n",
    "                                ax.add_patch(patch)\n",
    "                        #print(s_new)\n",
    "                print('made it through zdr centroids')\n",
    "                \n",
    "                if len(zdr_storm_lon) > 0:\n",
    "                    kdp_areas = []\n",
    "                    kdp_centroid_lon = []\n",
    "                    kdp_centroid_lat = []\n",
    "                    kdp_max = []\n",
    "                    kdp_storm_lon = []\n",
    "                    kdp_storm_lat = []\n",
    "                    for col in kdpc.collections:\n",
    "                        # Loop through all polygons that have the same intensity level\n",
    "                        for contour_path in col.get_paths(): \n",
    "                            # Create the polygon for this intensity level\n",
    "                            # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                            for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "                                print(1)\n",
    "                                cpa = np.asarray(cp[:])\n",
    "                                x = cpa[:,0]\n",
    "                                y = cpa[:,1]\n",
    "                                new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "                                if ncp == 0:\n",
    "                                    poly = new_shape\n",
    "                                else:\n",
    "                                    # Remove the holes if there are any\n",
    "                                    poly = poly.difference(new_shape)\n",
    "\n",
    "                            # do something with polygon\n",
    "                            #print(poly.area) \n",
    "                            print(poly.centroid.x)\n",
    "                            s_new = transform(proj, poly)\n",
    "                            projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "                            boundary = np.asarray(poly.boundary.xy)\n",
    "                            polypath = Path(boundary.transpose())\n",
    "                            coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "                            mask_kdp = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "                            #mean = np.mean(ZDRmasked[mask])\n",
    "                            #mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "                            #mean = np.mean(REFmasked[mask])\n",
    "                            if projected_area > 2 * units('km^2'):\n",
    "                                g = Geod(ellps='sphere')\n",
    "                                dist_kdp = np.zeros((np.asarray(max_lons_c).shape[0]))\n",
    "                                for i in range(dist_kdp.shape[0]):\n",
    "                                            distance_kdp = g.inv(poly.centroid.x, poly.centroid.y,\n",
    "                                                                   max_lons_c[i], max_lats_c[i])\n",
    "                                            #print(distance_1[2]/1000)\n",
    "                                            #print(\"KDP dist:\", distance_kdp)\n",
    "                                            dist_kdp[i] = distance_kdp[2]/1000.\n",
    "                                print(dist_kdp)\n",
    "                                if np.min(np.asarray(dist_kdp)) < 15.0:\n",
    "                                    #print('Got to KDP stuff')\n",
    "                                    kdp_areas.append((projected_area))\n",
    "                                    kdp_centroid_lon.append((poly.centroid.x))\n",
    "                                    kdp_centroid_lat.append((poly.centroid.y))\n",
    "                                    kdp_storm_lon.append((max_lons_c[np.where(dist_kdp == np.min(dist_kdp))[0][0]]))\n",
    "                                    kdp_storm_lat.append((max_lats_c[np.where(dist_kdp == np.min(dist_kdp))[0][0]]))\n",
    "                                    kdp_max.append((np.max(KDPmasked[mask_kdp])))\n",
    "                                    patch = PathPatch(polypath, facecolor='none', alpha=.5, edgecolor = 'grey', linewidth = 3)\n",
    "                                    ax.add_patch(patch)\n",
    "\n",
    "                    print('made it through kdp centroids')\n",
    "\n",
    "                    #Consolidating the arc objects associated with each storm:\n",
    "                    zdr_areas_arr = np.zeros((len(zdr_areas)))\n",
    "                    zdr_max_arr = np.zeros((len(zdr_max)))\n",
    "                    zdr_mean_arr = np.zeros((len(zdr_mean)))                    \n",
    "                    for i in range(len(zdr_areas)):\n",
    "                        zdr_areas_arr[i] = zdr_areas[i].magnitude\n",
    "                        zdr_max_arr[i] = zdr_max[i]\n",
    "                        zdr_mean_arr[i] = zdr_mean[i]\n",
    "\n",
    "                    zdr_centroid_lons = np.asarray(zdr_centroid_lon)\n",
    "                    zdr_centroid_lats = np.asarray(zdr_centroid_lat)\n",
    "                    zdr_con_areas = []\n",
    "                    zdr_con_maxes = []\n",
    "                    zdr_con_means = []\n",
    "                    zdr_con_centroid_lon = []\n",
    "                    zdr_con_centroid_lat = []\n",
    "                    zdr_con_max_lon = []\n",
    "                    zdr_con_max_lat = []\n",
    "                    zdr_con_storm_lon = []\n",
    "                    zdr_con_storm_lat = []\n",
    "\n",
    "                    #For KDP as well\n",
    "                    kdp_areas_arr = np.zeros((len(kdp_areas)))\n",
    "                    kdp_max_arr = np.zeros((len(kdp_max)))\n",
    "                    for i in range(len(kdp_areas)):\n",
    "                        kdp_areas_arr[i] = kdp_areas[i].magnitude\n",
    "                        kdp_max_arr[i] = kdp_max[i]\n",
    "                    kdp_centroid_lons = np.asarray(kdp_centroid_lon)\n",
    "                    kdp_centroid_lats = np.asarray(kdp_centroid_lat)\n",
    "                    kdp_con_areas = []\n",
    "                    kdp_con_maxes = []\n",
    "                    kdp_con_centroid_lon = []\n",
    "                    kdp_con_centroid_lat = []\n",
    "                    kdp_con_max_lon = []\n",
    "                    kdp_con_max_lat = []\n",
    "                    kdp_con_storm_lon = []\n",
    "                    kdp_con_storm_lat = []\n",
    "                    for i in enumerate(zdr_storm_lon):\n",
    "                        print(i[0])\n",
    "                        if i[0] != 0:\n",
    "                            if zdr_storm_lon[i[0]-1] == zdr_storm_lon[i[0]]:\n",
    "                                #print(\"Skipping this one\")\n",
    "                                continue\n",
    "                            else:\n",
    "                                print(zdr_storm_lon[i[0]])\n",
    "                                #Find the arc objects associated with this storm:\n",
    "                                zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                #print(\"zdr lons:\", zdr_objects_lons)\n",
    "                                #Get the sum of their areas\n",
    "                                print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                zdr_con_maxes.append(np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                zdr_con_means.append(np.mean(zdr_mean_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                zdr_con_max_lon.append(rlons_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                zdr_con_max_lat.append(rlats_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                #print(\"Areas sum:\", zdr_con_areas)\n",
    "                                #Find the actual centroids\n",
    "                                weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                zdr_con_centroid_lon.append(np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                zdr_con_centroid_lat.append(np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                zdr_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                zdr_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                                #print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "                                #print(\"lons in loop\", zdr_objects_lons)\n",
    "\n",
    "                                try:\n",
    "                                    #Find the kdp objects associated with this storm:\n",
    "                                    kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    #print(\"kdp lons:\", kdp_objects_lons)\n",
    "                                    kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    #Get the sum of their areas\n",
    "                                    print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                    kdp_con_areas.append(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    kdp_con_maxes.append(np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    kdp_con_max_lon.append(rlons_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                    kdp_con_max_lat.append(rlats_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                    #Find the actual centroids\n",
    "                                    weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    kdp_con_centroid_lon.append(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    #print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                    kdp_con_centroid_lat.append(np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    kdp_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                    kdp_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                                except:\n",
    "                                    print('storm missing kdp or zdr')\n",
    "                                    kdp_con_areas.append(0)\n",
    "                                    kdp_con_maxes.append(0)\n",
    "                                    kdp_con_max_lon.append(0)\n",
    "                                    kdp_con_max_lat.append(0)\n",
    "                                    kdp_con_centroid_lon.append(0)\n",
    "                                    kdp_con_centroid_lat.append(0)\n",
    "                                    kdp_con_storm_lon.append(0)\n",
    "                                    kdp_con_storm_lat.append(0)\n",
    "\n",
    "\n",
    "\n",
    "                        else:\n",
    "                            #print(zdr_storm_lon[i[0]])\n",
    "                            #Find the arc objects associated with this storm:\n",
    "                            zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                            zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                            #print(\"zdr lons:\", zdr_objects_lons)\n",
    "                            #print(\"arc lats:\", zdr_objects_lats)\n",
    "                            #Get the sum of their areas\n",
    "                            #print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                            zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                            zdr_con_maxes.append(np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                            zdr_con_means.append(np.mean(zdr_mean_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                            zdr_con_max_lon.append(rlons_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                            zdr_con_max_lat.append(rlats_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                            #print(\"Areas sum:\",zdr_con_areas)\n",
    "                            #Find the actual centroids\n",
    "                            weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                            zdr_con_centroid_lon.append(np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                            weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                            zdr_con_centroid_lat.append(np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                            zdr_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                            zdr_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                            #print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "                            #print(\"lons out of loop\", zdr_objects_lons)\n",
    "                            try:\n",
    "                                #Find the kdp objects associated with this storm:\n",
    "                                kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                #print(\"kdp lons:\", kdp_objects_lons)\n",
    "                                kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                #Get the sum of their areas\n",
    "                                #print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                kdp_con_areas.append(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                kdp_con_maxes.append(np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                kdp_con_max_lon.append(rlons_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                kdp_con_max_lat.append(rlats_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                #Find the actual centroids\n",
    "                                weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                kdp_con_centroid_lon.append(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                #print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                kdp_con_centroid_lat.append(np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                kdp_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                kdp_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                            except:\n",
    "                                print('storm missing kdp or zdr')\n",
    "                                kdp_con_areas.append(0)\n",
    "                                kdp_con_maxes.append(0)\n",
    "                                kdp_con_max_lon.append(0)\n",
    "                                kdp_con_max_lat.append(0)\n",
    "                                kdp_con_centroid_lon.append(0)\n",
    "                                kdp_con_centroid_lat.append(0)\n",
    "                                kdp_con_storm_lon.append(0)\n",
    "                                kdp_con_storm_lat.append(0)\n",
    "                                    \n",
    "                    #Calculate KDP-ZDR separation\n",
    "                    print('calculating separation')\n",
    "                    kdp_con_centroid_lons1 = np.asarray(kdp_con_centroid_lon)\n",
    "                    kdp_con_centroid_lats1 = np.asarray(kdp_con_centroid_lat)\n",
    "                    zdr_con_centroid_lons1 = np.asarray(zdr_con_centroid_lon)\n",
    "                    zdr_con_centroid_lats1 = np.asarray(zdr_con_centroid_lat)\n",
    "                    #Eliminate consolidated arcs smaller than a specified area\n",
    "                    area = 2\n",
    "                    zdr_con_areas_arr = np.asarray(zdr_con_areas)\n",
    "                    zdr_con_centroid_lats = zdr_con_centroid_lats1[zdr_con_areas_arr > area]\n",
    "                    zdr_con_centroid_lons = zdr_con_centroid_lons1[zdr_con_areas_arr > area]\n",
    "                    kdp_con_centroid_lats = kdp_con_centroid_lats1[zdr_con_areas_arr > area]\n",
    "                    kdp_con_centroid_lons = kdp_con_centroid_lons1[zdr_con_areas_arr > area]\n",
    "                    zdr_con_max_lons1 = np.asarray(zdr_con_max_lon)[zdr_con_areas_arr > area]\n",
    "                    zdr_con_max_lats1 = np.asarray(zdr_con_max_lat)[zdr_con_areas_arr > area]\n",
    "                    kdp_con_max_lons1 = np.asarray(kdp_con_max_lon)[zdr_con_areas_arr > area]\n",
    "                    kdp_con_max_lats1 = np.asarray(kdp_con_max_lat)[zdr_con_areas_arr > area]\n",
    "                    print('Boolean problem here')\n",
    "                    zdr_con_areas1 = zdr_con_areas_arr[zdr_con_areas_arr > area]\n",
    "\n",
    "                    kdp_inds = np.where(kdp_con_centroid_lats > 0)\n",
    "                    distance_kdp_zdr = g.inv(kdp_con_centroid_lons[kdp_inds], kdp_con_centroid_lats[kdp_inds], zdr_con_centroid_lons[kdp_inds], zdr_con_centroid_lats[kdp_inds])\n",
    "                    dist_kdp_zdr = distance_kdp_zdr[2] / 1000.\n",
    "                    #Now make an array for the distances which will have the same shape as the lats to prevent errors\n",
    "                    shaped_dist = np.zeros((np.shape(zdr_con_areas)))\n",
    "                    shaped_dist[kdp_inds] = dist_kdp_zdr\n",
    "                    print('maybe its here')\n",
    "                    #Do the same for the distances between the maxes\n",
    "                    distance_kdp_zdr_max = g.inv(kdp_con_max_lons1[kdp_inds], kdp_con_max_lats1[kdp_inds], zdr_con_max_lons1[kdp_inds], zdr_con_max_lats1[kdp_inds])\n",
    "                    dist_kdp_zdr_max = distance_kdp_zdr_max[2] / 1000.\n",
    "                    #Now make an array for the distances which will have the same shape as the lats to prevent errors\n",
    "                    shaped_dist_max = np.zeros((np.shape(zdr_con_areas)))\n",
    "                    shaped_dist_max[kdp_inds] = dist_kdp_zdr_max\n",
    "                    print('or not')\n",
    "                    \n",
    "                else:\n",
    "                    kdp_areas = []\n",
    "                    kdp_centroid_lon = []\n",
    "                    kdp_centroid_lat = []\n",
    "                    kdp_storm_lon = []\n",
    "                    kdp_storm_lat = []\n",
    "                    zdr_con_centroid_lats = []\n",
    "                    zdr_con_centroid_lons = []\n",
    "                    kdp_con_centroid_lats = []\n",
    "                    kdp_con_centroid_lons = []\n",
    "                    kdp_con_area = []\n",
    "                    zdr_con_areas1 = []\n",
    "                    \n",
    "                ###Now let's consolidate everything to fit the Pandas dataframe!\n",
    "                p_zdr_areas = []\n",
    "                p_zdr_maxes = []\n",
    "                p_zdr_means = []\n",
    "                p_separations = []\n",
    "                for storm in enumerate(max_lons_c):\n",
    "                    print(storm)\n",
    "                    print(np.flatnonzero(np.isclose(max_lons_c[storm[0]], zdr_con_storm_lon, rtol=1e-05)))\n",
    "                    matching_ind = np.flatnonzero(np.isclose(max_lons_c[storm[0]], zdr_con_storm_lon, rtol=1e-05))\n",
    "                    if matching_ind.shape[0] > 0:\n",
    "                        p_zdr_areas.append((zdr_con_areas[matching_ind[0]]))\n",
    "                        p_zdr_maxes.append((zdr_con_maxes[matching_ind[0]]))\n",
    "                        p_zdr_means.append((zdr_con_means[matching_ind[0]]))\n",
    "                        p_separations.append((shaped_dist[matching_ind[0]]))\n",
    "                    else:\n",
    "                        p_zdr_areas.append((0))\n",
    "                        p_zdr_maxes.append((0))\n",
    "                        p_zdr_means.append((0))\n",
    "                        p_separations.append((0))\n",
    "                    \n",
    "                #Now start plotting stuff!\n",
    "                print('made it through giant if statement')\n",
    "                if np.asarray(zdr_centroid_lon).shape[0] > 0:\n",
    "                    ax.scatter(zdr_centroid_lon, zdr_centroid_lat, marker = '*', s = 100, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                    #ax.scatter(zdr_con_max_lon, zdr_con_max_lat, marker = '*', s = 100, color = 'purple', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                if np.asarray(kdp_centroid_lon).shape[0] > 0:\n",
    "                    ax.scatter(kdp_centroid_lon, kdp_centroid_lat, marker = '^', s = 100, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                    #ax.scatter(kdp_con_max_lon, kdp_con_max_lat, marker = '^', s = 100, color = 'purple', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                print(\"plotted centroids\")\n",
    "                #Uncomment to print all object areas\n",
    "                #for i in enumerate(zdr_areas):\n",
    "                #    plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2\" %(zdr_areas[i[0]].magnitude), size = 23)\n",
    "                    #plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2 / %.2f km / %.2f dB\" %(zdr_areas[i[0]].magnitude, zdr_dist[i[0]], zdr_forw[i[0]]), size = 23)\n",
    "                    #plt.annotate(zdr_areas[i[0]], (zdr_centroid_lon[i[0]],zdr_centroid_lat[i[0]]))\n",
    "                #ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels1,linewide = .01, colors ='b', alpha = .5)\n",
    "                #plt.tight_layout()\n",
    "                #plt.savefig('ZDRarcannotated.png')\n",
    "                storm_times = []\n",
    "                for l in range(len(max_lons_c)):\n",
    "                    storm_times.append((time_start))\n",
    "                #Now record all data in a Pandas dataframe, using similar code as Tint (kind of borrowed from it)\n",
    "                print('making dataframe')\n",
    "                new_cells = pd.DataFrame({\n",
    "                    'scan': scan_index,\n",
    "                    'storm_id' : storm_ids,\n",
    "                    'storm_id1' : storm_ids,\n",
    "                    'storm_lon' : max_lons_c,\n",
    "                    'storm_lat' : max_lats_c,\n",
    "                    'zdr_area' : p_zdr_areas,\n",
    "                    'zdr_max' : p_zdr_maxes,\n",
    "                    'zdr_mean' : p_zdr_means,\n",
    "                    'kdp_zdr_sep' : p_separations,\n",
    "                    'times' : storm_times\n",
    "                })\n",
    "                print('setting index')\n",
    "                new_cells.set_index(['scan', 'storm_id'], inplace=True)\n",
    "                if scan_index == 0:\n",
    "                    print('first dataframe')\n",
    "                    tracks_dataframe = new_cells\n",
    "                else:\n",
    "                    tracks_dataframe = tracks_dataframe.append(new_cells)\n",
    "                n = n+1\n",
    "                scan_index = scan_index + 1\n",
    "                #max_lons_p = max_lons_c\n",
    "                #max_lats_p = max_lats_c\n",
    "                #storm_ids_p = storm_ids\n",
    "                #Plot the consolidated stuff!\n",
    "                if len(zdr_con_areas1) > 0:\n",
    "                    #ax.scatter(zdr_con_centroid_lon, zdr_con_centroid_lat, marker = '*', s = 500, color = 'orange', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                    try:\n",
    "                        for i in enumerate(zdr_con_centroid_lats):\n",
    "                            print(\"consolidated ZDR:\")\n",
    "                            ax.scatter(zdr_con_centroid_lons, zdr_con_centroid_lats, marker = '*', s = 500, color = 'orange', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                            try:\n",
    "                                plt.text(zdr_con_centroid_lons[i[0]]+.025, zdr_con_centroid_lats[i[0]]+.016, \"%.2f km^2 / %.2f dB\" %(zdr_con_areas1[i[0]], zdr_con_maxes[i[0]]), size = 23)\n",
    "                            except:\n",
    "                                print(\"oops zdr\")\n",
    "                        #plt.text(kdp_con_centroid_lon[i[0]]-.20, kdp_con_centroid_lat[i[0]]+.016, \"%.2f km\" %(dist_kdp_zdr[i[0]]), size = 23, color = 'red')                \n",
    "                    except:\n",
    "                        print('failed')\n",
    "                        try:\n",
    "                            plt.text(float(zdr_con_centroid_lons)+.016, float(zdr_con_centroid_lats)+.016, \"%.2f km^2 / %.2f dB\" %(zdr_con_areas1[i[0]], zdr_con_maxes[i[0]]), size = 23)\n",
    "                        except:\n",
    "                            print('no zdr centroids')\n",
    "                        #plt.text(float(kdp_con_centroid_lon)-.20, float(kdp_con_centroid_lat)+.016, \"%.2f km\" %(float(dist_kdp_zdr[0])), size = 23, color = 'red')\n",
    "                if len(kdp_con_areas) > 0:\n",
    "                    #ax.scatter(zdr_con_centroid_lon, zdr_con_centroid_lat, marker = '*', s = 500, color = 'orange', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                    try:\n",
    "                        for i in kdp_inds[0]:\n",
    "                            #plt.text(zdr_con_centroid_lon[i[0]]+.025, zdr_con_centroid_lat[i[0]]+.016, \"%.2f km^2\" %(zdr_con_areas[i[0]].magnitude), size = 23)\n",
    "                            try:\n",
    "                                plt.text(kdp_con_centroid_lons[i]-.20, kdp_con_centroid_lats[i]+.016, \"%.2f km\" %(shaped_dist[i]), size = 23, color = 'red')                \n",
    "                            except:\n",
    "                                print('oops kdp')\n",
    "                    except:\n",
    "                        print('failed')\n",
    "                        #plt.text(float(zdr_con_centroid_lon)+.016, float(zdr_con_centroid_lat)+.016, \"%.2f km^2\" %(float(zdr_con_areas[0])), size = 23)\n",
    "                        try:\n",
    "                            plt.text(float(kdp_con_centroid_lons)-.20, float(kdp_con_centroid_lats)+.016, \"%.2f km\" %(float(shaped_dist[0])), size = 23, color = 'red')\n",
    "                        except:\n",
    "                            print('no kdp centroids')\n",
    "                print(\"means there's a kdp problem\")\n",
    "                hour = time_start.hour\n",
    "                if hour < 10:\n",
    "                    hour = '0'+str(hour)\n",
    "                minute = time_start.minute\n",
    "                if minute < 10:\n",
    "                    minute = '0'+str(minute)\n",
    "                day = time_start.day\n",
    "                if day < 10:\n",
    "                    day = '0'+str(day)\n",
    "                title_plot = plt.title(station+' Radar Reflectivity, ZDR, and KDP '+str(time_start.year)+'-'+str(time_start.month)+'-'+str(time_start.day)+\n",
    "                                           ' '+str(hour)+':'+str(minute)+' UTC', size = 25)\n",
    "                #if np.asarray(zdr_storm_lon).shape[0] > 0:\n",
    "                #    ax.scatter(zdr_storm_lon,zdr_storm_lat, marker = \"o\", color = 'purple', s = 500)\n",
    "                #if np.asarray(kdp_storm_lon).shape[0] > 0:\n",
    "                #    ax.scatter(kdp_storm_lon,kdp_storm_lat, marker = \"o\", color = 'purple', s = 500)\n",
    "                try:\n",
    "                    ax.scatter(max_lons_c,max_lats_c, marker = \"o\", color = 'k', s = 500, alpha = .6)\n",
    "                except:\n",
    "                    \"No storm centroids found\"\n",
    "                try:\n",
    "                    plt.plot([zdr_con_centroid_lons[kdp_inds], kdp_con_centroid_lons[kdp_inds]], [zdr_con_centroid_lats[kdp_inds],kdp_con_centroid_lats[kdp_inds]], color = 'k', linewidth = 5, transform=ccrs.PlateCarree())\n",
    "                except:\n",
    "                    print('KDP-ZDR separation didt work')\n",
    "                ref_centroid_lon = max_lons_c\n",
    "                ref_centroid_lat = max_lats_c\n",
    "                for i in enumerate(ref_centroid_lon): \n",
    "                    plt.text(ref_centroid_lon[i[0]]+.016, ref_centroid_lat[i[0]]+.016, \"storm_id: %.1f\" %(storm_ids[i[0]]), size = 25)\n",
    "                #Comment out this line if not plotting tornado tracks\n",
    "                #plt.plot([start_torlons, end_torlons], [start_torlats, end_torlats], color = 'purple', linewidth = 5, transform=ccrs.PlateCarree())\n",
    "                #Add legend stuff\n",
    "                zdr_outline = mlines.Line2D([], [], color='blue', linewidth = 5, linestyle = 'solid', label='ZDR Arc Outline(Area/Max)')\n",
    "                kdp_outline = mlines.Line2D([], [], color='green', linewidth = 5,linestyle = 'solid', label='\"KDP Foot\" Outline')\n",
    "                separation_vector = mlines.Line2D([], [], color='black', linewidth = 5,linestyle = 'solid', label='KDP/ZDR Centroid Separation Vector (Red Text=Distance)')\n",
    "                tor_track = mlines.Line2D([], [], color='purple', linewidth = 5,linestyle = 'solid', label='Tornado Tracks')\n",
    "                elevation = mlines.Line2D([], [], color='grey', linewidth = 5,linestyle = 'solid', label='Height AGL (m)')\n",
    "\n",
    "                plt.legend(handles=[zdr_outline, kdp_outline, separation_vector, tor_track, elevation], loc = 3, fontsize = 25)\n",
    "                alt_levs = [1000, 2000]\n",
    "                cele = ax.contour(ungrid_lons,ungrid_lats,gate_altitude,alt_levs, linewidths = 7, alpha = .6, colors = 'grey')\n",
    "                plt.clabel(cele, fontsize=18, inline=1, inline_spacing=10, fmt='%i', rightside_up=True, use_clabeltext=True)\n",
    "                plt.tight_layout()\n",
    "                plt.savefig('Slide_ZDRArc_comp'+station+str(time_start.year)+str(time_start.month)+str(day)+str(hour)+str(minute)+'.png')\n",
    "                print('figure saved')\n",
    "                plt.close()\n",
    "\n",
    "    except:\n",
    "        continue\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Save off the dataframe as a pickle file\n",
    "tracks_dataframe.to_pickle('ExampleCase_7_11_17KMVX.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tracks_dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(max_lons_p)\n",
    "print(max_lats_p)\n",
    "print(storm_ids_p)\n",
    "print(max_lons_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "for i in range(max_lons_p.shape[0]):\n",
    "    distance_track = g.inv(poly.centroid.x, poly.centroid.y,\n",
    "                           max_lons_p[i], max_lats_p[i])\n",
    "    dist_track[i] = distance_track[2]/1000.\n",
    "print(dist_track)\n",
    "print(poly.centroid.x)\n",
    "print(max_lons_p)\n",
    "print(storm_ids_p)\n",
    "if np.min(dist_track) < 10.0:\n",
    "    storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "else:\n",
    "    storm_ids.append((storm_index))\n",
    "    storm_index = storm_index + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "radar = radar1.extract_sweeps([16])\n",
    "print(np.mean(radar.elevation['data']))\n",
    "print(netCDF4.num2date(radar.time['data'][0], radar.time['units']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "REFlev = [45]\n",
    "REFlev1 = [50]\n",
    "refc = ax.contour(rlons[0,:,:],rlats[0,:,:],smoothed_ref,REFlev, alpha=.01)\n",
    "\n",
    "ref_areas = []\n",
    "max_lons_c = []\n",
    "max_lats_c = []\n",
    "storm_ids = []\n",
    "\n",
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "           pyproj.Proj(init='epsg:3857'))\n",
    "\n",
    "#Look for reflectivity centroids\n",
    "for col in refc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print(1)\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        boundary = np.asarray(poly.boundary.xy)\n",
    "        polypath = Path(boundary.transpose())\n",
    "        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "        maskr = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "        meanr = np.mean(smoothed_ref[maskr])\n",
    "        print('past mask')\n",
    "        print(meanr)\n",
    "        if projected_area > 10 * units('km^2') and meanr > REFlev[0]:\n",
    "            print('found a storm')\n",
    "            #For big blobs with embedded supercells, find the embedded storm cores\n",
    "            if projected_area > 300 * units('km^2'):\n",
    "                print('found a big storm')\n",
    "                rlon_2 = rlons[0,:,:]\n",
    "                rlat_2 = rlats[0,:,:]\n",
    "                #smoothed_ref_m = ma.MaskedArray(smoothed_ref, mask=maskr)\n",
    "                smoothed_ref_m = ma.masked_where(maskr==False, smoothed_ref)\n",
    "                smoothed_ref_m = ma.filled(smoothed_ref_m, fill_value = -2)\n",
    "                rlon2m = ma.MaskedArray(rlon_2, mask=maskr)\n",
    "                rlat2m = ma.MaskedArray(rlat_2, mask=maskr)\n",
    "                refc1 = ax.contour(rlon2m,rlat2m,smoothed_ref_m,REFlev1, linewidths = 3, linestyle = '--', alpha=.01)\n",
    "                #refc1 = ax.contour(rlon_2[maskr],rlat_2[maskr],smoothed_ref[maskr],REFlev1, colors = 'g', linewidths = 3)\n",
    "                print('plotted a big storm')\n",
    "                #Look for reflectivity centroids\n",
    "                for col1 in refc1.collections:\n",
    "                    # Loop through all polygons that have the same intensity level\n",
    "                    print('made it to beginning of loop')\n",
    "                    for contour_path1 in col1.get_paths(): \n",
    "                        # Create the polygon for this intensity level\n",
    "                        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                        for ncp1,cp1 in enumerate(contour_path1.to_polygons()):\n",
    "                            print(2)\n",
    "                            cpa1 = np.asarray(cp1[:])\n",
    "                            x1 = cpa1[:,0]\n",
    "                            y1 = cpa1[:,1]\n",
    "                            new_shape1 = geometry.Polygon([(i[0], i[1]) for i in zip(x1,y1)])\n",
    "                            if ncp1 == 0:\n",
    "                                poly1 = new_shape1\n",
    "                            else:\n",
    "                                # Remove the holes if there are any\n",
    "                                poly1 = poly1.difference(new_shape)\n",
    "\n",
    "                        # do something with polygon\n",
    "                        #print(poly.area) \n",
    "                        print(poly1.centroid.x)\n",
    "                        s_new1 = transform(proj, poly1)\n",
    "                        projected_area1 = (transform(proj, poly1).area * units('m^2')).to('km^2')\n",
    "                        if projected_area1 > 10 * units('km^2'):\n",
    "                            ref_areas.append((projected_area1))\n",
    "                            max_lons_c.append((poly1.centroid.x))\n",
    "                            max_lats_c.append((poly1.centroid.y))\n",
    "                            if scan_index == 0:\n",
    "                                storm_ids.append((storm_index))\n",
    "                                storm_index = storm_index + 1\n",
    "                            else:\n",
    "                                #dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "                                max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "                                storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "                                dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                for i in range(max_lons_p.shape[0]):\n",
    "                                    distance_track = g.inv(poly1.centroid.x, poly1.centroid.y,\n",
    "                                                           max_lons_p[i], max_lats_p[i])\n",
    "                                    dist_track[i] = distance_track[2]/1000.\n",
    "                                print(dist_track)\n",
    "                                print('Poly lon', poly.centroid.x)\n",
    "                                print(max_lons_p)\n",
    "                                print(storm_ids_p)\n",
    "                                if np.min(dist_track) < 10.0:\n",
    "                                    storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "                                    print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                else:\n",
    "                                    storm_ids.append((storm_index))\n",
    "                                    print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                    storm_index = storm_index + 1\n",
    "                            print('added polygon')\n",
    "                        else:\n",
    "                            print('nope')\n",
    "            else:\n",
    "                ref_areas.append((projected_area))\n",
    "                max_lons_c.append((poly.centroid.x))\n",
    "                max_lats_c.append((poly.centroid.y))\n",
    "                if scan_index == 0:\n",
    "                    storm_ids.append((storm_index))\n",
    "                    storm_index = storm_index + 1\n",
    "                else:\n",
    "                    #dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                    max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "                    max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "                    storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "                    dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                    for i in range(max_lons_p.shape[0]):\n",
    "                        distance_track = g.inv(poly.centroid.x, poly.centroid.y,\n",
    "                                               max_lons_p[i], max_lats_p[i])\n",
    "                        dist_track[i] = distance_track[2]/1000.\n",
    "                    print(dist_track)\n",
    "                    print('Poly lon', poly.centroid.x)\n",
    "                    print(max_lons_p)\n",
    "                    print(storm_ids_p)\n",
    "                    if np.min(dist_track) < 10.0:\n",
    "                        storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "                        print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                    else:\n",
    "                        storm_ids.append((storm_index))\n",
    "                        print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                        storm_index = storm_index + 1\n",
    "                print('added polygon')\n",
    "\n",
    "    #print(s_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "storm_times = []\n",
    "for l in range(20):\n",
    "    storm_times.append((time_start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(storm_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(tracks_dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:]))\n",
    "max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "print(max_lons_p)\n",
    "max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "storm_ids_p = np.asarray(tracks_dataframe['storm_id'].loc[scan_index-1].iloc[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(storm_ids)\n",
    "for ids in enumerate(storm_ids):\n",
    "    print(ids[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(zdr_con_centroid_lons1)\n",
    "print(kdp_con_centroid_lons1)\n",
    "print(shaped_dist_max)\n",
    "#print(zdr_max_lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kdp_con_centroid_lons1 = np.asarray(kdp_con_centroid_lon)\n",
    "kdp_con_centroid_lats1 = np.asarray(kdp_con_centroid_lat)\n",
    "zdr_con_centroid_lons1 = np.asarray(zdr_con_centroid_lon)\n",
    "zdr_con_centroid_lats1 = np.asarray(zdr_con_centroid_lat)\n",
    "#Eliminate consolidated arcs smaller than a specified area\n",
    "area = 2\n",
    "zdr_con_areas_arr = np.asarray(zdr_con_areas)\n",
    "zdr_con_centroid_lats = zdr_con_centroid_lats1[zdr_con_areas_arr > area]\n",
    "zdr_con_centroid_lons = zdr_con_centroid_lons1[zdr_con_areas_arr > area]\n",
    "kdp_con_centroid_lats = kdp_con_centroid_lats1[zdr_con_areas_arr > area]\n",
    "kdp_con_centroid_lons = kdp_con_centroid_lons1[zdr_con_areas_arr > area]\n",
    "zdr_con_max_lons1 = np.asarray(zdr_con_max_lon)[zdr_con_areas_arr > area]\n",
    "zdr_con_max_lats1 = np.asarray(zdr_con_max_lat)[zdr_con_areas_arr > area]\n",
    "kdp_con_max_lons1 = np.asarray(kdp_con_max_lon)[zdr_con_areas_arr > area]\n",
    "kdp_con_max_lats1 = np.asarray(kdp_con_max_lat)[zdr_con_areas_arr > area]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zdr_con_max_lons1 = np.asarray(zdr_con_max_lon)[zdr_con_areas_arr > area]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(np.asarray(zdr_con_max_lon)[zdr_con_areas_arr > area])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(zdr_mean)\n",
    "print(zdr_cc_mean)\n",
    "print(zdr_con_maxes)\n",
    "rlons_2 = rlons[0,:,:]\n",
    "rlats_2 = rlats[0,:,:]\n",
    "print(np.where(ZDRmasked==zdr_con_maxes))\n",
    "print(rlons_2[np.where(ZDRmasked==zdr_con_maxes)])\n",
    "print(rlats_2[np.where(ZDRmasked==zdr_con_maxes)])\n",
    "zdr_con_max_lon = rlons_2[np.where(ZDRmasked==zdr_con_maxes)]\n",
    "zdr_con_max_lat = rlats_2[np.where(ZDRmasked==zdr_con_maxes)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(zdr_con_centroid_lats.shape)\n",
    "print(kdp_con_centroid_lats.shape)\n",
    "print(shaped_dist.shape)\n",
    "print(np.asarray(zdr_con_areas).shape)\n",
    "for i in enumerate(zdr_con_areas):\n",
    "        print(\"consolidated ZDR:\")\n",
    "        print(zdr_con_centroid_lon[i[0]]+.025, zdr_con_centroid_lat[i[0]]+.016, zdr_con_areas[i[0]])\n",
    "print(zdr_con_areas)\n",
    "\n",
    "zdr_con_areas_arr = np.asarray(zdr_con_areas)\n",
    "zdr_con_centroid_lats1 = zdr_con_centroid_lats[zdr_con_areas_arr > 10]\n",
    "zdr_con_centroid_lons1 = zdr_con_centroid_lons[zdr_con_areas_arr > 10]\n",
    "kdp_con_centroid_lats1 = kdp_con_centroid_lats[zdr_con_areas_arr > 10]\n",
    "kdp_con_centroid_lons1 = kdp_con_centroid_lons[zdr_con_areas_arr > 10]\n",
    "shaped_dist1 = shaped_dist[zdr_con_areas_arr > 10]\n",
    "zdr_con_areas1 = zdr_con_areas_arr[zdr_con_areas_arr > 10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#print(max_lons_c)\n",
    "#print(max_lats_c)\n",
    "#print(radar.fields)\n",
    "print(gate_altitude)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(shaped_dist)\n",
    "print(kdp_con_centroid_lats)\n",
    "print(kdp_inds)\n",
    "for i in kdp_inds[0]:\n",
    "    print(shaped_dist[i])\n",
    "    print(kdp_con_centroid_lats[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "REFlev = [40]\n",
    "refc = ax.contour(rlons[0,:,:],rlats[0,:,:],smoothed_ref,REFlev)\n",
    "\n",
    "ref_areas = []\n",
    "max_lons_c = []\n",
    "max_lats_c = []\n",
    "\n",
    "#Look for reflectivity centroids\n",
    "for col in refc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print(1)\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        if projected_area > 20 * units('km^2'):\n",
    "            ref_areas.append((projected_area))\n",
    "            max_lons_c.append((poly.centroid.x))\n",
    "            max_lats_c.append((poly.centroid.y))\n",
    "            print('added polygon')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(len(zdr_con_areas))\n",
    "print(zdr_con_areas)\n",
    "#print(dist_kdp_zdr)\n",
    "print(zdr_con_centroid_lat)\n",
    "print(kdp_inds)\n",
    "print(kdp_con_centroid_lats[kdp_inds])\n",
    "print(kdp_con_centroid_lats)\n",
    "print(kdp_con_centroid_lat)\n",
    "print(kdp_con_centroid_lats[kdp_inds])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig=plt.figure(n,figsize=(30.,25.))\n",
    "ax = plt.subplot(111,projection=ccrs.PlateCarree())\n",
    "ax.coastlines('50m',edgecolor='black',linewidth=0.75)\n",
    "#ax.add_feature(states_provinces,edgecolor='black',linewidth=0.5)\n",
    "ax.add_feature(counties, edgecolor = 'black', linewidth = 0.5)\n",
    "ax.add_feature(states, edgecolor = 'black', linewidth = 1.5)\n",
    "ax.set_extent([LL[0],UR[0],LL[1],UR[1]])\n",
    "# It is used within peak_local_max function\n",
    "smoothed_ref = ndi.gaussian_filter(REFmasked, sigma = 3, order = 0)\n",
    "image_max = ndi.maximum_filter(smoothed_ref, size=15, mode='constant')\n",
    "\n",
    "# Comparison between image_max and im to find the coordinates of local maxima\n",
    "coordinates = peak_local_max(smoothed_ref, min_distance=20)\n",
    "ref_maxes = REFmasked[coordinates[:,0], coordinates[:,1]]\n",
    "max_lons = rlons[0,coordinates[:,0], coordinates[:,1]]\n",
    "max_lats = rlats[0,coordinates[:,0], coordinates[:,1]]\n",
    "#max_lons_c = max_lons[ref_maxes > 45]\n",
    "#max_lats_c = max_lats[ref_maxes > 45]\n",
    "alt_levs = [2000]\n",
    "ax.contour(ungrid_lons,ungrid_lats,gate_altitude,alt_levs, linewidths = 7)\n",
    "REFlevels = np.arange(20,73,2)\n",
    "REFlev = [40]\n",
    "ax.contourf(rlons[0,:,:],rlats[0,:,:],smoothed_ref,REFlevels,cmap = plt.cm.gist_ncar)\n",
    "REFlev = [40]\n",
    "refc = ax.contour(rlons[0,:,:],rlats[0,:,:],smoothed_ref,REFlev)\n",
    "\n",
    "ref_areas = []\n",
    "max_lons_c = []\n",
    "max_lats_c = []\n",
    "\n",
    "#Look for reflectivity centroids\n",
    "for col in refc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print(1)\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        if projected_area > 20 * units('km^2'):\n",
    "            ref_areas.append((projected_area))\n",
    "            max_lons_c.append((poly.centroid.x))\n",
    "            max_lats_c.append((poly.centroid.y))\n",
    "\n",
    "    #print(s_new)\n",
    "\n",
    "#Uncomment to print all object areas\n",
    "#for i in enumerate(zdr_areas):\n",
    "#    plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2\" %(zdr_areas[i[0]].magnitude), size = 23)\n",
    "    #plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2 / %.2f km / %.2f dB\" %(zdr_areas[i[0]].magnitude, zdr_dist[i[0]], zdr_forw[i[0]]), size = 23)\n",
    "    #plt.annotate(zdr_areas[i[0]], (zdr_centroid_lon[i[0]],zdr_centroid_lat[i[0]]))\n",
    "#ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels1,linewide = .01, colors ='b', alpha = .5)\n",
    "#plt.tight_layout()\n",
    "#plt.savefig('ZDRarcannotated.png')\n",
    "\n",
    "\n",
    "hour = time_start.hour\n",
    "if hour < 10:\n",
    "    hour = '0'+str(hour)\n",
    "minute = time_start.minute\n",
    "if minute < 10:\n",
    "    minute = '0'+str(minute)\n",
    "day = time_start.day\n",
    "if day < 10:\n",
    "    day = '0'+str(day)\n",
    "title_plot = plt.title(station+' Radar Reflectivity, ZDR, and KDP '+str(time_start.year)+'-'+str(time_start.month)+'-'+str(time_start.day)+\n",
    "                           ' '+str(hour)+':'+str(minute)+' UTC', size = 25)\n",
    "#if np.asarray(zdr_storm_lon).shape[0] > 0:\n",
    "#    ax.scatter(zdr_storm_lon,zdr_storm_lat, marker = \"o\", color = 'purple', s = 500)\n",
    "#if np.asarray(kdp_storm_lon).shape[0] > 0:\n",
    "#    ax.scatter(kdp_storm_lon,kdp_storm_lat, marker = \"o\", color = 'purple', s = 500)\n",
    "ax.scatter(max_lons_c,max_lats_c, marker = \"o\", color = 'k', s = 500, alpha = .6)\n",
    "try:\n",
    "    plt.plot([zdr_con_centroid_lons[kdp_inds], kdp_con_centroid_lons[kdp_inds]], [zdr_con_centroid_lats[kdp_inds],kdp_con_centroid_lats[kdp_inds]], color = 'red', linewidth = 5, transform=ccrs.PlateCarree())\n",
    "except:\n",
    "    print('KDP-ZDR separation didt work')\n",
    "#plt.savefig('FFD_ZDRArc_comp'+station+str(time_start.year)+str(time_start.month)+str(day)+str(hour)+str(minute)+'.png')\n",
    "print('figure saved')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(len(kdp_con_areas))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#THIS CELL HAS THE PROBLEM MAKE DIST_KDP THE SAME SHAPE AS EVERYTHING ELSE\n",
    "\n",
    "kdp_inds = np.where(kdp_con_centroid_lats > 0)\n",
    "print(kdp_inds)\n",
    "print(kdp_con_centroid_lon)\n",
    "print(dist_kdp_zdr)\n",
    "shaped_dist = np.zeros((np.shape(kdp_con_centroid_lats)))\n",
    "print(shaped_dist.shape)\n",
    "shaped_dist[kdp_inds] = dist_kdp_zdr\n",
    "for i in kdp_inds:\n",
    "    print(i)\n",
    "    print(kdp_con_centroid_lon[i[0]])\n",
    "    print(kdp_con_centroid_lat[i[0]])\n",
    "    print(shaped_dist[i[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Calculate KDP-ZDR separation\n",
    "kdp_con_centroid_lons = np.asarray(kdp_con_centroid_lon)\n",
    "kdp_con_centroid_lats = np.asarray(kdp_con_centroid_lat)\n",
    "zdr_con_centroid_lons = np.asarray(zdr_con_centroid_lon)\n",
    "zdr_con_centroid_lats = np.asarray(zdr_con_centroid_lat)\n",
    "distance_kdp_zdr = g.inv(kdp_con_centroid_lons[kdp_inds], kdp_con_centroid_lats[kdp_inds], zdr_con_centroid_lons[kdp_inds], zdr_con_centroid_lats[kdp_inds])\n",
    "dist_kdp_zdr = distance_kdp_zdr[2] / 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kdp_inds = np.where(kdp_con_centroid_lats > 0)\n",
    "print(kdp_inds[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(zdr_storm_lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in enumerate(zdr_storm_lon):\n",
    "    print(i[0])\n",
    "    if i[0] != 0:\n",
    "        if zdr_storm_lon[i[0]-1] == zdr_storm_lon[i[0]]:\n",
    "            print(\"Skipping this one\")\n",
    "            #print(\"last lon:\",zdr_storm_lon[i[0]-1]\n",
    "            #print(\"this lon:\",zdr_storm_lon[i[0]])\n",
    "\n",
    "            continue\n",
    "        else:\n",
    "            print(\"using lon:\",zdr_storm_lon[i[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zdr_areas_arr = np.zeros((len(zdr_areas)))\n",
    "for i in range(len(zdr_areas)):\n",
    "    zdr_areas_arr[i] = zdr_areas[i].magnitude\n",
    "zdr_centroid_lons = np.asarray(zdr_centroid_lon)\n",
    "zdr_centroid_lats = np.asarray(zdr_centroid_lat)\n",
    "zdr_con_areas = []\n",
    "zdr_con_centroid_lon = []\n",
    "zdr_con_centroid_lat = []\n",
    "\n",
    "#For KDP as well\n",
    "kdp_areas_arr = np.zeros((len(kdp_areas)))\n",
    "for i in range(len(kdp_areas)):\n",
    "    kdp_areas_arr[i] = kdp_areas[i].magnitude\n",
    "kdp_centroid_lons = np.asarray(kdp_centroid_lon)\n",
    "kdp_centroid_lats = np.asarray(kdp_centroid_lat)\n",
    "kdp_con_areas = []\n",
    "kdp_con_centroid_lon = []\n",
    "kdp_con_centroid_lat = []\n",
    "dist_kdp_zdr = []\n",
    "for i in enumerate(zdr_storm_lon):\n",
    "    print(i[0])\n",
    "    if i[0] != 0:\n",
    "        if zdr_storm_lon[i[0]-1] == zdr_storm_lon[i[0]]:\n",
    "            print(\"Skipping this one\")\n",
    "            continue\n",
    "        else:\n",
    "            print(zdr_storm_lon[i[0]])\n",
    "            #Find the arc objects associated with this storm:\n",
    "            zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            print(\"zdr lons:\", zdr_objects_lons)\n",
    "            #Get the sum of their areas\n",
    "            print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            print(\"Areas sum:\", zdr_con_areas)\n",
    "            #Find the actual centroids\n",
    "            weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_con_centroid_lon.append(np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_con_centroid_lat.append(np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "            print(\"lons in loop\", zdr_objects_lons)\n",
    "\n",
    "            try:\n",
    "                #Find the kdp objects associated with this storm:\n",
    "                kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                print(\"kdp lons:\", kdp_objects_lons)\n",
    "                kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                #Get the sum of their areas\n",
    "                print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                kdp_con_areas.append(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                #Find the actual centroids\n",
    "                weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                kdp_con_centroid_lon.append(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                kdp_con_centroid_lat.append(np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                #Get the distance between the KDP and ZDR centroids\n",
    "                print(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "                distance_kdp_zdr = g.inv(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "                #print(distance_1[2]/1000)\n",
    "                print(\"separation:\", distance_kdp_zdr[2])\n",
    "                dist_kdp_zdr.append(distance_kdp_zdr[2]/1000.)\n",
    "            except:\n",
    "                print('storm missing kdp or zdr')\n",
    "    else:\n",
    "        print(zdr_storm_lon[i[0]])\n",
    "        #Find the arc objects associated with this storm:\n",
    "        zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "        zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "        print(\"zdr lons:\", zdr_objects_lons)\n",
    "        print(\"arc lats:\", zdr_objects_lats)\n",
    "        #Get the sum of their areas\n",
    "        print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "        zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "        print(\"Areas sum:\",zdr_con_areas)\n",
    "        #Find the actual centroids\n",
    "        weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "        zdr_con_centroid_lon.append(np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "        weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "        zdr_con_centroid_lat.append(np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "        print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "        print(\"lons out of loop\", zdr_objects_lons)\n",
    "        try:\n",
    "            #Find the kdp objects associated with this storm:\n",
    "            kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            print(\"kdp lons:\", kdp_objects_lons)\n",
    "            kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            #Get the sum of their areas\n",
    "            print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            kdp_con_areas.append(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            #Find the actual centroids\n",
    "            weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            kdp_con_centroid_lon.append(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            kdp_con_centroid_lat.append(np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            #Get the distance between the KDP and ZDR centroids\n",
    "            print(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "            distance_kdp_zdr = g.inv(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "            #print(distance_1[2]/1000)\n",
    "            print(\"separation:\", distance_kdp_zdr[2])\n",
    "            dist_kdp_zdr.append(distance_kdp_zdr[2]/1000.)\n",
    "        except:\n",
    "            print('storm missing kdp or zdr')\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kdp_con_centroid_lons = np.asarray(kdp_con_centroid_lon)\n",
    "kdp_con_centroid_lats = np.asarray(kdp_con_centroid_lat)\n",
    "zdr_con_centroid_lons = np.asarray(zdr_con_centroid_lon)\n",
    "zdr_con_centroid_lats = np.asarray(zdr_con_centroid_lat)\n",
    "kdp_inds = \n",
    "distance_kdp_zdr = g.inv(kdp_con_centroid_lons[kdp_inds], kdp_con_centroid_lats[kdp_inds], zdr_con_centroid_lons[kdp_inds], zdr_con_centroid_lats[kdp_inds])\n",
    "dist_kdp_zdr = distance_kdp_zdr[2] / 1000.\n",
    "print(dist_kdp_zdr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(zdr_storm_lon)\n",
    "print(kdp_storm_lon)\n",
    "print(zdr_areas)\n",
    "zdr_areas_arr = np.zeros((len(zdr_areas)))\n",
    "for i in range(len(zdr_areas)):\n",
    "    zdr_areas_arr[i] = zdr_areas[i].magnitude\n",
    "zdr_centroid_lons = np.asarray(zdr_centroid_lon)\n",
    "zdr_centroid_lats = np.asarray(zdr_centroid_lat)\n",
    "zdr_con_areas = []\n",
    "zdr_con_centroid_lon = []\n",
    "zdr_con_centroid_lat = []\n",
    "\n",
    "#For KDP as well\n",
    "kdp_areas_arr = np.zeros((len(kdp_areas)))\n",
    "for i in range(len(kdp_areas)):\n",
    "    kdp_areas_arr[i] = kdp_areas[i].magnitude\n",
    "kdp_centroid_lons = np.asarray(kdp_centroid_lon)\n",
    "kdp_centroid_lats = np.asarray(kdp_centroid_lat)\n",
    "kdp_con_areas = []\n",
    "kdp_con_centroid_lon = []\n",
    "kdp_con_centroid_lat = []\n",
    "dist_kdp_zdr = []\n",
    "for i in enumerate(zdr_storm_lon):\n",
    "    print(i[0])\n",
    "    if i[0] != 0:\n",
    "        if zdr_storm_lon[i[0]-1] == zdr_storm_lon[i[0]]:\n",
    "            print(\"Skipping this one\")\n",
    "            continue\n",
    "        else:\n",
    "            print(zdr_storm_lon[i[0]])\n",
    "            #Find the arc objects associated with this storm:\n",
    "            zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            print(\"zdr lons:\", zdr_objects_lons)\n",
    "            #Get the sum of their areas\n",
    "            print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            print(\"Areas sum:\", zdr_con_areas)\n",
    "            #Find the actual centroids\n",
    "            weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_con_centroid_lon = np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_con_centroid_lat = np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "            print(\"lons in loop\", zdr_objects_lons)\n",
    "            \n",
    "            try:\n",
    "                #Find the kdp objects associated with this storm:\n",
    "                kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                print(\"kdp lons:\", kdp_objects_lons)\n",
    "                kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                #Find the actual centroids\n",
    "                weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                kdp_con_centroid_lon = np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                kdp_con_centroid_lat = np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                #Get the distance between the KDP and ZDR centroids\n",
    "                print(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "                distance_kdp_zdr = g.inv(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "                #print(distance_1[2]/1000)\n",
    "                print(\"separation:\", distance_kdp_zdr[2])\n",
    "                dist_kdp_zdr.append(distance_kdp_zdr[2]/1000.)\n",
    "            except:\n",
    "                print('storm missing kdp or zdr')\n",
    "    else:\n",
    "            print(zdr_storm_lon[i[0]])\n",
    "            #Find the arc objects associated with this storm:\n",
    "            zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            print(\"zdr lons:\", zdr_objects_lons)\n",
    "            print(\"arc lats:\", zdr_objects_lats)\n",
    "            #Get the sum of their areas\n",
    "            print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "            print(\"Areas sum:\",zdr_con_areas)\n",
    "            #Find the actual centroids\n",
    "            weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_con_centroid_lon = np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "            zdr_con_centroid_lat = np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "            print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "            print(\"lons out of loop\", zdr_objects_lons)\n",
    "            try:\n",
    "                #Find the kdp objects associated with this storm:\n",
    "                kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                print(\"kdp lons:\", kdp_objects_lons)\n",
    "                kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                #Find the actual centroids\n",
    "                weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                kdp_con_centroid_lon = np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                kdp_con_centroid_lat = np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                #Get the distance between the KDP and ZDR centroids\n",
    "                print(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "                distance_kdp_zdr = g.inv(kdp_con_centroid_lon, kdp_con_centroid_lat, zdr_con_centroid_lon, zdr_con_centroid_lat)\n",
    "                #print(distance_1[2]/1000)\n",
    "                print(\"separation:\", distance_kdp_zdr[2])\n",
    "                dist_kdp_zdr.append(distance_kdp_zdr[2]/1000.)\n",
    "            except:\n",
    "                print('storm missing kdp or zdr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(dist_kdp_zdr)\n",
    "print(kdp_objects_lons)\n",
    "print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "print(weighted_lons_kdp)\n",
    "print(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "print(np.sum(weighted_lons_kdp))\n",
    "print(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#ZDRmasked3 = ma.masked_where(CCmasked < .85, ZDRmasked2)\n",
    "fig=plt.figure(2,figsize=(30.,25.))\n",
    "ax = plt.subplot(111,projection=ccrs.PlateCarree())\n",
    "ax.coastlines('50m',edgecolor='black',linewidth=0.75)\n",
    "#ax.add_feature(states_provinces,edgecolor='black',linewidth=0.5)\n",
    "#Add states/counties and set extent\n",
    "ax.add_feature(counties, edgecolor = 'black', linewidth = 0.5)\n",
    "ax.add_feature(states, edgecolor = 'black', linewidth = 1.5)\n",
    "ax.set_extent([LL[0],UR[0],LL[1],UR[1]])\n",
    "zdlevs = np.arange(1.75,8,1)\n",
    "#cf = ax.contourf(rlons[0,:,:],rlats[0,:,:],grad_ffd,grad_levs,cmap = plt.cm.gist_ncar)\n",
    "#ax.barbs(rlons[0,:,:],rlats[0,:,:],REFgradient[0,:,:],REFgradient[1,:,:],length=10,regrid_shape=12, color = 'k')\n",
    "ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked1,zdlevs,cmap = plt.cm.viridis)\n",
    "ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked2,zdlevs,cmap = plt.cm.plasma_r, alpha = .7)\n",
    "for i in enumerate(zdr_areas):\n",
    "        plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2\" %(zdr_areas[i[0]].magnitude), size = 23)\n",
    "if np.asarray(zdr_centroid_lon).shape[0] > 0:\n",
    "        ax.scatter(zdr_centroid_lon, zdr_centroid_lat, marker = '*', s = 500, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "#ax.scatter(zdr_con_centroid_lon, zdr_con_centroid_lat, marker = '*', s = 500, color = 'orange', zorder = 10, transform=ccrs.PlateCarree())\n",
    "if len(zdr_con_areas) > 0:\n",
    "    ax.scatter(zdr_con_centroid_lon, zdr_con_centroid_lat, marker = '*', s = 500, color = 'orange', zorder = 10, transform=ccrs.PlateCarree())\n",
    "    try:\n",
    "        plt.text(zdr_con_centroid_lon[i[0]]+.016, zdr_con_centroid_lat[i[0]]+.016, \"%.2f km^2\" %(zdr_con_areas[i[0]].magnitude), size = 23)\n",
    "    except:\n",
    "        print('failed')\n",
    "        plt.text(float(zdr_con_centroid_lon)+.016, float(zdr_con_centroid_lat)+.016, \"%.2f km^2\" %(float(zdr_con_areas[0])), size = 23)\n",
    "\n",
    "#plt.colorbar(cf)\n",
    "#plt.savefig(\"DDCbroughttoyoubytechnicolor1.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#ZDRmasked3 = ma.masked_where(CCmasked < .85, ZDRmasked2)\n",
    "fig=plt.figure(2,figsize=(30.,25.))\n",
    "ax = plt.subplot(111,projection=ccrs.PlateCarree())\n",
    "ax.coastlines('50m',edgecolor='black',linewidth=0.75)\n",
    "#ax.add_feature(states_provinces,edgecolor='black',linewidth=0.5)\n",
    "#Add states/counties and set extent\n",
    "ax.add_feature(counties, edgecolor = 'black', linewidth = 0.5)\n",
    "ax.add_feature(states, edgecolor = 'black', linewidth = 1.5)\n",
    "ax.set_extent([LL[0],UR[0],LL[1],UR[1]])\n",
    "zdlevs = np.arange(1.75,8,1)\n",
    "zdl = [4.0]\n",
    "kdlevs = np.arange(.75,3, .25)\n",
    "#cf = ax.contourf(rlons[0,:,:],rlats[0,:,:],grad_ffd,grad_levs,cmap = plt.cm.gist_ncar)\n",
    "#ax.barbs(rlons[0,:,:],rlats[0,:,:],REFgradient[0,:,:],REFgradient[1,:,:],length=10,regrid_shape=12, color = 'k')\n",
    "#ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked,zdlevs,cmap = plt.cm.viridis, zorder = 3, alpha = .76)\n",
    "#ax.contour(rlons[0,:,:],rlats[0,:,:],ZDRmasked,zdl, colors = 'red', linewidths = 3, zorder = 4)\n",
    "#ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,kdlevs,cmap = plt.cm.gist_ncar)\n",
    "ax.scatter(-94.4, 32.79, s=2000, marker = '^', color = 'y', zorder = 5)\n",
    "#ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked2,zdlevs,cmap = plt.cm.plasma_r, alpha = .7)\n",
    "zdrp = ax.pcolormesh(ungrid_lons, ungrid_lats, zdr_c, cmap=plt.cm.nipy_spectral, vmin = -2, vmax = 6)\n",
    "plt.colorbar(zdrp)\n",
    "#plt.savefig(\"DDCbroughttoyoubytechnicolor1.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(radar.fields)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "               pyproj.Proj(init='epsg:3857'))\n",
    "zdr_areas = []\n",
    "zdr_centroid_lon = []\n",
    "zdr_centroid_lat = []\n",
    "zdr_mean = []\n",
    "zdr_max = []\n",
    "\n",
    "for col in zdrc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "                print(new_shape)\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        boundary = np.asarray(poly.boundary.xy)\n",
    "        polypath = Path(boundary.transpose())\n",
    "        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "        mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "        mean = np.mean(ZDRmasked[mask])\n",
    "        if projected_area > 20 * units('km^2') and mean > 3.0:\n",
    "            zdr_areas.append((projected_area))\n",
    "            zdr_centroid_lon.append((poly.centroid.x))\n",
    "            zdr_centroid_lat.append((poly.centroid.y))\n",
    "            zdr_mean.append((mean))\n",
    "            zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "        #print(s_new)\n",
    "print(zdr_areas)\n",
    "print(zdr_centroid_lon)\n",
    "print(zdr_centroid_lat)\n",
    "print(zdr_mean)\n",
    "print(zdr_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(len(zdrc.collections)):\n",
    "    p = zdrc.collections[i].get_paths()[0]\n",
    "    v = p.vertices\n",
    "    x = v[:,0]\n",
    "    y = v[:,1]\n",
    "    if len(x)>2:\n",
    "        poly = sp.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "        print(i, poly.area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(len(zdrc.collections))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(enumerate(contour_path.to_polygons()))\n",
    "for cp in enumerate(contour_path.to_polygons()):\n",
    "    print(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "               pyproj.Proj(init='epsg:3857'))\n",
    "kdp_areas = []\n",
    "kdp_centroid_lon = []\n",
    "kdp_centroid_lat = []\n",
    "\n",
    "from shapely import geometry\n",
    "for col in kdpc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print(1)\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        if projected_area > 20 * units('km^2'):\n",
    "            kdp_areas.append((projected_area))\n",
    "            kdp_centroid_lon.append((poly.centroid.x))\n",
    "            kdp_centroid_lat.append((poly.centroid.y))\n",
    "\n",
    "        #print(s_new)\n",
    "print(kdp_areas)\n",
    "print(kdp_centroid_lon)\n",
    "print(kdp_centroid_lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "               pyproj.Proj(init='epsg:3857'))\n",
    "zdr_areas = []\n",
    "zdr_centroid_lon = []\n",
    "zdr_centroid_lat = []\n",
    "zdr_mean = []\n",
    "zdr_max = []\n",
    "\n",
    "for col in zdrc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    print('hi')\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        print('hi')\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print('hi')\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "                print('hi')\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "                print('hi')\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        #print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        boundary = np.asarray(poly.boundary.xy)\n",
    "        polypath = Path(boundary.transpose())\n",
    "        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "        mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "        mean = np.mean(ZDRmasked[mask])\n",
    "        if projected_area > 20 * units('km^2') and mean > 3.0:\n",
    "            zdr_areas.append((projected_area))\n",
    "            zdr_centroid_lon.append((poly.centroid.x))\n",
    "            zdr_centroid_lat.append((poly.centroid.y))\n",
    "            zdr_mean.append((mean))\n",
    "            zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "        #print(s_new)\n",
    "print(zdr_areas)\n",
    "print(zdr_centroid_lon)\n",
    "print(zdr_centroid_lat)\n",
    "print(zdr_mean)\n",
    "print(zdr_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "query = rs.query()\n",
    "#Here, set the initial time of the archived radar loop you want.\n",
    "dt = datetime(2017, 7, 11, 23, 30) # Our specified time\n",
    "station = 'KMVX'\n",
    "query.stations(station).time_range(dt, dt + timedelta(hours=.1))\n",
    "cat = rs.get_catalog(query)\n",
    "cat.datasets\n",
    "f = 27\n",
    "n = 1\n",
    "for item in sorted(cat.datasets.items()):\n",
    "    # After looping over the list of sorted datasets, pull the actual Dataset object out\n",
    "    # of our list of items and access over CDMRemote\n",
    "    try:\n",
    "        ds = item[1]\n",
    "        radar = pyart.io.nexrad_cdm.read_nexrad_cdm(ds.access_urls['OPENDAP'])\n",
    "    except:\n",
    "        continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "radar1 = radar.extract_sweeps([1])\n",
    "print(np.mean(radar1.elevation['data']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(radar.nsweeps)\n",
    "print(np.max(np.asarray(radar1.fields['differential_reflectivity']['data'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(radar.nsweeps):\n",
    "        radar2 = radar.extract_sweeps([i])\n",
    "        print(np.mean(radar2.elevation['data']))\n",
    "        if ((np.mean(radar2.elevation['data']) < .52) and (np.max(np.asarray(radar2.fields['differential_reflectivity']['data'])) != 1.0)):\n",
    "                print(np.mean(radar2.elevation['data']))\n",
    "                time_start = netCDF4.num2date(radar2.time['data'][0], radar2.time['units'])\n",
    "                print(time_start)\n",
    "                kdp_dict = pyart.retrieve.kdp_proc.kdp_maesaka(radar2)\n",
    "                print('kdp worked')\n",
    "        "
   ]
  },
  {
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}