Dual Pol Algorithm Working Copy.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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "in loop\n",
      "19\n",
      "0.524727\n",
      "2015-12-26 23:34:48.789000\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",
      "(-98.803055555555559, 31.073055555555555, <cartopy.crs.PlateCarree object at 0x0000021F8983E620>)\n",
      "plotting\n",
      "1\n",
      "-96.8812747155981\n",
      "past mask\n",
      "45.3401\n",
      "1\n",
      "-97.09931133932156\n",
      "past mask\n",
      "47.6567\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.80368420553583\n",
      "past mask\n",
      "46.459\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-97.82341500821965\n",
      "past mask\n",
      "45.4161\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "1\n",
      "-97.08435086794545\n",
      "[  6.16578444  54.92213553]\n",
      "1\n",
      "-97.06696926465543\n",
      "[  4.74433517  52.06813301]\n",
      "1\n",
      "-97.0973105691487\n",
      "[  0.35274701  50.13663533]\n",
      "1\n",
      "-97.1377901091425\n",
      "1\n",
      "-97.33104506612752\n",
      "1\n",
      "-96.80582236513114\n",
      "[ 49.11017592   1.53838995]\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "-96.8036842055\n",
      "[ 17.14834354]\n",
      "[ 11.44369472]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mwilson41\\AppData\\Local\\Continuum\\Anaconda3\\envs\\radar\\lib\\site-packages\\ipykernel_launcher.py:566: DeprecationWarning: assignment will raise an error in the future, most likely because your index result shape does not match the value array shape. You can use `arr.flat[index] = values` to keep the old behaviour.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "figure saved\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "0.529442\n",
      "2015-12-26 23:37:29.218000\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",
      "(-98.803055555555559, 31.073055555555555, <cartopy.crs.PlateCarree object at 0x0000021FDAB640F8>)\n",
      "plotting\n",
      "1\n",
      "-97.07897001478904\n",
      "past mask\n",
      "47.9916\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.78880385849048\n",
      "past mask\n",
      "46.0825\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-97.00575314317037\n",
      "past mask\n",
      "45.1375\n",
      "1\n",
      "-97.81194229435708\n",
      "past mask\n",
      "45.3443\n",
      "1\n",
      "-97.31549413341564\n",
      "past mask\n",
      "45.8873\n",
      "found a storm\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "1\n",
      "-97.07368038325687\n",
      "[   1.59625477   51.73642142  142.47360506]\n",
      "1\n",
      "-97.31235997025345\n",
      "1\n",
      "-96.78283485923819\n",
      "[  50.80391358    0.60972503  108.59338903]\n",
      "1\n",
      "-97.31103755339093\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "-96.7888038585\n",
      "[ 6.49542211]\n",
      "[ 7.46305585]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mwilson41\\AppData\\Local\\Continuum\\Anaconda3\\envs\\radar\\lib\\site-packages\\ipykernel_launcher.py:566: DeprecationWarning: assignment will raise an error in the future, most likely because your index result shape does not match the value array shape. You can use `arr.flat[index] = values` to keep the old behaviour.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "figure saved\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "0.51136\n",
      "2015-12-26 23:40:04.960000\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",
      "(-98.803055555555559, 31.073055555555555, <cartopy.crs.PlateCarree object at 0x0000021F8983E780>)\n",
      "plotting\n",
      "1\n",
      "-97.0619425855878\n",
      "past mask\n",
      "47.7816\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.78245265011705\n",
      "past mask\n",
      "45.9558\n",
      "1\n",
      "-96.98917428427791\n",
      "past mask\n",
      "45.633\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-97.8025981376634\n",
      "past mask\n",
      "45.3761\n",
      "1\n",
      "-97.3030775866881\n",
      "past mask\n",
      "45.5234\n",
      "found a storm\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "1\n",
      "-97.04861330871557\n",
      "[   3.29126671   62.40217603  146.35637879]\n",
      "1\n",
      "-97.08115671880486\n",
      "[   1.86799752   60.21107341  143.35928163]\n",
      "1\n",
      "-97.06781688546596\n",
      "1\n",
      "-96.78217662738858\n",
      "[  48.94498529   26.27514548  111.1721159 ]\n",
      "1\n",
      "-96.99810531515024\n",
      "made it through kdp centroids\n",
      "0\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n",
      "figure saved\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "0.529495\n",
      "2015-12-26 23:42:45.163000\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",
      "(-98.803055555555559, 31.073055555555555, <cartopy.crs.PlateCarree object at 0x0000021FD679C990>)\n",
      "plotting\n",
      "1\n",
      "-97.29277272941718\n",
      "past mask\n",
      "45.3465\n",
      "1\n",
      "-97.04521075736095\n",
      "past mask\n",
      "47.0505\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.76361460771906\n",
      "past mask\n",
      "46.2314\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.9734327142891\n",
      "past mask\n",
      "46.095\n",
      "found a storm\n",
      "added polygon\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "1\n",
      "-97.04556358846564\n",
      "[  3.66736611  52.95876109  63.71482088]\n",
      "1\n",
      "-97.07465052148878\n",
      "1\n",
      "-97.05369676850844\n",
      "[  0.96413694  49.76373392  59.62103532]\n",
      "1\n",
      "-97.37430275468071\n",
      "1\n",
      "-96.75928538039842\n",
      "[ 49.92605811   0.4168489   26.61648228]\n",
      "1\n",
      "-96.98419928500202\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "2\n",
      "-96.7636146077\n",
      "[ 3.30522554]\n",
      "[ 6.82846596]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mwilson41\\AppData\\Local\\Continuum\\Anaconda3\\envs\\radar\\lib\\site-packages\\ipykernel_launcher.py:566: DeprecationWarning: assignment will raise an error in the future, most likely because your index result shape does not match the value array shape. You can use `arr.flat[index] = values` to keep the old behaviour.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "figure saved\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "0.512291\n",
      "2015-12-26 23:45:21.300000\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",
      "(-98.803055555555559, 31.073055555555555, <cartopy.crs.PlateCarree object at 0x0000021FD678C150>)\n",
      "plotting\n",
      "1\n",
      "-97.02727821706574\n",
      "past mask\n",
      "47.5909\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.75106234810114\n",
      "past mask\n",
      "46.3169\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.96779988365333\n",
      "past mask\n",
      "46.453\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-97.79758308299844\n",
      "past mask\n",
      "45.0595\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "1\n",
      "-97.02381401892502\n",
      "[  3.49438284  51.83783579  64.31098214]\n",
      "1\n",
      "-97.04306321361463\n",
      "1\n",
      "-97.02747147750843\n",
      "[  1.037575    48.13997963  59.84136954]\n",
      "1\n",
      "-96.98514103737544\n",
      "[  4.39870175  45.34881444  58.71411237]\n",
      "1\n",
      "-97.01673100326894\n",
      "1\n",
      "-96.74526123694845\n",
      "[ 49.20441332   0.55344821  28.22798968]\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "2\n",
      "-96.7510623481\n",
      "[ 3.01648754]\n",
      "[ 10.35736734]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mwilson41\\AppData\\Local\\Continuum\\Anaconda3\\envs\\radar\\lib\\site-packages\\ipykernel_launcher.py:566: DeprecationWarning: assignment will raise an error in the future, most likely because your index result shape does not match the value array shape. You can use `arr.flat[index] = values` to keep the old behaviour.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "figure saved\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "0.528652\n",
      "2015-12-26 23:48:01.584000\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",
      "(-98.803055555555559, 31.073055555555555, <cartopy.crs.PlateCarree object at 0x0000021FABAFB728>)\n",
      "plotting\n",
      "1\n",
      "-97.01381030912015\n",
      "past mask\n",
      "47.6577\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.74907423852726\n",
      "past mask\n",
      "45.8764\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-96.96296994652306\n",
      "past mask\n",
      "46.2803\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-97.7861061040401\n",
      "past mask\n",
      "45.1121\n",
      "made it here\n",
      "made it through zdr centroids\n",
      "1\n",
      "-96.99617619854614\n",
      "[  2.95513946  49.70084354  64.28892588]\n",
      "1\n",
      "-97.03362997916497\n",
      "1\n",
      "-97.01881238090935\n",
      "[  1.30811735  47.57335373  60.76940405]\n",
      "1\n",
      "-96.99035500276212\n",
      "1\n",
      "-96.7775422075718\n",
      "1\n",
      "-96.747234627678\n",
      "[ 48.23600678   0.31316392  28.7305522 ]\n",
      "1\n",
      "-96.96529974167248\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "-96.7490742385\n",
      "[ 7.92308517  3.24081478]\n",
      "[ 6.36751769]\n",
      "2\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mwilson41\\AppData\\Local\\Continuum\\Anaconda3\\envs\\radar\\lib\\site-packages\\ipykernel_launcher.py:566: DeprecationWarning: assignment will raise an error in the future, most likely because your index result shape does not match the value array shape. You can use `arr.flat[index] = values` to keep the old behaviour.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "figure saved\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n",
      "in loop\n",
      "19\n"
     ]
    }
   ],
   "source": [
    "storm_relative_dir = 140\n",
    "srdir = 140\n",
    "query = rs.query()\n",
    "#Here, set the initial time of the archived radar loop you want.\n",
    "dt = datetime(2015,12, 26, 23, 30) # Our specified time\n",
    "station = 'KFWS'\n",
    "query.stations(station).time_range(dt, dt + timedelta(hours=.3))\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",
    "        radar1 = pyart.io.nexrad_cdm.read_nexrad_cdm(ds.access_urls['OPENDAP'])\n",
    "        #Now let's calculate and plot specific differential phase.\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']) < .60) 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",
    "                \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",
    "                                            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",
    "                                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.0]\n",
    "                kdplev = [.75]\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) < 15.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",
    "\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",
    "                    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",
    "                                #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",
    "                                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",
    "                                    \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",
    "                            #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",
    "                            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",
    "                    \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(kdp_con_centroid_lats)))\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(kdp_con_centroid_lats)))\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 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",
    "                \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 = 'red', linewidth = 5, transform=ccrs.PlateCarree())\n",
    "                except:\n",
    "                    print('KDP-ZDR separation didt work')\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 = n+1\n",
    "    except:\n",
    "        continue\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": {},
   "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
}