SLS ZDR Arc Algorithm NewNumpy.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\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\pyart\\graph\\cm.py:104: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if 'red' in spec:\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\pyart\\graph\\cm_colorblind.py:32: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if 'red' in spec:\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\babel\\plural.py:301: DeprecationWarning: Flags not at the start of the expression \\s+(?u)\n",
      "  (None, re.compile(r'\\s+(?u)')),\n"
     ]
    }
   ],
   "source": [
    "### This script will automatically detect ZDR arcs (and KDP feet) in WSR-88D radar data\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, get_wind_speed, get_wind_components\n",
    "import matplotlib.lines as mlines\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import csv\n",
    "import pickle\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "import nexradaws\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Bring in Homeyer rainbow colormap. This line will not be needed once the next version of PyART comes out, \n",
    "#but is commented out here since I'm defaulting it to run without it at the moment.\n",
    "#%run code_colormaps_CVD/colormap_generator.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Adding arrays containing the algorithm inputs for several tornadic supercell cases as examples\n",
    "#FFD angle, determined visually as the angle of the reflectivity gradient vector along the FFD\n",
    "storm_relative_dirstm = np.asarray([180, 200, 225, 190, 225, 190, 180, 190, 160, 170, 190, 165, 150])\n",
    "#ZDR value defining the ZDR arc 'core', here set to 3.25 dB\n",
    "zdrlevstm = np.asarray([3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25, 3.25])\n",
    "#KDP value defining the edge of the KDP foot. All sectiond dealing with KDP are still in the early stages of development\n",
    "kdplevstm = np.asarray([1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5])\n",
    "#First reflectivity contour used in the storm tracking algorithm\n",
    "REFlevstm = np.asarray([43, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45])\n",
    "#Second reflectivity contour used in the storm tracking algorithm\n",
    "REFlev1stm = np.asarray([48, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50])\n",
    "#Storm area threshold at which the tracking algorithm starts looking for more intense embedded cores\n",
    "big_stormstm = np.asarray([300, 300, 300, 300, 300, 300, 300, 600, 300, 300, 300, 300, 300])\n",
    "#Value to get rid of bad data files from THREDDS server, not really needed anymore\n",
    "zero_z_triggerstm = np.asarray([25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25])\n",
    "#Number of storm to track for earlier machine learning algorithm validation experiment, \n",
    "#not needed for just running the algorithm (just set to any integer)\n",
    "storm_to_trackstm = np.asarray([6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6])\n",
    "#Year of supercell case start time\n",
    "yearstm = np.asarray([2017, 2016, 2017, 2017, 2015, 2018, 2016, 2014, 2014, 2017, 2016, 2018, 2015])\n",
    "#Month of case start time\n",
    "monthstm = np.asarray([4, 8, 2, 7, 7, 6, 6, 6, 5, 5, 5, 5, 11])\n",
    "#Day of case start time\n",
    "daystm = np.asarray([2, 27, 7, 11, 14, 11, 13, 18, 11, 18, 24, 29, 16])\n",
    "#Hour of case start time (in UTC)\n",
    "hourstm = np.asarray([17, 19, 20, 22, 0, 21, 22, 2, 19, 18, 22, 20, 22])\n",
    "#Minute of case start time\n",
    "start_minstm = np.asarray([30, 30, 30, 0, 0, 30, 0, 30, 55, 17, 18, 35, 25])\n",
    "#Duration of analysis period, in hours\n",
    "durationstm = np.asarray([2.0, 2.1, 2.2, 3.3, 3.0, 1.7, 3.6, 2.6, 1.7, 1.9, 1.5, 2.1, 2.7])\n",
    "#WSR-88D radar site\n",
    "stationstm = ['KPOE', 'KMVX', 'KDGX', 'KMVX', 'KIND', 'KOAX', 'KAMA', 'KFSD', 'KUEX', 'KFDR', 'KDDC', 'KDDC', 'KAMA']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checklist of things to go through to make sure the algorithm will run on a particular computer:\n",
    "1. Specify a folder to direct the radar downloads to (line 34 in the cell below)\n",
    "2. Make sure the saved Random Forest file is in the directory this script is running in\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Actual code for the ZDR arc algorithm\n",
    "def multi_case_algorithm_ML1(storm_relative_dir, zdrlev, kdplev, REFlev, REFlev1, big_storm, zero_z_trigger, storm_to_track, year, month, day, hour, start_min, duration, station):    \n",
    "    #Settings\n",
    "    #Set vector perpendicular to FFD Z gradient\n",
    "    storm_relative_dir = storm_relative_dir\n",
    "    #Set ZDR Threshold for outlining arcs\n",
    "    zdrlev = [zdrlev]\n",
    "    #Set KDP Threshold for finding KDP feet\n",
    "    kdplev = [kdplev]\n",
    "    #Set reflectivity thresholds for storm tracking algorithm\n",
    "    REFlev = [REFlev]\n",
    "    REFlev1 = [REFlev1]\n",
    "    #Set storm size threshold that triggers subdivision of big storms\n",
    "    big_storm = big_storm #km^2\n",
    "    #Set search radii around storm centroids for ZDR arc objects. If fixed threshold is needed, uncomment these and comment\n",
    "    #out the lines where these variables are set in the algorithm\n",
    "    Outer_r = 30 #km\n",
    "    Inner_r = 6 #km\n",
    "    #Set trigger to ignore strangely-formatted files right before 00Z\n",
    "    #Pre-SAILS #: 17\n",
    "    #SAILS #: 25\n",
    "    zero_z_trigger = zero_z_trigger\n",
    "\n",
    "    storm_to_track = storm_to_track\n",
    "\n",
    "    #Here, set the initial time of the archived radar loop you want.\n",
    "    dt = datetime(year,month, day, hour, start_min) # Our specified time\n",
    "    station = station\n",
    "    end_dt = dt + timedelta(hours=duration)\n",
    "    #Set up nexrad interface\n",
    "    conn = nexradaws.NexradAwsInterface()\n",
    "    scans = conn.get_avail_scans_in_range(dt,end_dt,station)\n",
    "    #Change this to whatever folder you want the radar data to download into.\n",
    "    results = conn.download(scans, 'RadarFolder')\n",
    "    #query = rs.query()\n",
    "    #Set the duration of the loop in hours\n",
    "    #query.stations(station).time_range(dt, dt + timedelta(hours=duration))\n",
    "    #cat = rs.get_catalog(query)\n",
    "    #cat.datasets\n",
    "    \n",
    "    #Create an option for just reading all files in a folder\n",
    "    #folder = 'May27KDDC'\n",
    "\n",
    "    #Setting counters for figures and Pandas indices\n",
    "    f = 27\n",
    "    n = 1\n",
    "    storm_index = 0\n",
    "    scan_index = 0\n",
    "    #Create geod object for later distance and area calculations\n",
    "    g = Geod(ellps='sphere')\n",
    "    #Open the placefile\n",
    "    f = open(\"ARCexample\"+station+str(dt.year)+str(dt.month)+str(dt.day)+str(dt.hour)+str(dt.minute)+\"_Placefile.txt\", \"w+\")\n",
    "    f.write(\"Title: ZDR Arc Placefile \\n\")\n",
    "    f.write(\"Refresh: 8 \\n \\n\")\n",
    "    ####\n",
    "    #Load ML algorithm\n",
    "    #Make surt that you have this file in the directory that this script is in.\n",
    "    forest_loaded = pickle.load(open('BestRandomForest.pkl', 'rb'))\n",
    "    #Create file for ML algorithm test/training data\n",
    "    #Tornadic filename\n",
    "    #with open('Machine_Learning/ML_test'+station+str(dt.year)+str(dt.month)+str(dt.day)+str(dt.hour)+str(dt.minute)+'.csv', 'w') as csvfile:\n",
    "    #Nontornadic filename\n",
    "    #with open('Machine_Learning/NT2_ML_test'+station+str(dt.year)+str(dt.month)+str(dt.day)+str(dt.hour)+str(dt.minute)+'.csv', 'w') as csvfile:\n",
    "    #    fieldnames = ['number', 'hour', 'minute','area','distance','angle','mean','max','mean_cc','mean_kdp','mean_Z','mean_graddir','mean_grad', 'raw_angle']\n",
    "    #    writer = csv.DictWriter(csvfile, fieldnames=fieldnames)\n",
    "    #    writer.writeheader()\n",
    "    #print(datetime.utcnow())\n",
    "    ####\n",
    "    #Actual algorithm code starts here\n",
    "    #for item in sorted(cat.datasets.items()):\n",
    "    for i,scan in enumerate(results.iter_success(),start=1):\n",
    "    #Local file option:\n",
    "    #for radar_file in os.listdir(folder):\n",
    "        #print(radar_file)\n",
    "        #Loop over all files in the dataset and pull out each 0.5 degree tilt for analysis\n",
    "        try:\n",
    "            #ds = item[1]\n",
    "            radar1 = scan.open_pyart()\n",
    "            #Local file option\n",
    "            #radar1 = pyart.io.nexrad_archive.read_nexrad_archive(folder+'/'+radar_file)\n",
    "            #print('file read')\n",
    "            #Make sure the file isn't a strange format\n",
    "            if radar1.nsweeps > zero_z_trigger:\n",
    "                continue\n",
    "            for i in range(radar1.nsweeps):\n",
    "                print('in loop')\n",
    "                print(radar1.nsweeps)\n",
    "                radar = radar1.extract_sweeps([i])\n",
    "                #Checking to make sure the tilt in question has all needed data and is the right elevation\n",
    "                if ((np.mean(radar.elevation['data']) < .65) and (np.max(np.asarray(radar.fields['differential_reflectivity']['data'])) != np.min(np.asarray(radar.fields['differential_reflectivity']['data'])))):\n",
    "                    n = n+1\n",
    "                    print(np.mean(radar.elevation['data']))\n",
    "                    time_start = netCDF4.num2date(radar.time['data'][0], radar.time['units'])\n",
    "                    object_number=0.0\n",
    "                    print(time_start)\n",
    "                    month = time_start.month\n",
    "                    if month < 10:\n",
    "                        month = '0'+str(month)\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",
    "                    time_beg = time_start - timedelta(minutes=0.5)\n",
    "                    time_end = time_start + timedelta(minutes=0.5)\n",
    "                    sec_beg = time_beg.second\n",
    "                    sec_end = time_end.second\n",
    "                    min_beg = time_beg.minute\n",
    "                    min_end = time_end.minute\n",
    "                    h_beg = time_beg.hour\n",
    "                    h_end = time_end.hour\n",
    "                    d_beg = time_beg.day\n",
    "                    d_end = time_end.day\n",
    "                    if sec_beg < 10:\n",
    "                        sec_beg = '0'+str(sec_beg)\n",
    "                    if sec_end < 10:\n",
    "                        sec_end = '0'+str(sec_end)\n",
    "                    if min_beg < 10:\n",
    "                        min_beg = '0'+str(min_beg)\n",
    "                    if min_end < 10:\n",
    "                        min_end = '0'+str(min_end)\n",
    "                    if h_beg < 10:\n",
    "                        h_beg = '0'+str(h_beg)\n",
    "                    if h_end < 10:\n",
    "                        h_end = '0'+str(h_end)\n",
    "                    if d_beg < 10:\n",
    "                        d_beg = '0'+str(d_beg)\n",
    "                    if d_end < 10:\n",
    "                        d_end = '0'+str(d_end)\n",
    "                    #Add KDP to the dataset\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",
    "                    print(datetime.utcnow())\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",
    "                    #Mask out everything with reflectivity below Z=20 dBZ for Z and ZDR\n",
    "                    refl_c = np.copy(ref_ungridded)\n",
    "                    ref_c = ma.masked_where(refl_c < 20., refl_c)\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",
    "                    #Get ungridded lats/lons\n",
    "                    #kdp_ungridded = radar.fields['KDP']['data']\n",
    "                    phidp_ungridded = radar.fields['differential_phase']['data']\n",
    "                    cc_ungridded = radar.fields['cross_correlation_ratio']['data']\n",
    "                    print(datetime.utcnow())\n",
    "                    #NOW add KDP to the dataset\n",
    "                    #Now calculate KDP manually following NWS methodology\n",
    "                    #First, get the phidp gradient\n",
    "                    phidp_gradient = np.asarray(np.gradient(phidp_ungridded))/0.50\n",
    "                    kdp_raw = phidp_gradient[1,:,:]\n",
    "                    kdp1raw_c = np.copy(kdp_raw)\n",
    "                    kdpraw_c = ma.masked_where(kdp_raw > 40., kdp1raw_c)\n",
    "                    kdp_NWS = np.zeros((kdp_raw.shape[0], kdp_raw.shape[1]))\n",
    "                    #Do NWS smoothing process\n",
    "                    for i in range(kdp_raw.shape[0]):\n",
    "                        for j in range(kdp_raw.shape[1]):\n",
    "                            if cc_ungridded[i, j] > 0.90:\n",
    "                                #print(j)\n",
    "                                if ref_ungridded[i, j] > 40:\n",
    "                                    try:\n",
    "                                        kdp_new = np.mean(kdpraw_c[i, j-4:j+4])\n",
    "                                        kdp_NWS[i, j] = kdp_new\n",
    "                                        #print(np.mean(kdpraw_c[i, j-4:j+4]))\n",
    "                                    except:\n",
    "                                        kdp_NWS[i, j] = 0\n",
    "                                else:\n",
    "                                    try:\n",
    "                                        kdp_new = np.mean(kdpraw_c[i, j-12:j+12])\n",
    "                                        kdp_NWS[i, j] = kdp_new\n",
    "                                    except:\n",
    "                                        kdp_NWS[i, j] = 0  \n",
    "                    #Mask w/Z\n",
    "                    kdp1nws_c = np.copy(kdp_NWS)\n",
    "                    kdpnws_c = ma.masked_where(refl_c < 20., kdp1nws_c)\n",
    "                    #Create dictionary\n",
    "                    kdp_nwsdict = {}\n",
    "                    kdp_nwsdict['units'] = 'degrees/km'\n",
    "                    kdp_nwsdict['standard_name'] = 'specific_differential_phase_hv'\n",
    "                    kdp_nwsdict['long_name'] = 'Specific Differential Phase (KDP)'\n",
    "                    kdp_nwsdict['coordinates'] = 'elevation azimuth range'\n",
    "                    kdp_nwsdict['data'] = kdpnws_c\n",
    "                    kdp_nwsdict['valid_min'] = 0.0\n",
    "                    kdp_nwsdict['Clipf'] = 3906250000.0\n",
    "                    #Add field to radar\n",
    "                    radar.add_field('KDP', kdp_nwsdict)\n",
    "                    #Test data\n",
    "                    kdp_ungridded_nws = radar.fields['KDP']['data']\n",
    "                    ungrid_lons = radar.gate_longitude['data']\n",
    "                    ungrid_lats = radar.gate_latitude['data']\n",
    "                    print(datetime.utcnow())\n",
    "                    #Get ungridded gate altitudes\n",
    "                    gate_altitude = radar.gate_altitude['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 on a ~250 m x 250 m grid\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 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",
    "                    print(datetime.utcnow())\n",
    "                    #Mask everything below 20dbz on the grid\n",
    "                    ZDRmasked1 = ma.masked_where(REF < 20, ZDR)\n",
    "                    REFmasked = ma.masked_where(REF < 20, REF)\n",
    "                    #Use a 50 dBZ mask for KDP to only get areas in the storm core. This threshold should be considered more closely\n",
    "                    KDPmasked = ma.masked_where(REF < 50, KDP)\n",
    "                    KDPmasked = ma.filled(KDPmasked, fill_value = -2)\n",
    "                    #Try to filter out spots not in forward flank using Z gradient direction\n",
    "                    print('made it to smoothing')\n",
    "                    #First, smooth Z, take the gradient, and find its direction\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_mag = get_wind_speed(REFgradient[1,:,:] * units('m/s'), REFgradient[0,:,:] * units('m/s'))\n",
    "                    grad_dir = ma.masked_where(REF < 20, grad_dir1)\n",
    "                    #Get difference between the gradient direction and the FFD gradient direction calculated earlier\n",
    "                    srdir = storm_relative_dir\n",
    "                    srirad = np.copy(srdir)*units('degrees').to('radian')\n",
    "                    grad_dir = grad_dir*units('degrees').to('radian')\n",
    "                    grad_ffd = np.abs(np.arctan2(np.sin(grad_dir-srirad), np.cos(grad_dir-srirad)))\n",
    "                    print(grad_ffd)\n",
    "                    print(np.max(grad_ffd))\n",
    "                    print(np.min(grad_ffd))\n",
    "                    grad_ffd = np.asarray(grad_ffd)*units('radian')\n",
    "                    grad_ex = np.copy(grad_ffd)\n",
    "                    grad_ffd = grad_ffd.to('degrees')\n",
    "                    print(grad_ffd)\n",
    "                    print(np.max(grad_ffd))\n",
    "                    print(np.min(grad_ffd))\n",
    "                    #Mask out areas where the difference between the two is too large and the ZDR is likely not in the forward flank\n",
    "                    ZDRmasked2 = ma.masked_where(grad_ffd > 120 * units('degrees'), ZDRmasked1)\n",
    "                    ZDRmasked = ma.masked_where(CC < .60, ZDRmasked2)\n",
    "                    #Add a fill value for the ZDR mask so that contours will be closed\n",
    "                    ZDRmasked = ma.filled(ZDRmasked, fill_value = -2)\n",
    "                    #Extract the gridded lats and lons\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 our data 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",
    "                    #Make sure these shapefiles are in the same directory as the script\n",
    "                    #Lines with county/state shapefiles are commented out here, so you can add your own\n",
    "                    #or use the shapefiles provided on the github page\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",
    "                    #Create a figure and plot up the initial data and contours for the algorithm\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(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",
    "                    refp = ax.pcolormesh(ungrid_lons, ungrid_lats, ref_c, cmap=plt.cm.gist_ncar, vmin = 10, vmax = 73)\n",
    "                    #Homeyer rainbow colormap commented out for now, but can be added back in if you have it\n",
    "                    #refp = ax.pcolormesh(ungrid_lons, ungrid_lats, ref_c, cmap='HomeyerRainbow', vmin = 10, vmax = 73)\n",
    "                    #Option to have a ZDR background instead of Z:\n",
    "                    #zdrp = ax.pcolormesh(ungrid_lons, ungrid_lats, zdr_c, cmap=plt.cm.nipy_spectral, vmin = -2, vmax = 6)\n",
    "                    ##\n",
    "                    #Storm tracking algorithm starts here\n",
    "                    ##\n",
    "                    #Reflectivity smoothed for storm tracker\n",
    "                    smoothed_ref = ndi.gaussian_filter(REFmasked, sigma = 3, order = 0)\n",
    "                    #REFlev = [45]\n",
    "                    #REFlev1 = [50]\n",
    "                    #1st Z contour plotted\n",
    "                    refc = ax.contour(rlons[0,:,:],rlats[0,:,:],smoothed_ref,REFlev, alpha=.4)\n",
    "                    #Empty arrays made for storm characteristics\n",
    "                    ref_areas = []\n",
    "                    max_lons_c = []\n",
    "                    max_lats_c = []\n",
    "                    storm_ids = []\n",
    "                    #Set up projection for area calculations\n",
    "                    proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "                               pyproj.Proj(init='epsg:3857'))\n",
    "\n",
    "                    #Main part of storm tracking algorithm starts by looping through all contours looking for Z centroids\n",
    "                    #This method for breaking contours into polygons based on this stack overflow tutorial:\n",
    "                    #https://gis.stackexchange.com/questions/99917/converting-matplotlib-contour-objects-to-shapely-objects\n",
    "                    for level in refc.collections:\n",
    "                        #Loops through each closed polygon in the contour \n",
    "                        for contour_poly in level.get_paths(): \n",
    "                            for n_contour,contour in enumerate(contour_poly.to_polygons()):\n",
    "                                print(1)\n",
    "                                contour_a = np.asarray(contour[:])\n",
    "                                xa = contour_a[:,0]\n",
    "                                ya = contour_a[:,1]\n",
    "                                polygon_new = geometry.Polygon([(i[0], i[1]) for i in zip(xa,ya)])\n",
    "                                #Eliminates 'holes' in the polygons\n",
    "                                if n_contour == 0:\n",
    "                                    polygon = polygon_new\n",
    "                                else:\n",
    "                                    polygon = polygon.difference(polygon_new)\n",
    "\n",
    "                            print(polygon.centroid.x)\n",
    "                            #Transform the polygon's coordinates to the proper projection and calculate area\n",
    "                            pr_area = (transform(proj, polygon).area * units('m^2')).to('km^2')\n",
    "                            #Use the polygon boundary to select all points within the polygon via a mask\n",
    "                            boundary = np.asarray(polygon.boundary.xy)\n",
    "                            polypath = Path(boundary.transpose())\n",
    "                            coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T \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 pr_area > 10 * units('km^2') and meanr > REFlev[0]:\n",
    "                                print('found a storm')\n",
    "                                #For big blobs with embedded supercells, find the embedded storm cores\n",
    "                                #Normal 'big storm' cutoff 300 km^2\n",
    "                                if pr_area > big_storm * 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",
    "                                    #This section uses the 2nd reflectivity threshold to subdivide big storms, in a similar manner to the \n",
    "                                    #previous section's method for finding storms\n",
    "                                    refc1 = ax.contour(rlon2m,rlat2m,smoothed_ref_m,REFlev1, linewidths = 3, linestyle = '--', alpha=.4)\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 level1 in refc1.collections:\n",
    "                                        print('made it to beginning of loop')\n",
    "                                        for contour_poly1 in level1.get_paths(): \n",
    "                                            for n_contour1,contour1 in enumerate(contour_poly1.to_polygons()):\n",
    "                                                print(2)\n",
    "                                                contour_a1 = np.asarray(contour1[:])\n",
    "                                                xa1 = contour_a1[:,0]\n",
    "                                                ya1 = contour_a1[:,1]\n",
    "                                                polygon_new1 = geometry.Polygon([(i[0], i[1]) for i in zip(xa1,ya1)])\n",
    "                                                if n_contour1 == 0:\n",
    "                                                    polygon1 = polygon_new1\n",
    "                                                else:\n",
    "                                                    polygon1 = polygon1.difference(polygon_new1)\n",
    "\n",
    "                                            print(polygon1.centroid.x)\n",
    "                                            pr_area1 = (transform(proj, polygon1).area * units('m^2')).to('km^2')\n",
    "                                            boundary1 = np.asarray(polygon1.boundary.xy)\n",
    "                                            polypath1 = Path(boundary1.transpose())\n",
    "                                            maskr1 = polypath1.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "                                            meanr1 = np.mean(smoothed_ref[maskr1])\n",
    "                                            #Add objects that fit requirements to the list of storm objects\n",
    "                                            if pr_area1 > 10 * units('km^2') and meanr1 > REFlev1[0]:\n",
    "                                                ref_areas.append((pr_area1.magnitude*2))\n",
    "                                                max_lons_c.append((polygon1.centroid.x))\n",
    "                                                max_lats_c.append((polygon1.centroid.y))\n",
    "                                                #For tracking, assign ID numbers and match current storms to any previous storms that are close enough to \n",
    "                                                #be the same\n",
    "                                                if scan_index == 0:\n",
    "                                                    storm_ids.append((storm_index))\n",
    "                                                    storm_index = storm_index + 1\n",
    "                                                else:\n",
    "                                                    #dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                                    max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "                                                    max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "                                                    storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "                                                    dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                                    for i in range(max_lons_p.shape[0]):\n",
    "                                                        distance_track = g.inv(polygon1.centroid.x, polygon1.centroid.y,\n",
    "                                                                               max_lons_p[i], max_lats_p[i])\n",
    "                                                        dist_track[i] = distance_track[2]/1000.\n",
    "                                                    print(dist_track)\n",
    "                                                    print('Poly lon', polygon1.centroid.x)\n",
    "                                                    print(max_lons_p)\n",
    "                                                    print(storm_ids_p)\n",
    "                                                    if np.min(dist_track) < 10.0:\n",
    "                                                        storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "                                                        print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                                    else:\n",
    "                                                        storm_ids.append((storm_index))\n",
    "                                                        print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                                        storm_index = storm_index + 1\n",
    "                                                print('added polygon')\n",
    "                                            else:\n",
    "                                                print('nope')\n",
    "                                #Do the same thing for objects from the 1st reflectivity threshold\n",
    "                                else:\n",
    "                                    ref_areas.append((pr_area.magnitude))\n",
    "                                    max_lons_c.append((polygon.centroid.x))\n",
    "                                    max_lats_c.append((polygon.centroid.y))\n",
    "                                    if scan_index == 0:\n",
    "                                        storm_ids.append((storm_index))\n",
    "                                        storm_index = storm_index + 1\n",
    "                                    else:\n",
    "                                        #dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                        max_lons_p = np.asarray(tracks_dataframe['storm_lon'].loc[scan_index-1].iloc[:])\n",
    "                                        max_lats_p = np.asarray(tracks_dataframe['storm_lat'].loc[scan_index-1].iloc[:])\n",
    "                                        storm_ids_p = np.asarray(tracks_dataframe['storm_id1'].loc[scan_index-1].iloc[:])\n",
    "                                        dist_track = np.zeros((np.asarray(max_lons_p).shape[0]))\n",
    "                                        for i in range(max_lons_p.shape[0]):\n",
    "                                            distance_track = g.inv(polygon.centroid.x, polygon.centroid.y,\n",
    "                                                                   max_lons_p[i], max_lats_p[i])\n",
    "                                            dist_track[i] = distance_track[2]/1000.\n",
    "                                        print(dist_track)\n",
    "                                        print('Poly lon', polygon.centroid.x)\n",
    "                                        print(max_lons_p)\n",
    "                                        print(storm_ids_p)\n",
    "                                        if np.min(dist_track) < 10.0:\n",
    "                                            storm_ids.append((storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]]))\n",
    "                                            print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                        else:\n",
    "                                            storm_ids.append((storm_index))\n",
    "                                            print('storm id', storm_ids_p[np.where(dist_track == np.min(dist_track))[0][0]])\n",
    "                                            storm_index = storm_index + 1\n",
    "                                    print('added polygon')\n",
    "\n",
    "                        #print(s_new)\n",
    "                    #Setup tracking index for storm of interest\n",
    "                    tracking_ind=np.where(np.asarray(storm_ids)==storm_to_track)[0]\n",
    "                    print('tracking id')\n",
    "                    print(tracking_ind)\n",
    "                    max_lons_c = np.asarray(max_lons_c)\n",
    "                    max_lats_c = np.asarray(max_lats_c)\n",
    "                    ref_areas = np.asarray(ref_areas)\n",
    "                    #Create the ZDR and KDP contours which will later be broken into polygons\n",
    "                    zdrc = ax.contour(rlons[0,:,:],rlats[0,:,:],ZDRmasked,zdrlev,linewidths = 2, colors='purple', alpha = .7)\n",
    "                    kdpc = ax.contour(rlons[0,:,:],rlats[0,:,:],KDPmasked,kdplev,linewidths = 2, colors='green', alpha = .8)\n",
    "\n",
    "                    print('made it here')\n",
    "                    plt.savefig('testfig.png')\n",
    "\n",
    "                    #Create ZDR arc objects using a similar method as employed in making the storm objects\n",
    "                    if len(max_lons_c) > 0:\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",
    "                        zdr_masks = []\n",
    "                        #print(\"here too\")\n",
    "                        #Break contours into polygons using the same method as for reflectivity\n",
    "                        for level in zdrc.collections:\n",
    "                            for contour_poly in level.get_paths(): \n",
    "                                for n_contour,contour in enumerate(contour_poly.to_polygons()):\n",
    "                                    #print('hi')\n",
    "                                    contour_a = np.asarray(contour[:])\n",
    "                                    xa = contour_a[:,0]\n",
    "                                    ya = contour_a[:,1]\n",
    "                                    polygon_new = geometry.Polygon([(i[0], i[1]) for i in zip(xa,ya)])\n",
    "                                    if n_contour == 0:\n",
    "                                        polygon = polygon_new\n",
    "                                        #print('hi')\n",
    "                                    else:\n",
    "                                        polygon = polygon.difference(polygon_new)\n",
    "                                        #print('hi')\n",
    "                                pr_area = (transform(proj, polygon).area * units('m^2')).to('km^2')\n",
    "                                boundary = np.asarray(polygon.boundary.xy)\n",
    "                                polypath = Path(boundary.transpose())\n",
    "                                coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T \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",
    "                                mean_Z = np.mean(REF[mask])\n",
    "                                mean_graddir = np.mean(grad_ffd[mask])\n",
    "                                mean_grad = np.mean(grad_mag[mask])\n",
    "                                mean_kdp = np.mean(KDP[mask])\n",
    "                                #Only select out objects larger than 1 km^2 with high enough CC\n",
    "                                if pr_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",
    "                                    rawangle = np.zeros((np.asarray(max_lons_c).shape[0]))\n",
    "                                    back = np.zeros((np.asarray(max_lons_c).shape[0]))\n",
    "                                    zdr_polypath = polypath\n",
    "                                    #Assign ZDR arc objects to the nearest acceptable storm object\n",
    "                                    for i in range(dist.shape[0]):\n",
    "                                                distance_1 = g.inv(polygon.centroid.x, polygon.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",
    "                                                #print('Raw back angle', back[i])\n",
    "                                                if distance_1[1] < 0:\n",
    "                                                    back[i] = distance_1[1] + 360\n",
    "                                                #print('fixed back', back[i])\n",
    "                                                forw[i] = np.abs(back[i] - storm_relative_dir)\n",
    "                                                #print('raw forw', forw[i])\n",
    "                                                rawangle[i] = back[i] - storm_relative_dir\n",
    "                                                #Account for weird angles\n",
    "                                                if forw[i] > 180:\n",
    "                                                    #print('Big angle')\n",
    "                                                    forw[i] = 360 - forw[i]\n",
    "                                                    rawangle[i] = (360-forw[i])*(-1)\n",
    "                                                    #print(rawangle[i])\n",
    "                                                dist[i] = distance_1[2]/1000.\n",
    "                                                rawangle[i] = rawangle[i]*(-1)\n",
    "                                                #print('fixed forw', forw[i])\n",
    "                                                #print('raw angle', rawangle[i])\n",
    "                                    #print(dist.shape)\n",
    "                                    #Set search radii around storms for ZDR arcs objects\n",
    "                                    #Outer Radius-variable, 5km past the square root of the storm object's area\n",
    "                                    #Outer_r = (1.0*np.sqrt(ref_areas[np.where(dist == np.min(dist))[0][0]]))\n",
    "                                    #Outer_r = 15.0\n",
    "                                    #Inner Radius-variable, 2km pask the square root of 1/4 of the storm's area\n",
    "                                    #Inner_r = (0.25*np.sqrt(ref_areas[np.where(dist == np.min(dist))[0][0]]))\n",
    "                                    #Inner_r = 6.0\n",
    "                                    #Pick out only ZDR arc objects with a reasonable probability of actually being in the FFD region\n",
    "                                    #using their location relative to the storm centroid\n",
    "                                    if (forw[np.where(dist == np.min(dist))[0][0]] < 180 and np.min(dist) < Outer_r) or (forw[np.where(dist == np.min(dist))[0][0]] < 140 and np.min(dist) < Inner_r):\n",
    "                                        #Use ML algorithm to eliminate non-arc objects\n",
    "                                        #Get x and y components\n",
    "                                        if (rawangle[np.where(dist == np.min(dist))[0][0]] > 0):\n",
    "                                            directions_raw = 360 - rawangle[np.where(dist == np.min(dist))[0][0]]\n",
    "                                        else:\n",
    "                                            directions_raw = (-1) * rawangle[np.where(dist == np.min(dist))[0][0]]\n",
    "                                        \n",
    "                                        xc, yc = get_wind_components(np.min(dist), directions_raw * units('degree'))\n",
    "                                        print('got xc')\n",
    "                                        ARC_X = np.zeros((1, 12))\n",
    "                                        print('got array1')\n",
    "                                        ARC_X[:,0] = pr_area.magnitude\n",
    "                                        print('got array2')\n",
    "                                        ARC_X[:,1] = np.min(dist)\n",
    "                                        print('got array3')\n",
    "                                        ARC_X[:,2] = np.max(ZDRmasked[mask]) / mean\n",
    "                                        print('got array4')\n",
    "                                        ARC_X[:,3] = (np.max(ZDRmasked[mask]) / mean) * pr_area.magnitude\n",
    "                                        print('got array5')\n",
    "                                        ARC_X[:,4] = mean_cc\n",
    "                                        print('got array6')\n",
    "                                        ARC_X[:,5] = mean_kdp\n",
    "                                        print('got array7')\n",
    "                                        ARC_X[:,6] = mean_Z\n",
    "                                        print('got array8')\n",
    "                                        ARC_X[:,7] = mean_graddir\n",
    "                                        print('got array9')\n",
    "                                        ARC_X[:,8] = mean_grad\n",
    "                                        print('got array10')\n",
    "                                        ARC_X[:,9] = rawangle[np.where(dist == np.min(dist))[0][0]]\n",
    "                                        print('got array11')\n",
    "                                        ARC_X[:,10] = xc\n",
    "                                        print('got array12')\n",
    "                                        ARC_X[:,11] = yc\n",
    "                                        print('got to prediction')\n",
    "                                        pred_zdr = forest_loaded.predict(ARC_X)\n",
    "                                        print(pred_zdr)\n",
    "                                        if (pred_zdr[0]==1):\n",
    "                                            print('arc')\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((pr_area))\n",
    "                                            zdr_centroid_lon.append((polygon.centroid.x))\n",
    "                                            zdr_centroid_lat.append((polygon.centroid.y))\n",
    "                                            zdr_mean.append((mean))\n",
    "                                            zdr_cc_mean.append((mean_cc))\n",
    "                                            zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "                                            zdr_masks.append(mask)\n",
    "                                            patch = PathPatch(polypath, facecolor='none', alpha=.8, edgecolor = 'blue', linewidth = 3)\n",
    "                                            ax.add_patch(patch)\n",
    "                                            #Add polygon to placefile\n",
    "                                            f.write('TimeRange: '+str(time_start.year)+'-'+str(month)+'-'+str(d_beg)+'T'+str(h_beg)+':'+str(min_beg)+':'+str(sec_beg)+'Z '+str(time_start.year)+'-'+str(month)+'-'+str(d_end)+'T'+str(h_end)+':'+str(min_end)+':'+str(sec_end)+'Z')\n",
    "                                            f.write('\\n')\n",
    "                                            f.write(\"Color: 000 000 139 \\n\")\n",
    "                                            f.write('Line: 3, 0, \"ZDR Arc Outline\" \\n')\n",
    "                                            for i in range(len(zdr_polypath.vertices)):\n",
    "                                                f.write(\"%.5f\" %(zdr_polypath.vertices[i][1]))\n",
    "                                                f.write(\", \")\n",
    "                                                f.write(\"%.5f\" %(zdr_polypath.vertices[i][0]))\n",
    "                                                f.write('\\n')\n",
    "                                            #f.write(str(zdr_polypath.vertices[0][0])+', '+str(zdr_polypath.vertices[0][1])+'\\n')\n",
    "                                            f.write(\"End: \\n \\n\")\n",
    "                                        #if (((max_lons_c[np.where(dist == np.min(dist))[0][0]]) in max_lons_c[tracking_ind]) and ((max_lats_c[np.where(dist == np.min(dist))[0][0]]) in max_lats_c[tracking_ind])):\n",
    "                                        #    plt.text(float(polygon.centroid.x), float(polygon.centroid.y), \"%.1f\" %(float(object_number)), size = 23, color = 'purple')\n",
    "                                            #Add a line that writes all of the attibutes of each object to a csv\n",
    "                                            #Tornadic filename\n",
    "                                            #with open('Machine_Learning/ML_test'+station+str(dt.year)+str(dt.month)+str(dt.day)+str(dt.hour)+str(dt.minute)+'.csv', 'a') as csvfile:\n",
    "                                            #Nontornadic filename\n",
    "                                           # with open('Machine_Learning/NT2_ML_test'+station+str(dt.year)+str(dt.month)+str(dt.day)+str(dt.hour)+str(dt.minute)+'.csv', 'a') as csvfile:\n",
    "                                            #    fieldnames = ['number', 'hour', 'minute','area','distance','angle','mean','max','mean_cc','mean_kdp','mean_Z','mean_graddir','mean_grad', 'raw_angle']\n",
    "                                            #    writer = csv.DictWriter(csvfile, fieldnames=fieldnames)\n",
    "                                            #    writer.writerow({'number': object_number, 'hour': hour, 'minute': minute, 'area': pr_area.magnitude, 'distance': np.min(dist), 'angle': forw[np.where(dist == np.min(dist))[0][0]], 'mean': mean, 'max': np.max(ZDRmasked[mask]), 'mean_cc': mean_cc, 'mean_kdp': mean_kdp, 'mean_Z': mean_Z, 'mean_graddir': mean_graddir.magnitude, 'mean_grad': mean_grad.magnitude, 'raw_angle': rawangle[np.where(dist == np.min(dist))[0][0]]})\n",
    "                                            #object_number=object_number+1\n",
    "                                                     #print(s_new)\n",
    "                        print('made it through zdr centroids')\n",
    "                        #Identify KDP foot objects in a similar way to the ZDR arc objects\n",
    "                        if len(max_lons_c) > 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 level in kdpc.collections:\n",
    "                                for contour_poly in level.get_paths(): \n",
    "                                    for n_contour,contour in enumerate(contour_poly.to_polygons()):\n",
    "                                        print(1)\n",
    "                                        contour_a = np.asarray(contour[:])\n",
    "                                        xa = contour_a[:,0]\n",
    "                                        ya = contour_a[:,1]\n",
    "                                        polygon_new = geometry.Polygon([(i[0], i[1]) for i in zip(xa,ya)])\n",
    "                                        if n_contour == 0:\n",
    "                                            polygon = polygon_new\n",
    "                                        else:\n",
    "                                            polygon = polygon.difference(polygon_new)\n",
    "\n",
    "                                    print(polygon.centroid.x)\n",
    "                                    pr_area = (transform(proj, polygon).area * units('m^2')).to('km^2')\n",
    "                                    boundary = np.asarray(polygon.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 pr_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(polygon.centroid.x, polygon.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_path = polypath\n",
    "                                            kdp_areas.append((pr_area))\n",
    "                                            kdp_centroid_lon.append((polygon.centroid.x))\n",
    "                                            kdp_centroid_lat.append((polygon.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",
    "                                            #Add polygon to placefile\n",
    "                                            #f.write('TimeRange: '+str(time_start.year)+'-'+str(month)+'-'+str(d_beg)+'T'+str(h_beg)+':'+str(min_beg)+':00Z '+str(time_start.year)+'-'+str(month)+'-'+str(d_end)+'T'+str(h_end)+':'+str(min_end)+':00Z')\n",
    "                                            f.write('TimeRange: '+str(time_start.year)+'-'+str(month)+'-'+str(d_beg)+'T'+str(h_beg)+':'+str(min_beg)+':'+str(sec_beg)+'Z '+str(time_start.year)+'-'+str(month)+'-'+str(d_end)+'T'+str(h_end)+':'+str(min_end)+':'+str(sec_end)+'Z')\n",
    "                                            f.write('\\n')\n",
    "                                            f.write(\"Color: 000 139 000 \\n\")\n",
    "                                            f.write('Line: 3, 0, \"KDP Foot Outline\" \\n')\n",
    "                                            for i in range(len(kdp_path.vertices)):\n",
    "                                                f.write(\"%.5f\" %(kdp_path.vertices[i][1]))\n",
    "                                                f.write(\", \")\n",
    "                                                f.write(\"%.5f\" %(kdp_path.vertices[i][0]))\n",
    "                                                f.write('\\n')\n",
    "                                            #f.write(str(zdr_polypath.vertices[0][0])+', '+str(zdr_polypath.vertices[0][1])+'\\n')\n",
    "                                            f.write(\"End: \\n \\n\")\n",
    "\n",
    "                            print('made it through kdp centroids')\n",
    "\n",
    "                            #Consolidating the arc objects associated with each storm:\n",
    "                            zdr_areas_arr = np.zeros((len(zdr_areas)))\n",
    "                            zdr_max_arr = np.zeros((len(zdr_max)))\n",
    "                            zdr_mean_arr = np.zeros((len(zdr_mean)))                    \n",
    "                            for i in range(len(zdr_areas)):\n",
    "                                zdr_areas_arr[i] = zdr_areas[i].magnitude\n",
    "                                zdr_max_arr[i] = zdr_max[i]\n",
    "                                zdr_mean_arr[i] = zdr_mean[i]\n",
    "\n",
    "                            zdr_centroid_lons = np.asarray(zdr_centroid_lon)\n",
    "                            zdr_centroid_lats = np.asarray(zdr_centroid_lat)\n",
    "                            zdr_con_areas = []\n",
    "                            zdr_con_maxes = []\n",
    "                            zdr_con_means = []\n",
    "                            zdr_con_centroid_lon = []\n",
    "                            zdr_con_centroid_lat = []\n",
    "                            zdr_con_max_lon = []\n",
    "                            zdr_con_max_lat = []\n",
    "                            zdr_con_storm_lon = []\n",
    "                            zdr_con_storm_lat = []\n",
    "                            zdr_con_masks = []\n",
    "                            zdr_con_dev = []\n",
    "                            zdr_con_10max = []\n",
    "                            zdr_con_mode = []\n",
    "                            zdr_con_median = []\n",
    "                            zdr_masks = np.asarray(zdr_masks)\n",
    "                            #Consolidate KDP objects as well\n",
    "                            kdp_areas_arr = np.zeros((len(kdp_areas)))\n",
    "                            kdp_max_arr = np.zeros((len(kdp_max)))\n",
    "                            for i in range(len(kdp_areas)):\n",
    "                                kdp_areas_arr[i] = kdp_areas[i].magnitude\n",
    "                                kdp_max_arr[i] = kdp_max[i]\n",
    "                            kdp_centroid_lons = np.asarray(kdp_centroid_lon)\n",
    "                            kdp_centroid_lats = np.asarray(kdp_centroid_lat)\n",
    "                            kdp_con_areas = []\n",
    "                            kdp_con_maxes = []\n",
    "                            kdp_con_centroid_lon = []\n",
    "                            kdp_con_centroid_lat = []\n",
    "                            kdp_con_max_lon = []\n",
    "                            kdp_con_max_lat = []\n",
    "                            kdp_con_storm_lon = []\n",
    "                            kdp_con_storm_lat = []\n",
    "                            for i in enumerate(zdr_storm_lon):\n",
    "                                print(i[0])\n",
    "                                if i[0] != 0:\n",
    "                                    if zdr_storm_lon[i[0]-1] == zdr_storm_lon[i[0]]:\n",
    "                                        #print(\"Skipping this one\")\n",
    "                                        continue\n",
    "                                    else:\n",
    "                                        print(zdr_storm_lon[i[0]])\n",
    "                                        #Find the arc objects associated with this storm:\n",
    "                                        zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        #print(\"zdr lons:\", zdr_objects_lons)\n",
    "                                        #Get the sum of their areas\n",
    "                                        print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                        zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        zdr_con_maxes.append(np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        zdr_con_means.append(np.mean(zdr_mean_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        zdr_con_max_lon.append(rlons_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                        zdr_con_max_lat.append(rlats_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                        #print(\"Areas sum:\", zdr_con_areas)\n",
    "                                        #Find the actual centroids\n",
    "                                        weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        zdr_con_centroid_lon.append(np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        zdr_con_centroid_lat.append(np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        zdr_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                        zdr_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                                        zdr_con_masks.append(np.sum(zdr_masks[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])],axis=0, dtype=bool))\n",
    "                                        mask_con = np.sum(zdr_masks[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])], axis=0, dtype=bool)\n",
    "                                        zdr_con_dev.append(np.std(ZDRmasked[mask_con]))\n",
    "                                        ZDRsorted = np.sort(ZDRmasked[mask_con])[::-1]\n",
    "                                        zdr_con_10max.append(np.mean(ZDRsorted[0:10]))\n",
    "                                        zdr_con_mode.append(stats.mode(ZDRmasked[mask_con]))\n",
    "                                        zdr_con_median.append(np.median(ZDRmasked[mask_con]))\n",
    "                                        #print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "                                        #print(\"lons in loop\", zdr_objects_lons)\n",
    "\n",
    "                                        try:\n",
    "                                            #Find the kdp objects associated with this storm:\n",
    "                                            kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                            #print(\"kdp lons:\", kdp_objects_lons)\n",
    "                                            kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                            #Get the sum of their areas\n",
    "                                            print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                            kdp_con_areas.append(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                            kdp_con_maxes.append(np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                            kdp_con_max_lon.append(rlons_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                            kdp_con_max_lat.append(rlats_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                            #Find the actual centroids\n",
    "                                            weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                            kdp_con_centroid_lon.append(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                            weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                            #print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                            kdp_con_centroid_lat.append(np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                            kdp_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                            kdp_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                                        except:\n",
    "                                            print('storm missing kdp or zdr')\n",
    "                                            kdp_con_areas.append(0)\n",
    "                                            kdp_con_maxes.append(0)\n",
    "                                            kdp_con_max_lon.append(0)\n",
    "                                            kdp_con_max_lat.append(0)\n",
    "                                            kdp_con_centroid_lon.append(0)\n",
    "                                            kdp_con_centroid_lat.append(0)\n",
    "                                            kdp_con_storm_lon.append(0)\n",
    "                                            kdp_con_storm_lat.append(0)\n",
    "\n",
    "\n",
    "\n",
    "                                else:\n",
    "                                    #print(zdr_storm_lon[i[0]])\n",
    "                                    #Find the arc objects associated with this storm:\n",
    "                                    zdr_objects_lons = zdr_centroid_lons[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    zdr_objects_lats = zdr_centroid_lats[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    #print(\"zdr lons:\", zdr_objects_lons)\n",
    "                                    #print(\"arc lats:\", zdr_objects_lats)\n",
    "                                    #Get the sum of their areas\n",
    "                                    #print(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                    zdr_con_areas.append(np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    zdr_con_maxes.append(np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    zdr_con_means.append(np.mean(zdr_mean_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    zdr_con_max_lon.append(rlons_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                    zdr_con_max_lat.append(rlats_2d[np.where(ZDRmasked==np.max(zdr_max_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                    #print(\"Areas sum:\",zdr_con_areas)\n",
    "                                    #Find the actual centroids\n",
    "                                    weighted_lons = zdr_objects_lons * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    zdr_con_centroid_lon.append(np.sum(weighted_lons) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    weighted_lats = zdr_objects_lats * zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                    zdr_con_centroid_lat.append(np.sum(weighted_lats) / np.sum(zdr_areas_arr[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                    zdr_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                    zdr_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                                    zdr_con_masks.append(np.sum(zdr_masks[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])], axis=0, dtype=bool))\n",
    "                                    mask_con = np.sum(zdr_masks[np.where(zdr_storm_lon == zdr_storm_lon[i[0]])], axis=0, dtype=bool)\n",
    "                                    zdr_con_dev.append(np.std(ZDRmasked[mask_con]))\n",
    "                                    ZDRsorted = np.sort(ZDRmasked[mask_con])[::-1]\n",
    "                                    zdr_con_10max.append(np.mean(ZDRsorted[0:10]))\n",
    "                                    zdr_con_mode.append(stats.mode(ZDRmasked[mask_con]))\n",
    "                                    zdr_con_median.append(np.median(ZDRmasked[mask_con]))\n",
    "                                    #print(\"New centroid lon:\", zdr_con_centroid_lon, \"New centroid lat:\", zdr_con_centroid_lat)\n",
    "                                    #print(\"lons out of loop\", zdr_objects_lons)\n",
    "                                    try:\n",
    "                                        #Find the kdp objects associated with this storm:\n",
    "                                        kdp_objects_lons = kdp_centroid_lons[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        #print(\"kdp lons:\", kdp_objects_lons)\n",
    "                                        kdp_objects_lats = kdp_centroid_lats[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        #Get the sum of their areas\n",
    "                                        #print(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                        kdp_con_areas.append(np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        kdp_con_maxes.append(np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        kdp_con_max_lon.append(rlons_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                        kdp_con_max_lat.append(rlats_2d[np.where(KDPmasked==np.max(kdp_max_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))])\n",
    "                                        #Find the actual centroids\n",
    "                                        weighted_lons_kdp = kdp_objects_lons * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        kdp_con_centroid_lon.append(np.sum(weighted_lons_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        weighted_lats_kdp = kdp_objects_lats * kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]\n",
    "                                        #print(\"Could be it:\",\"weighted lons:\",weighted_lons_kdp, \"object lons\",kdp_objects_lons, \"areas:\",kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])])\n",
    "                                        kdp_con_centroid_lat.append(np.sum(weighted_lats_kdp) / np.sum(kdp_areas_arr[np.where(kdp_storm_lon == zdr_storm_lon[i[0]])]))\n",
    "                                        kdp_con_storm_lon.append(zdr_storm_lon[i[0]])\n",
    "                                        kdp_con_storm_lat.append(zdr_storm_lat[i[0]])\n",
    "                                    except:\n",
    "                                        print('storm missing kdp or zdr')\n",
    "                                        kdp_con_areas.append(0)\n",
    "                                        kdp_con_maxes.append(0)\n",
    "                                        kdp_con_max_lon.append(0)\n",
    "                                        kdp_con_max_lat.append(0)\n",
    "                                        kdp_con_centroid_lon.append(0)\n",
    "                                        kdp_con_centroid_lat.append(0)\n",
    "                                        kdp_con_storm_lon.append(0)\n",
    "                                        kdp_con_storm_lat.append(0)\n",
    "\n",
    "                            #Calculate KDP-ZDR separation\n",
    "                            print('calculating separation')\n",
    "                            kdp_con_centroid_lons1 = np.asarray(kdp_con_centroid_lon)\n",
    "                            kdp_con_centroid_lats1 = np.asarray(kdp_con_centroid_lat)\n",
    "                            zdr_con_centroid_lons1 = np.asarray(zdr_con_centroid_lon)\n",
    "                            zdr_con_centroid_lats1 = np.asarray(zdr_con_centroid_lat)\n",
    "                            #Eliminate consolidated arcs smaller than a specified area\n",
    "                            area = 2 #km*2\n",
    "                            zdr_con_areas_arr = np.asarray(zdr_con_areas)\n",
    "                            zdr_con_centroid_lats = zdr_con_centroid_lats1[zdr_con_areas_arr > area]\n",
    "                            zdr_con_centroid_lons = zdr_con_centroid_lons1[zdr_con_areas_arr > area]\n",
    "                            kdp_con_centroid_lats = kdp_con_centroid_lats1[zdr_con_areas_arr > area]\n",
    "                            kdp_con_centroid_lons = kdp_con_centroid_lons1[zdr_con_areas_arr > area]\n",
    "                            zdr_con_max_lons1 = np.asarray(zdr_con_max_lon)[zdr_con_areas_arr > area]\n",
    "                            zdr_con_max_lats1 = np.asarray(zdr_con_max_lat)[zdr_con_areas_arr > area]\n",
    "                            kdp_con_max_lons1 = np.asarray(kdp_con_max_lon)[zdr_con_areas_arr > area]\n",
    "                            kdp_con_max_lats1 = np.asarray(kdp_con_max_lat)[zdr_con_areas_arr > area]\n",
    "                            print('Boolean problem here')\n",
    "                            zdr_con_areas1 = zdr_con_areas_arr[zdr_con_areas_arr > area]\n",
    "\n",
    "                            kdp_inds = np.where(kdp_con_centroid_lats > 0)\n",
    "                            distance_kdp_zdr = g.inv(kdp_con_centroid_lons[kdp_inds], kdp_con_centroid_lats[kdp_inds], zdr_con_centroid_lons[kdp_inds], zdr_con_centroid_lats[kdp_inds])\n",
    "                            dist_kdp_zdr = distance_kdp_zdr[2] / 1000.\n",
    "                            #Now make an array for the distances which will have the same shape as the lats to prevent errors\n",
    "                            shaped_dist = np.zeros((np.shape(zdr_con_areas)))\n",
    "                            shaped_dist[kdp_inds] = dist_kdp_zdr\n",
    "                            print('maybe its here')\n",
    "                            #Do the same for the distances between the maxes\n",
    "                            distance_kdp_zdr_max = g.inv(kdp_con_max_lons1[kdp_inds], kdp_con_max_lats1[kdp_inds], zdr_con_max_lons1[kdp_inds], zdr_con_max_lats1[kdp_inds])\n",
    "                            dist_kdp_zdr_max = distance_kdp_zdr_max[2] / 1000.\n",
    "                            #Now make an array for the distances which will have the same shape as the lats to prevent errors\n",
    "                            shaped_dist_max = np.zeros((np.shape(zdr_con_areas)))\n",
    "                            shaped_dist_max[kdp_inds] = dist_kdp_zdr_max\n",
    "                            print('or not')\n",
    "                            print('made it to angles')\n",
    "                            #Get separation angle for KDP-ZDR centroids\n",
    "                            back_k = distance_kdp_zdr[1]\n",
    "                            #print('Raw back angle', back[i])\n",
    "                            for i in range(back_k.shape[0]):\n",
    "                                print('loop is ok')\n",
    "                                if distance_kdp_zdr[1][i] < 0:\n",
    "                                    print('if is ok')\n",
    "                                    back_k[i] = distance_kdp_zdr[1][i] + 360\n",
    "                            print('through loop 1')\n",
    "                            #print('fixed back', back[i])\n",
    "                            forw_k = np.abs(back_k - storm_relative_dir)\n",
    "                            #print('raw forw', forw[i])\n",
    "                            rawangle_k = back_k - storm_relative_dir\n",
    "                            #Account for weird angles\n",
    "                            for i in range(back_k.shape[0]):\n",
    "                                if forw_k[i] > 180:\n",
    "                                    #print('Big angle')\n",
    "                                    forw_k[i] = 360 - forw_k[i]\n",
    "                                    rawangle_k[i] = (360-forw_k[i])*(-1)\n",
    "                                #print(rawangle[i])\n",
    "                            print('through loop 2')\n",
    "                            rawangle_k = rawangle_k*(-1)\n",
    "                            \n",
    "                            #Now make an array for the distances which will have the same shape as the lats to prevent errors\n",
    "                            shaped_ang = np.zeros((np.shape(zdr_con_areas)))\n",
    "                            shaped_ang[kdp_inds] = rawangle_k\n",
    "                            shaped_ang = (180-np.abs(shaped_ang))*(shaped_ang/np.abs(shaped_ang))\n",
    "\n",
    "                        else:\n",
    "                            print('No ZDR arcs')\n",
    "                            kdp_areas = []\n",
    "                            kdp_centroid_lon = []\n",
    "                            kdp_centroid_lat = []\n",
    "                            kdp_storm_lon = []\n",
    "                            kdp_storm_lat = []\n",
    "                            zdr_con_centroid_lats = []\n",
    "                            zdr_con_centroid_lons = []\n",
    "                            kdp_con_centroid_lats = []\n",
    "                            kdp_con_centroid_lons = []\n",
    "                            kdp_con_area = []\n",
    "                            zdr_con_areas1 = []\n",
    "\n",
    "                        ###Now let's consolidate everything to fit the Pandas dataframe!\n",
    "                        p_zdr_areas = []\n",
    "                        p_zdr_maxes = []\n",
    "                        p_zdr_means = []\n",
    "                        p_zdr_devs = []\n",
    "                        p_zdr_10max = []\n",
    "                        p_zdr_mode = []\n",
    "                        p_zdr_median = []\n",
    "                        p_separations = []\n",
    "                        p_sp_angle = []\n",
    "                        if len(zdr_storm_lon) > 0:\n",
    "                            for storm in enumerate(max_lons_c):\n",
    "                                print(storm)\n",
    "                                print(np.flatnonzero(np.isclose(max_lons_c[storm[0]], zdr_con_storm_lon, rtol=1e-05)))\n",
    "                                matching_ind = np.flatnonzero(np.isclose(max_lons_c[storm[0]], zdr_con_storm_lon, rtol=1e-05))\n",
    "                                if matching_ind.shape[0] > 0:\n",
    "                                    p_zdr_areas.append((zdr_con_areas[matching_ind[0]]))\n",
    "                                    p_zdr_maxes.append((zdr_con_maxes[matching_ind[0]]))\n",
    "                                    p_zdr_means.append((zdr_con_means[matching_ind[0]]))\n",
    "                                    p_zdr_devs.append((zdr_con_dev[matching_ind[0]]))\n",
    "                                    p_zdr_10max.append((zdr_con_10max[matching_ind[0]]))\n",
    "                                    p_zdr_mode.append((zdr_con_mode[matching_ind[0]]))\n",
    "                                    p_zdr_median.append((zdr_con_median[matching_ind[0]]))\n",
    "                                    p_separations.append((shaped_dist[matching_ind[0]]))\n",
    "                                    p_sp_angle.append((shaped_ang[matching_ind[0]]))\n",
    "                                else:\n",
    "                                    p_zdr_areas.append((0))\n",
    "                                    p_zdr_maxes.append((0))\n",
    "                                    p_zdr_means.append((0))\n",
    "                                    p_zdr_devs.append((0))\n",
    "                                    p_zdr_10max.append((0))\n",
    "                                    p_zdr_mode.append((0))\n",
    "                                    p_zdr_median.append((0))\n",
    "                                    p_separations.append((0))\n",
    "                                    p_sp_angle.append((0))\n",
    "                        else:\n",
    "                            for storm in enumerate(max_lons_c):\n",
    "                                p_zdr_areas.append((0))\n",
    "                                p_zdr_maxes.append((0))\n",
    "                                p_zdr_means.append((0))\n",
    "                                p_zdr_devs.append((0))\n",
    "                                p_zdr_10max.append((0))\n",
    "                                p_zdr_mode.append((0))\n",
    "                                p_zdr_median.append((0))\n",
    "                                p_separations.append((0))\n",
    "                                p_sp_angle.append((0))\n",
    "\n",
    "                        #Now start plotting stuff!\n",
    "                        print('made it through giant if statement')\n",
    "                        if np.asarray(zdr_centroid_lon).shape[0] > 0:\n",
    "                            ax.scatter(zdr_centroid_lon, zdr_centroid_lat, marker = '*', s = 100, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                            #ax.scatter(zdr_con_max_lon, zdr_con_max_lat, marker = '*', s = 100, color = 'purple', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                        if np.asarray(kdp_centroid_lon).shape[0] > 0:\n",
    "                            ax.scatter(kdp_centroid_lon, kdp_centroid_lat, marker = '^', s = 100, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                            #ax.scatter(kdp_con_max_lon, kdp_con_max_lat, marker = '^', s = 100, color = 'purple', zorder = 10, transform=ccrs.PlateCarree())\n",
    "                        print(\"plotted centroids\")\n",
    "                        #Uncomment to print all object areas\n",
    "                        #for i in enumerate(zdr_areas):\n",
    "                        #    plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2\" %(zdr_areas[i[0]].magnitude), size = 23)\n",
    "                            #plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f km^2 / %.2f km / %.2f dB\" %(zdr_areas[i[0]].magnitude, zdr_dist[i[0]], zdr_forw[i[0]]), size = 23)\n",
    "                            #plt.annotate(zdr_areas[i[0]], (zdr_centroid_lon[i[0]],zdr_centroid_lat[i[0]]))\n",
    "                        #ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels1,linewide = .01, colors ='b', alpha = .5)\n",
    "                        #plt.tight_layout()\n",
    "                        #plt.savefig('ZDRarcannotated.png')\n",
    "                        storm_times = []\n",
    "                        for l in range(len(max_lons_c)):\n",
    "                            storm_times.append((time_start))\n",
    "\n",
    "                    #If there are no storms, set everything to empty arrays!\n",
    "                    else:\n",
    "                        print('Filling arrays with no storms')\n",
    "                        storm_ids = []\n",
    "                        storm_ids = []\n",
    "                        max_lons_c = []\n",
    "                        max_lats_c = []\n",
    "                        p_zdr_areas = []\n",
    "                        p_zdr_maxes = []\n",
    "                        p_zdr_means = []\n",
    "                        p_zdr_devs = []\n",
    "                        p_zdr_10max = []\n",
    "                        p_zdr_mode = []\n",
    "                        p_zdr_median = []\n",
    "                        p_separations = []\n",
    "                        p_sp_angle = []\n",
    "                        zdr_con_areas1 = []\n",
    "                        storm_times = time_start\n",
    "                    #Now record all data in a Pandas dataframe.\n",
    "                    print('making dataframe')\n",
    "                    new_cells = pd.DataFrame({\n",
    "                        'scan': scan_index,\n",
    "                        'storm_id' : storm_ids,\n",
    "                        'storm_id1' : storm_ids,\n",
    "                        'storm_lon' : max_lons_c,\n",
    "                        'storm_lat' : max_lats_c,\n",
    "                        'zdr_area' : p_zdr_areas,\n",
    "                        'zdr_max' : p_zdr_maxes,\n",
    "                        'zdr_mean' : p_zdr_means,\n",
    "                        'zdr_std' : p_zdr_devs,\n",
    "                        'zdr_10max' : p_zdr_10max,\n",
    "                        'zdr_mode' : p_zdr_mode,\n",
    "                        'zdr_median' : p_zdr_median,\n",
    "                        'kdp_zdr_sep' : p_separations,\n",
    "                        'kdp_zdr_angle' : p_sp_angle,\n",
    "                        'times' : storm_times\n",
    "                    })\n",
    "                    print('setting index')\n",
    "                    new_cells.set_index(['scan', 'storm_id'], inplace=True)\n",
    "                    if scan_index == 0:\n",
    "                        print('first dataframe')\n",
    "                        tracks_dataframe = new_cells\n",
    "                    else:\n",
    "                        tracks_dataframe = tracks_dataframe.append(new_cells)\n",
    "                    n = n+1\n",
    "                    scan_index = scan_index + 1\n",
    "                    #max_lons_p = max_lons_c\n",
    "                    #max_lats_p = max_lats_c\n",
    "                    #storm_ids_p = storm_ids\n",
    "                    #Plot the consolidated stuff!\n",
    "                    #Write some text objects for the ZDR arc attributes to add to the placefile\n",
    "                    f.write(\"Color: 139 000 000 \\n\")\n",
    "                    f.write('Font: 1, 30, 1,\"Arial\" \\n')\n",
    "                    print('wrote font')\n",
    "                    for y in range(len(p_zdr_areas)):\n",
    "                        print(y)\n",
    "                        #f.write('Text: '+str(max_lats_c[y])+','+str(max_lons_c[y])+', 1, \"X\",\" Arc Area: '+str(p_zdr_areas[y])+'\\\\n Arc Mean: '+str(p_zdr_means[y])+'\\\\n KDP-ZDR Separation: '+str(p_separations[y])+'\\\\n Separation Angle: '+str(p_sp_angle[y])+'\" \\n')\n",
    "                        f.write('Text: '+str(max_lats_c[y])+','+str(max_lons_c[y])+', 1, \"X\",\" Arc Area: %.2f km^2 \\\\n Arc Mean: %.2f dB \\\\n Arc 10 Max Mean: %.2f dB \\\\n KDP-ZDR Separation: %.2f km \\\\n Separation Angle: %.2f degrees\" \\n' %(p_zdr_areas[y], p_zdr_means[y], p_zdr_10max[y], p_separations[y], p_sp_angle[y]))\n",
    "                    print('made it to plotting')\n",
    "                    if ((len(zdr_con_areas1) > 0) & (len(max_lons_c) > 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",
    "                        else:\n",
    "                            print('No kdp')\n",
    "                    else:\n",
    "                        print('No zdr arcs')\n",
    "                    print(\"means there's a kdp problem\")\n",
    "                    #hour = time_start.hour\n",
    "                    #if hour < 10:\n",
    "                    #    hour = '0'+str(hour)\n",
    "                    #minute = time_start.minute\n",
    "                    #if minute < 10:\n",
    "                    #    minute = '0'+str(minute)\n",
    "                    #day = time_start.day\n",
    "                    #if day < 10:\n",
    "                    #    day = '0'+str(day)\n",
    "                    title_plot = plt.title(station+' Radar Reflectivity, ZDR, and KDP '+str(time_start.year)+'-'+str(time_start.month)+'-'+str(time_start.day)+\n",
    "                                               ' '+str(hour)+':'+str(minute)+' UTC', size = 25)\n",
    "                    #if np.asarray(zdr_storm_lon).shape[0] > 0:\n",
    "                    #    ax.scatter(zdr_storm_lon,zdr_storm_lat, marker = \"o\", color = 'purple', s = 500)\n",
    "                    #if np.asarray(kdp_storm_lon).shape[0] > 0:\n",
    "                    #    ax.scatter(kdp_storm_lon,kdp_storm_lat, marker = \"o\", color = 'purple', s = 500)\n",
    "                    #try:\n",
    "                    #    ax.scatter(max_lons_c,max_lats_c, marker = \"o\", color = 'k', s = 500, alpha = .6)\n",
    "                    #except:\n",
    "                    #    \"No storm centroids found\"\n",
    "                    try:\n",
    "                        plt.plot([zdr_con_centroid_lons[kdp_inds], kdp_con_centroid_lons[kdp_inds]], [zdr_con_centroid_lats[kdp_inds],kdp_con_centroid_lats[kdp_inds]], color = 'k', linewidth = 5, transform=ccrs.PlateCarree())\n",
    "                    except:\n",
    "                        print('KDP-ZDR separation didt work')\n",
    "                    ref_centroid_lon = max_lons_c\n",
    "                    ref_centroid_lat = max_lats_c\n",
    "                    if len(max_lons_c) > 0:\n",
    "                        ax.scatter(max_lons_c,max_lats_c, marker = \"o\", color = 'k', s = 500, alpha = .6)\n",
    "                        for i in enumerate(ref_centroid_lon): \n",
    "                            plt.text(ref_centroid_lon[i[0]]+.016, ref_centroid_lat[i[0]]+.016, \"storm_id: %.1f\" %(storm_ids[i[0]]), size = 25)\n",
    "                    #Comment out this line if not plotting tornado tracks\n",
    "                    #plt.plot([start_torlons, end_torlons], [start_torlats, end_torlats], color = 'purple', linewidth = 5, transform=ccrs.PlateCarree())\n",
    "                    #Add legend stuff\n",
    "                    zdr_outline = mlines.Line2D([], [], color='blue', linewidth = 5, linestyle = 'solid', label='ZDR Arc Outline(Area/Max)')\n",
    "                    kdp_outline = mlines.Line2D([], [], color='green', linewidth = 5,linestyle = 'solid', label='\"KDP Foot\" Outline')\n",
    "                    separation_vector = mlines.Line2D([], [], color='black', linewidth = 5,linestyle = 'solid', label='KDP/ZDR Centroid Separation Vector (Red Text=Distance)')\n",
    "                    #tor_track = mlines.Line2D([], [], color='purple', linewidth = 5,linestyle = 'solid', label='Tornado Tracks')\n",
    "                    elevation = mlines.Line2D([], [], color='grey', linewidth = 5,linestyle = 'solid', label='Height AGL (m)')\n",
    "\n",
    "                    plt.legend(handles=[zdr_outline, kdp_outline, separation_vector, elevation], loc = 3, fontsize = 25)\n",
    "                    alt_levs = [1000, 2000]\n",
    "                    cele = ax.contour(ungrid_lons,ungrid_lats,gate_altitude-radar.altitude['data'][0],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",
    "                    print('made it to saving')\n",
    "                    plt.savefig('Machine_Learning/ZSlideR325P_ZDRArc_example'+station+str(time_start.year)+str(time_start.month)+str(day)+str(hour)+str(minute)+'.png')\n",
    "                    print('figure saved')\n",
    "                    plt.close()\n",
    "\n",
    "        except:\n",
    "            continue\n",
    "    f.close()\n",
    "    plt.show()\n",
    "    return tracks_dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:55:15.826813\n",
      "Downloaded KPOE20170402_173137_V06\n",
      "1 out of 1 files downloaded...0 errors\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\sklearn\\base.py:306: UserWarning: Trying to unpickle estimator DecisionTreeClassifier from version 0.19.0 when using version 0.21.2. This might lead to breaking code or invalid results. Use at your own risk.\n",
      "  UserWarning)\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\sklearn\\base.py:306: UserWarning: Trying to unpickle estimator RandomForestClassifier from version 0.19.0 when using version 0.21.2. This might lead to breaking code or invalid results. Use at your own risk.\n",
      "  UserWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "in loop\n",
      "23\n",
      "0.41806\n",
      "2017-04-02 17:31:37.003000\n",
      "its this line\n",
      "heres the problem\n",
      "2019-08-01 19:55:37.023559\n",
      "2019-08-01 19:55:37.077529\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:4185: UserWarning: Warning: converting a masked element to nan.\n",
      "  warnings.warn(\"Warning: converting a masked element to nan.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:56:19.377360\n",
      "almost gridding\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\pyart\\map\\gates_to_grid.py:177: DeprecationWarning: Barnes weighting function is deprecated. Please use Barnes 2 to be consistent with Pauley and Wu 1990.\n",
      "  \" Pauley and Wu 1990.\", DeprecationWarning)\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:852: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  return umath.absolute(a) * self.tolerance >= umath.absolute(b)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:56:20.292484\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "[[-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] ...,  [-- -- -- ..., 1.635841965675354 1.7074615955352783 --] [-- -- -- ..., 1.6195695400238037 -- --] [-- -- -- ..., 1.607196569442749 -- --]] radian\n",
      "3.141324758529663 radian\n",
      "3.218650817871094e-05 radian\n",
      "[[  0.           0.           0.         ...,   0.           0.           0.        ] [  0.           0.           0.         ...,   0.           0.           0.        ] [  0.           0.           0.         ...,   0.           0.           0.        ] ...,  [  0.           0.           0.         ...,  93.72684479  97.83034515    0.        ] [  0.           0.           0.         ...,  92.79450226   0.           0.        ] [  0.           0.           0.         ...,  92.08557892   0.           0.        ]] degree\n",
      "179.98464965820312 degree\n",
      "0.0 degree\n",
      "Set up our projection\n",
      "(-94.476112365722656, 29.655277252197266, <cartopy.crs.PlateCarree object at 0x0000027A4F68E360>)\n",
      "plotting\n",
      "1\n",
      "-93.8745061742109\n",
      "past mask\n",
      "45.3799\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-92.76402999250013\n",
      "past mask\n",
      "44.4531\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-93.14733677509034\n",
      "past mask\n",
      "46.1871\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-92.80884957297816\n",
      "past mask\n",
      "45.1691\n",
      "found a storm\n",
      "added polygon\n",
      "1\n",
      "-92.20884305336274\n",
      "past mask\n",
      "43.5488\n",
      "1\n",
      "-94.14059021875602\n",
      "past mask\n",
      "45.9471\n",
      "found a storm\n",
      "added polygon\n",
      "tracking id\n",
      "[]\n",
      "made it here\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "made it through zdr centroids\n",
      "1\n",
      "-93.15218283145637\n",
      "[ 179.10538729   37.42479042    0.66738752   34.96627706  174.04462136]\n",
      "1\n",
      "-93.18248955085376\n",
      "1\n",
      "-92.83322305160556\n",
      "1\n",
      "-92.84916036210242\n",
      "[ 182.09246792   18.20958293   30.91082037    3.89537309  182.37774105]\n",
      "1\n",
      "-92.78868254023875\n",
      "[ 182.25692872   19.99120125   37.62775589    3.53639361  183.70113182]\n",
      "1\n",
      "-94.18453567144114\n",
      "[  38.65634868  198.03392677  171.7078876   182.89182255    7.6391842 ]\n",
      "made it through kdp centroids\n",
      "0\n",
      "storm missing kdp or zdr\n",
      "1\n",
      "-93.1473367751\n",
      "[ 22.68821826]\n",
      "[ 3.66465527]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it to angles\n",
      "loop is ok\n",
      "if is ok\n",
      "through loop 1\n",
      "through loop 2\n",
      "(0, -93.874506174210893)\n",
      "[]\n",
      "(1, -92.764029992500127)\n",
      "[0]\n",
      "(2, -93.147336775090338)\n",
      "[1]\n",
      "(3, -92.808849572978161)\n",
      "[]\n",
      "(4, -94.140590218756017)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "first dataframe\n",
      "wrote font\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "made it to plotting\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\ipykernel_launcher.py:931: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "made it to saving\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "0.376236\n",
      "2017-04-02 17:33:33.729000\n",
      "its this line\n",
      "heres the problem\n",
      "2019-08-01 19:56:32.503830\n",
      "2019-08-01 19:56:32.547770\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:4185: UserWarning: Warning: converting a masked element to nan.\n",
      "  warnings.warn(\"Warning: converting a masked element to nan.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:57:07.356881\n",
      "almost gridding\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\pyart\\map\\gates_to_grid.py:177: DeprecationWarning: Barnes weighting function is deprecated. Please use Barnes 2 to be consistent with Pauley and Wu 1990.\n",
      "  \" Pauley and Wu 1990.\", DeprecationWarning)\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:852: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  return umath.absolute(a) * self.tolerance >= umath.absolute(b)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:57:08.280354\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "[[-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] ...,  [-- -- -- ..., 1.514047622680664 1.6215323209762573 1.7420603036880493] [-- -- -- ..., 1.524146318435669 1.586686134338379 1.657728910446167] [-- -- -- ..., 1.5338221788406372 1.574331283569336 1.6204575300216675]] radian\n",
      "3.141566514968872 radian\n",
      "5.53131103515625e-05 radian\n",
      "[[  0.           0.           0.         ...,   0.           0.           0.        ] [  0.           0.           0.         ...,   0.           0.           0.        ] [  0.           0.           0.         ...,   0.           0.           0.        ] ...,  [  0.           0.           0.         ...,  86.74854279  92.90695953   99.81270599] [  0.           0.           0.         ...,  87.32715607  90.91042328   94.98087311] [  0.           0.           0.         ...,  87.88153839  90.20253754   92.84537506]] degree\n",
      "179.99850463867188 degree\n",
      "0.0 degree\n",
      "Set up our projection\n",
      "(-94.476112365722656, 29.655277252197266, <cartopy.crs.PlateCarree object at 0x0000027A4F68E410>)\n",
      "plotting\n",
      "1\n",
      "-94.26333108322784\n",
      "past mask\n",
      "43.0158\n",
      "1\n",
      "-93.87220796826519\n",
      "past mask\n",
      "45.1118\n",
      "found a storm\n",
      "[   0.86205848  200.82118087  180.41242592  184.2380096    32.06987413]\n",
      "Poly lon -93.87220796826519\n",
      "[-93.87450617 -92.76402999 -93.14733678 -92.80884957 -94.14059022]\n",
      "[0 1 2 3 4]\n",
      "storm id 0\n",
      "added polygon\n",
      "1\n",
      "-92.79434118270257\n",
      "past mask\n",
      "43.5764\n",
      "1\n",
      "-92.75084859584598\n",
      "past mask\n",
      "44.4781\n",
      "found a storm\n",
      "[ 198.43167199    3.15875002   37.93656125   15.04741003  198.63255101]\n",
      "Poly lon -92.75084859584598\n",
      "[-93.87450617 -92.76402999 -93.14733678 -92.80884957 -94.14059022]\n",
      "[0 1 2 3 4]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-93.138589232426\n",
      "past mask\n",
      "46.7769\n",
      "found a storm\n",
      "[ 177.4317758    36.51867538    2.95637155   32.97628663  172.76710859]\n",
      "Poly lon -93.138589232426\n",
      "[-93.87450617 -92.76402999 -93.14733678 -92.80884957 -94.14059022]\n",
      "[0 1 2 3 4]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-92.79770197655365\n",
      "past mask\n",
      "44.8597\n",
      "found a storm\n",
      "[ 182.41823876   19.34782376   36.52619217    2.4373326   183.66433764]\n",
      "Poly lon -92.79770197655365\n",
      "[-93.87450617 -92.76402999 -93.14733678 -92.80884957 -94.14059022]\n",
      "[0 1 2 3 4]\n",
      "storm id 3\n",
      "added polygon\n",
      "1\n",
      "-94.14960380611369\n",
      "past mask\n",
      "46.0018\n",
      "found a storm\n",
      "[  31.46349022  201.34883032  176.04280616  185.88039734    1.35174342]\n",
      "Poly lon -94.14960380611369\n",
      "[-93.87450617 -92.76402999 -93.14733678 -92.80884957 -94.14059022]\n",
      "[0 1 2 3 4]\n",
      "storm id 4\n",
      "added polygon\n",
      "tracking id\n",
      "[]\n",
      "made it here\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "made it through zdr centroids\n",
      "1\n",
      "-93.15274922880074\n",
      "[ 177.28706882   38.69364952    1.3932376    35.84846668  173.0981922 ]\n",
      "1\n",
      "-92.91215522706175\n",
      "1\n",
      "-92.82205128005559\n",
      "1\n",
      "-92.77266677786731\n",
      "[ 181.28355509   19.97209262   38.28184599    4.38874092  184.03710201]\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "-92.7977019766\n",
      "[ 14.54251109]\n",
      "[ 3.36312261]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it to angles\n",
      "loop is ok\n",
      "if is ok\n",
      "loop is ok\n",
      "through loop 1\n",
      "through loop 2\n",
      "(0, -93.872207968265187)\n",
      "[]\n",
      "(1, -92.75084859584598)\n",
      "[]\n",
      "(2, -93.138589232426)\n",
      "[0]\n",
      "(3, -92.797701976553654)\n",
      "[1]\n",
      "(4, -94.149603806113689)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "wrote font\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "made it to plotting\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\ipykernel_launcher.py:902: 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": [
      "made it to saving\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "0.376942\n",
      "2017-04-02 17:34:58.888000\n",
      "its this line\n",
      "heres the problem\n",
      "2019-08-01 19:57:21.036479\n",
      "2019-08-01 19:57:21.088448\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:4185: UserWarning: Warning: converting a masked element to nan.\n",
      "  warnings.warn(\"Warning: converting a masked element to nan.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:57:59.088729\n",
      "almost gridding\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\pyart\\map\\gates_to_grid.py:177: DeprecationWarning: Barnes weighting function is deprecated. Please use Barnes 2 to be consistent with Pauley and Wu 1990.\n",
      "  \" Pauley and Wu 1990.\", DeprecationWarning)\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:852: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  return umath.absolute(a) * self.tolerance >= umath.absolute(b)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:58:00.038184\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "[[-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] ...,  [-- -- -- ..., 1.2244861125946045 1.4257316589355469 1.6307663917541504] [-- -- -- ..., 1.5008009672164917 1.67899489402771 1.8651154041290283] [-- -- -- ..., 1.5640323162078857 1.6976892948150635 1.8421666622161865]] radian\n",
      "3.1414265632629395 radian\n",
      "1.9073486328125e-06 radian\n",
      "[[   0.            0.            0.         ...,    0.            0.            0.        ] [   0.            0.            0.         ...,    0.            0.            0.        ] [   0.            0.            0.         ...,    0.            0.            0.        ] ...,  [   0.            0.            0.         ...,   70.15789032   81.6884079    93.43603516] [   0.            0.            0.         ...,   85.98956299    96.19932556  106.8632431 ] [   0.            0.            0.         ...,   89.61244965    97.27043152  105.54837799]] degree\n",
      "179.990478515625 degree\n",
      "0.0 degree\n",
      "Set up our projection\n",
      "(-94.476112365722656, 29.655277252197266, <cartopy.crs.PlateCarree object at 0x0000027A4F68E2B0>)\n",
      "plotting\n",
      "1\n",
      "-94.25005905713937\n",
      "past mask\n",
      "43.7952\n",
      "found a storm\n",
      "[ 148.77411151  145.16173674  107.92825364  138.82601962  126.07647466]\n",
      "Poly lon -94.25005905713937\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-93.86354682673245\n",
      "past mask\n",
      "44.7165\n",
      "found a storm\n",
      "[   0.82326031  198.6779047   177.90129147  182.62818949   32.8527638 ]\n",
      "Poly lon -93.86354682673245\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 0\n",
      "added polygon\n",
      "1\n",
      "-92.78232973925768\n",
      "past mask\n",
      "44.4124\n",
      "found a storm\n",
      "[ 204.50290158    8.75111768   36.27957605   24.44291745  204.31763311]\n",
      "Poly lon -92.78232973925768\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-92.74426863338388\n",
      "past mask\n",
      "44.1951\n",
      "found a storm\n",
      "[ 197.6714157     2.08521408   37.76213476   15.08415863  198.91049262]\n",
      "Poly lon -92.74426863338388\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-93.13321896736787\n",
      "past mask\n",
      "46.8808\n",
      "found a storm\n",
      "[ 176.76546786   37.03002393    1.75797027   33.65072507  172.99227175]\n",
      "Poly lon -93.13321896736787\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-92.79318143651271\n",
      "past mask\n",
      "44.6454\n",
      "found a storm\n",
      "[ 182.41997194   17.62136497   35.42997828    1.06060731  184.6078289 ]\n",
      "Poly lon -92.79318143651271\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 3\n",
      "added polygon\n",
      "1\n",
      "-94.04683879255242\n",
      "past mask\n",
      "43.83\n",
      "found a storm\n",
      "[  35.36293158  184.12515984  158.24423242  169.21981338   15.9084651 ]\n",
      "Poly lon -94.04683879255242\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 4\n",
      "added polygon\n",
      "1\n",
      "-94.14977034870486\n",
      "past mask\n",
      "45.4633\n",
      "found a storm\n",
      "[  31.07173202  201.39794798  175.69276644  186.37607555    1.88492854]\n",
      "Poly lon -94.14977034870486\n",
      "[-93.87220797 -92.7508486  -93.13858923 -92.79770198 -94.14960381]\n",
      "[0 1 2 3 4]\n",
      "storm id 4\n",
      "added polygon\n",
      "tracking id\n",
      "[6]\n",
      "made it here\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "made it through zdr centroids\n",
      "1\n",
      "-92.79680962125609\n",
      "1\n",
      "-93.14589344987668\n",
      "[ 106.79624362  175.28979809   37.8643819    38.68659768    1.4836312\n",
      "   35.24869723  155.73463685  173.18366754]\n",
      "1\n",
      "-92.80446598895757\n",
      "1\n",
      "-92.74916234270121\n",
      "1\n",
      "-92.7852474271316\n",
      "made it through kdp centroids\n",
      "0\n",
      "storm missing kdp or zdr\n",
      "1\n",
      "-93.1332189674\n",
      "[ 35.12899241]\n",
      "[ 15.17303117]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it to angles\n",
      "loop is ok\n",
      "if is ok\n",
      "through loop 1\n",
      "through loop 2\n",
      "(0, -94.250059057139367)\n",
      "[]\n",
      "(1, -93.863546826732446)\n",
      "[]\n",
      "(2, -92.782329739257676)\n",
      "[0]\n",
      "(3, -92.744268633383882)\n",
      "[]\n",
      "(4, -93.133218967367867)\n",
      "[1]\n",
      "(5, -92.793181436512711)\n",
      "[]\n",
      "(6, -94.046838792552421)\n",
      "[]\n",
      "(7, -94.149770348704862)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "wrote font\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "made it to plotting\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\ipykernel_launcher.py:931: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "made it to saving\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "0.421162\n",
      "2017-04-02 17:36:36.927000\n",
      "its this line\n",
      "heres the problem\n",
      "2019-08-01 19:58:12.235216\n",
      "2019-08-01 19:58:12.281191\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:4185: UserWarning: Warning: converting a masked element to nan.\n",
      "  warnings.warn(\"Warning: converting a masked element to nan.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:58:49.774767\n",
      "almost gridding\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\pyart\\map\\gates_to_grid.py:177: DeprecationWarning: Barnes weighting function is deprecated. Please use Barnes 2 to be consistent with Pauley and Wu 1990.\n",
      "  \" Pauley and Wu 1990.\", DeprecationWarning)\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\numpy\\ma\\core.py:852: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  return umath.absolute(a) * self.tolerance >= umath.absolute(b)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-08-01 19:58:50.664259\n",
      "made it to smoothing\n",
      "made it through gradient\n",
      "[[-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] [-- -- -- ..., -- -- --] ...,  [-- -- -- ..., 0.1866157054901123 0.12651562690734863 0.09039306640625] [-- -- -- ..., 0.2915985584259033 0.20787476003170013 0.15343499183654785] [-- -- -- ..., 0.41399598121643066 0.30831098556518555 0.23376131057739258]] radian\n",
      "3.1413204669952393 radian\n",
      "5.4836273193359375e-06 radian\n",
      "[[  0.           0.           0.         ...,   0.           0.           0.        ] [  0.           0.           0.         ...,   0.           0.           0.        ] [  0.           0.           0.         ...,   0.           0.           0.        ] ...,  [  0.           0.           0.         ...,  10.69229221   7.24881172    5.17914104] [  0.           0.           0.         ...,  16.70736694  11.91034698    8.79117775] [  0.           0.           0.         ...,  23.72022247  17.6649189   13.39353657]] degree\n",
      "179.98440551757812 degree\n",
      "0.0 degree\n",
      "Set up our projection\n",
      "(-94.476112365722656, 29.655277252197266, <cartopy.crs.PlateCarree object at 0x0000027A42C86BA0>)\n",
      "plotting\n",
      "1\n",
      "-94.2489021345034\n",
      "past mask\n",
      "43.5084\n",
      "found a storm\n",
      "[   2.56791302  146.57291771  144.32167012  145.73800366  108.46628967\n",
      "  139.22015288  112.18642803  125.38813054]\n",
      "Poly lon -94.2489021345034\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 5\n",
      "added polygon\n",
      "1\n",
      "-94.25903726309772\n",
      "past mask\n",
      "42.851\n",
      "1\n",
      "-93.86010566945454\n",
      "past mask\n",
      "44.4785\n",
      "found a storm\n",
      "[ 149.37113178    0.38612731  204.20687306  197.32773737  176.61443538\n",
      "  182.05050427   36.19332603   32.20879166]\n",
      "Poly lon -93.86010566945454\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 0\n",
      "added polygon\n",
      "1\n",
      "-92.76908487983212\n",
      "past mask\n",
      "44.6184\n",
      "found a storm\n",
      "[ 144.58982793  202.87204145    2.59234759    8.29025039   36.79036441\n",
      "   23.2205723   187.49080192  204.82357921]\n",
      "Poly lon -92.76908487983212\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-92.73096094243144\n",
      "past mask\n",
      "44.2043\n",
      "found a storm\n",
      "[ 146.29086977  195.56797819   13.99779364    3.16919488   38.45914354\n",
      "   13.62486401  181.79235348  198.99940321]\n",
      "Poly lon -92.73096094243144\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 1\n",
      "added polygon\n",
      "1\n",
      "-93.123072080905\n",
      "past mask\n",
      "46.5797\n",
      "found a storm\n",
      "[ 108.73753944  174.89594407   36.45729789   36.72357831    2.41663226\n",
      "   32.77088467  155.80958785  173.25744631]\n",
      "Poly lon -93.123072080905\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 2\n",
      "added polygon\n",
      "1\n",
      "-92.78599203859187\n",
      "past mask\n",
      "44.6552\n",
      "found a storm\n",
      "[ 139.85169243  181.88415691   26.03690886   16.32398887   35.23329362\n",
      "    0.95706835  168.83630239  185.95671273]\n",
      "Poly lon -92.78599203859187\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 3\n",
      "added polygon\n",
      "1\n",
      "-92.77115905085584\n",
      "past mask\n",
      "43.0276\n",
      "1\n",
      "-94.03907355282945\n",
      "past mask\n",
      "45.1467\n",
      "found a storm\n",
      "[ 117.38497066   33.25744955  189.92995848  184.46069581  158.85191682\n",
      "  170.13311936    2.66645419   15.8564616 ]\n",
      "Poly lon -94.03907355282945\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 6\n",
      "added polygon\n",
      "1\n",
      "-94.14886975843238\n",
      "past mask\n",
      "45.2848\n",
      "found a storm\n",
      "[ 129.28950275   31.06007019  206.75480922  201.24810926  175.62519274\n",
      "  186.86616832   18.52625055    1.33479777]\n",
      "Poly lon -94.14886975843238\n",
      "[-94.25005906 -93.86354683 -92.78232974 -92.74426863 -93.13321897\n",
      " -92.79318144 -94.04683879 -94.14977035]\n",
      "[5 0 1 1 2 3 6 4]\n",
      "storm id 4\n",
      "added polygon\n",
      "tracking id\n",
      "[6]\n",
      "made it here\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 1.]\n",
      "arc\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "got xc\n",
      "got array1\n",
      "got array2\n",
      "got array3\n",
      "got array4\n",
      "got array5\n",
      "got array6\n",
      "got array7\n",
      "got array8\n",
      "got array9\n",
      "got array10\n",
      "got array11\n",
      "got array12\n",
      "got to prediction\n",
      "[ 0.]\n",
      "made it through zdr centroids\n",
      "1\n",
      "-92.76995898243875\n",
      "[ 145.09131374  203.4298348     0.68481851   12.11396344   36.85390861\n",
      "   24.50076626  189.44308134  206.26346177]\n",
      "1\n",
      "-93.13645869636825\n",
      "[ 107.5995297   173.61152924   38.19689893   38.97235036    1.57328861\n",
      "   34.59819802  156.05386042  172.83846857]\n",
      "1\n",
      "-93.10731276830617\n",
      "1\n",
      "-92.82347053060471\n",
      "1\n",
      "-92.74539369272487\n",
      "made it through kdp centroids\n",
      "0\n",
      "1\n",
      "-93.1230720809\n",
      "[ 38.641176]\n",
      "[ 13.24486993]\n",
      "calculating separation\n",
      "Boolean problem here\n",
      "maybe its here\n",
      "or not\n",
      "made it to angles\n",
      "loop is ok\n",
      "if is ok\n",
      "loop is ok\n",
      "if is ok\n",
      "through loop 1\n",
      "through loop 2\n",
      "(0, -94.248902134503396)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devnew\\lib\\site-packages\\ipykernel_launcher.py:902: 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": [
      "[]\n",
      "(1, -93.860105669454541)\n",
      "[]\n",
      "(2, -92.769084879832121)\n",
      "[0]\n",
      "(3, -92.730960942431437)\n",
      "[]\n",
      "(4, -93.123072080905004)\n",
      "[1]\n",
      "(5, -92.785992038591871)\n",
      "[]\n",
      "(6, -94.03907355282945)\n",
      "[]\n",
      "(7, -94.148869758432383)\n",
      "[]\n",
      "made it through giant if statement\n",
      "plotted centroids\n",
      "making dataframe\n",
      "setting index\n",
      "wrote font\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "made it to plotting\n",
      "consolidated ZDR:\n",
      "consolidated ZDR:\n",
      "means there's a kdp problem\n",
      "made it to saving\n",
      "figure saved\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "in loop\n",
      "23\n",
      "2019-08-01 19:59:01.816580\n"
     ]
    }
   ],
   "source": [
    "#Loop to run the actual algorithm\n",
    "print(datetime.utcnow())\n",
    "#You can either run all example cases (commented out line) or one at a time (runing line)\n",
    "#for i in range(len(durationstm)):\n",
    "for i in [0]:\n",
    "    tracks_dataframe = pd.DataFrame()\n",
    "    tracks_dataframe = multi_case_algorithm_ML1(storm_relative_dirstm[i], zdrlevstm[i], kdplevstm[i], REFlevstm[i], REFlev1stm[i], big_stormstm[i], zero_z_triggerstm[i], storm_to_trackstm[i], yearstm[i], monthstm[i], daystm[i], hourstm[i], start_minstm[i], durationstm[i], stationstm[i])\n",
    "    #Save the dataframe off as a pickle file\n",
    "    tracks_dataframe.to_pickle('MLexample_valid_tor'+str(yearstm[i])+str(monthstm[i])+str(daystm[i])+str(stationstm[i])+'.pkl')\n",
    "print(datetime.utcnow())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now have three different forms of algorithm output:\n",
    "    1. Saved .png images of each radar scan.\n",
    "    2. A looping GR2 placefile of all algorithm output.\n",
    "    3. A Pickle file containing tracks for all storms and their associated ZDR arc characteristics."
   ]
  },
  {
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}