ZDR Arc Detector Tropical.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### This script will plot some graphics showing the relative locations of the ZDR and KDP maxes in the storm.\n",
    "import matplotlib.pyplot as plt\n",
    "import pyart\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "from metpy.units import atleast_1d, check_units, concatenate, units\n",
    "from matplotlib.patches import PathPatch\n",
    "from matplotlib.path import Path\n",
    "from siphon.radarserver import RadarServer\n",
    "#rs = RadarServer('http://thredds-aws.unidata.ucar.edu/thredds/radarServer/nexrad/level2/S3/')\n",
    "rs = RadarServer('http://thredds.ucar.edu/thredds/radarServer/nexrad/level2/IDD/')\n",
    "from datetime import datetime, timedelta\n",
    "from siphon.cdmr import Dataset\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from cartopy.io.shapereader import Reader\n",
    "from cartopy.feature import ShapelyFeature\n",
    "from metpy.units import atleast_1d, check_units, concatenate, units\n",
    "from shapely.geometry import polygon as sp\n",
    "import pyproj \n",
    "import shapely.ops as ops\n",
    "from shapely.ops import transform\n",
    "from shapely.geometry.polygon import Polygon\n",
    "from functools import partial\n",
    "from shapely import geometry\n",
    "import netCDF4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matts\\Miniconda3\\envs\\devel\\lib\\site-packages\\pyart\\io\\nexrad_cdm.py:141: RuntimeWarning: invalid value encountered in greater\n",
      "  azim_data[ray_i:end] = dvars[azimuth_var][var_index][:nradials]\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devel\\lib\\site-packages\\pyart\\io\\nexrad_cdm.py:142: RuntimeWarning: invalid value encountered in greater\n",
      "  elev_data[ray_i:end] = dvars[elevation_var][var_index][:nradials]\n",
      "C:\\Users\\matts\\Miniconda3\\envs\\devel\\lib\\site-packages\\pyart\\io\\nexrad_cdm.py:144: RuntimeWarning: invalid value encountered in greater\n",
      "  dvars[elevation_var][var_index][:nradials])\n"
     ]
    }
   ],
   "source": [
    "query = rs.query()\n",
    "#Here, set the initial time of the archived radar loop you want.\n",
    "dt = datetime(2017, 8, 26, 4, 10) # Our specified time\n",
    "query.stations('KHGX').time_range(dt, dt + timedelta(hours=.8))\n",
    "cat = rs.get_catalog(query)\n",
    "cat.datasets\n",
    "f = 27\n",
    "n = 1\n",
    "for item in sorted(cat.datasets.items()):\n",
    "    # After looping over the list of sorted datasets, pull the actual Dataset object out\n",
    "    # of our list of items and access over CDMRemote\n",
    "    try:\n",
    "        ds = item[1]\n",
    "        radar = pyart.io.nexrad_cdm.read_nexrad_cdm(ds.access_urls['OPENDAP'])\n",
    "        time_start = netCDF4.num2date(radar.time['data'][0], radar.time['units'])\n",
    "        #Now let's calculate and plot specific differential phase.\n",
    "        radar = radar.extract_sweeps([0])\n",
    "        kdp_dict = pyart.retrieve.kdp_proc.kdp_maesaka(radar)\n",
    "        radar.add_field('KDP', kdp_dict[0])\n",
    "\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",
    "\n",
    "        # exclude masked gates from the gridding\n",
    "        gatefilter = pyart.filters.GateFilter(radar)\n",
    "        gatefilter.exclude_masked('differential_reflectivity')\n",
    "\n",
    "        #Now let's grid the data\n",
    "        grid = pyart.map.grid_from_radars(\n",
    "            (radar,), gatefilters=(gatefilter, ),\n",
    "            grid_shape=(1, 500, 500),\n",
    "            grid_limits=((200, 200), (-123000.0, 123000.0), (-123000.0, 123000.0)),\n",
    "            fields=['differential_reflectivity','reflectivity','KDP'])\n",
    "\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",
    "\n",
    "    #Mask everything below 20dbz\n",
    "    #import numpy.ma as ma\n",
    "        ZDRmasked = ma.masked_where(REF < 20, ZDR)\n",
    "        REFmasked = ma.masked_where(REF < 20, REF)\n",
    "        KDPmasked = ma.masked_where(REF < 20, KDP)\n",
    "\n",
    "        rlons = grid.point_longitude['data']\n",
    "        rlats = grid.point_latitude['data']\n",
    "        cenlat = radar.latitude['data'][0]\n",
    "        cenlon = radar.longitude['data'][0]\n",
    "        #Let's set up the map projection!\n",
    "        # Set up our projection\n",
    "        crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)\n",
    "\n",
    "        # Set up our array of latitude and longitude values and transform to \n",
    "        # the desired projection.\n",
    "\n",
    "        tlatlons = crs.transform_points(ccrs.LambertConformal(central_longitude=265, central_latitude=25, standard_parallels=(25.,25.)),rlons[0,:,:],rlats[0,:,:])\n",
    "        tlons = tlatlons[:,:,0]\n",
    "        tlats = tlatlons[:,:,1]\n",
    "\n",
    "        # Limit the extent of the map area, must convert to proper coords.\n",
    "        LL = (cenlon-1,cenlat-1,ccrs.PlateCarree())\n",
    "        UR = (cenlon+1,cenlat+1,ccrs.PlateCarree())\n",
    "        print(LL)\n",
    "\n",
    "        # Get data to plot state and province boundaries\n",
    "        states_provinces = cfeature.NaturalEarthFeature(\n",
    "                category='cultural',\n",
    "                name='admin_1_states_provinces_lakes',\n",
    "                scale='50m',\n",
    "                facecolor='none')\n",
    "        fname = 'cb_2016_us_county_20m/cb_2016_us_county_20m.shp'\n",
    "        fname2 = 'cb_2016_us_state_20m/cb_2016_us_state_20m.shp'\n",
    "        counties = ShapelyFeature(Reader(fname).geometries(),ccrs.PlateCarree(), facecolor = 'none', edgecolor = 'black')\n",
    "        states = ShapelyFeature(Reader(fname2).geometries(),ccrs.PlateCarree(), facecolor = 'none', edgecolor = 'black')\n",
    "        fig=plt.figure(n,figsize=(30.,25.))\n",
    "        ax = plt.subplot(111,projection=ccrs.PlateCarree())\n",
    "        ax.coastlines('50m',edgecolor='black',linewidth=0.75)\n",
    "        #ax.add_feature(states_provinces,edgecolor='black',linewidth=0.5)\n",
    "        ax.add_feature(counties, edgecolor = 'black', linewidth = 0.5)\n",
    "        ax.add_feature(states, edgecolor = 'black', linewidth = 1.5)\n",
    "        ax.set_extent([LL[0],UR[0],LL[1],UR[1]])\n",
    "        REFlevels = np.arange(20,73,2)\n",
    "        ax.contourf(rlons[0,:,:],rlats[0,:,:],REFmasked,REFlevels,cmap = plt.cm.gist_ncar)\n",
    "        ZDRlevels = np.arange(3, np.max(ZDRmasked)+((np.max(ZDRmasked))-3)/2, (np.max(ZDRmasked))-3)\n",
    "        zdrlev = [1.5]\n",
    "        kdplev = [.75]\n",
    "        #ZDRlevels = np.arange(3,5.5,.5)\n",
    "        #ZDRlevels1 = np.arange(5,10,.5)\n",
    "        KDPlevels = np.arange(.75, np.max(KDPmasked)+((np.max(KDPmasked))-1.5)/2, (np.max(KDPmasked))-1.5)\n",
    "        #KDPlevels = np.arange(.75,1.75,.25)\n",
    "        #KDPlevels1 = np.arange(1.5,10,.25)\n",
    "        zdrc = ax.contour(rlons[0,:,:],rlats[0,:,:],ZDRmasked,zdrlev,linewidths = 5, colors='blue', alpha = .8)\n",
    "        #zrdc = ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked,ZDRlevels,linewide = .01, colors='pink', alpha = .8)\n",
    "        #ax.contourf(rlons[0,:,:],rlats[0,:,:],ZDRmasked,ZDRlevels1,linewide = .01, colors='crimson', alpha = .8)\n",
    "\n",
    "        #kdpc = ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels,linewide = .01, colors ='green', alpha = .5)\n",
    "        kdpc = ax.contour(rlons[0,:,:],rlats[0,:,:],KDPmasked,kdplev,linewidths = 5, colors='green', alpha = .8)\n",
    "        #ax.contourf(rlons[0,:,:],rlats[0,:,:],KDPmasked,KDPlevels1,linewide = .01, colors ='b', alpha = .5)\n",
    "        plt.savefig('testfig.png')\n",
    "\n",
    "        proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "                   pyproj.Proj(init='epsg:3857'))\n",
    "        zdr_areas = []\n",
    "        zdr_centroid_lon = []\n",
    "        zdr_centroid_lat = []\n",
    "        zdr_mean = []\n",
    "        zdr_max = []\n",
    "\n",
    "        for col in zdrc.collections:\n",
    "            # Loop through all polygons that have the same intensity level\n",
    "            print('hi')\n",
    "            for contour_path in col.get_paths(): \n",
    "                # Create the polygon for this intensity level\n",
    "                # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                print('hi')\n",
    "                for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "                    print('hi')\n",
    "                    cpa = np.asarray(cp[:])\n",
    "                    x = cpa[:,0]\n",
    "                    y = cpa[:,1]\n",
    "                    new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "                    if ncp == 0:\n",
    "                        poly = new_shape\n",
    "                        #print('hi')\n",
    "                    else:\n",
    "                        # Remove the holes if there are any\n",
    "                        poly = poly.difference(new_shape)\n",
    "                        #print('hi')\n",
    "\n",
    "                # do something with polygon\n",
    "                #print(poly.area) \n",
    "                #print(poly.centroid.x)\n",
    "                s_new = transform(proj, poly)\n",
    "                projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "                boundary = np.asarray(poly.boundary.xy)\n",
    "                polypath = Path(boundary.transpose())\n",
    "                coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "                mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "                mean = np.mean(ZDRmasked[mask])\n",
    "                if projected_area > 20 * units('km^2') and mean > 1.5:\n",
    "                    zdr_areas.append((projected_area))\n",
    "                    zdr_centroid_lon.append((poly.centroid.x))\n",
    "                    zdr_centroid_lat.append((poly.centroid.y))\n",
    "                    zdr_mean.append((mean))\n",
    "                    zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "                #print(s_new)\n",
    "\n",
    "        kdp_areas = []\n",
    "        kdp_centroid_lon = []\n",
    "        kdp_centroid_lat = []\n",
    "\n",
    "        for col in kdpc.collections:\n",
    "            # Loop through all polygons that have the same intensity level\n",
    "            for contour_path in col.get_paths(): \n",
    "                # Create the polygon for this intensity level\n",
    "                # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "                for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "                    print(1)\n",
    "                    cpa = np.asarray(cp[:])\n",
    "                    x = cpa[:,0]\n",
    "                    y = cpa[:,1]\n",
    "                    new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "                    if ncp == 0:\n",
    "                        poly = new_shape\n",
    "                    else:\n",
    "                        # Remove the holes if there are any\n",
    "                        poly = poly.difference(new_shape)\n",
    "\n",
    "                # do something with polygon\n",
    "                #print(poly.area) \n",
    "                print(poly.centroid.x)\n",
    "                s_new = transform(proj, poly)\n",
    "                projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "                if projected_area > 20 * units('km^2'):\n",
    "                    kdp_areas.append((projected_area))\n",
    "                    kdp_centroid_lon.append((poly.centroid.x))\n",
    "                    kdp_centroid_lat.append((poly.centroid.y))\n",
    "\n",
    "            #print(s_new)\n",
    "\n",
    "        if np.asarray(zdr_centroid_lon).shape[0] > 0:\n",
    "            ax.scatter(zdr_centroid_lon, zdr_centroid_lat, marker = '*', s = 500, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "        if np.asarray(kdp_centroid_lon).shape[0] > 0:\n",
    "            ax.scatter(kdp_centroid_lon, kdp_centroid_lat, marker = '^', s = 500, color = 'black', zorder = 10, transform=ccrs.PlateCarree())\n",
    "        for i in enumerate(zdr_areas):\n",
    "            plt.text(zdr_centroid_lon[i[0]]+.016, zdr_centroid_lat[i[0]]+.016, \"%.2f / %.2f / %.2f\" %(zdr_areas[i[0]].magnitude, zdr_mean[i[0]], zdr_max[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",
    "        title_plot = plt.title('Radar Reflectivity, ZDR, and KDP '+str(time_start.year)+'-'+str(time_start.month)+'-'+str(time_start.day)+\n",
    "                                   ' '+str(time_start.hour)+':'+str(time_start.minute)+' UTC', size = 25)\n",
    "        plt.savefig('ZDRArcFinderTropicalKHGX'+str(time_start.year)+str(time_start.month)+str(time_start.day)+str(time_start.hour)+'0'+str(time_start.minute)+'.png')\n",
    "        n = n+1\n",
    "    except:\n",
    "        continue\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POLYGON ((-87.14650269865113 32.82473736538209, -87.14654868829419 32.82917140303287, -87.14700313010604 32.8336093530376, -87.14618932442924 32.83803515542336, -87.14426907178432 32.84095094193749, -87.14348134826658 32.84244278025931, -87.14320003794761 32.84687367235605, -87.14400880658687 32.85131505357682, -87.14412653817635 32.85152746337874, -87.14523285494377 32.85576038189856, -87.14504025932081 32.86019213395635, -87.14588613693429 32.86463383537282, -87.14721510949792 32.86908016611732, -87.1484274403183 32.87352538031153, -87.14909014961773 32.87491189447746, -87.15042181754164 32.87797803027838, -87.15396458257069 32.88244542354081, -87.15426466374069 32.8827220207276, -87.15949797148092 32.88618561136618, -87.16476890509351 32.88684218515808, -87.1658554374341 32.88699164220916, -87.17003894155998 32.88757239816948, -87.17529764957479 32.88917123609809, -87.18057662112363 32.88923522403788, -87.18585960480704 32.88899139722734, -87.19115226210907 32.88799858675585, -87.19333091160763 32.88724750682975, -87.19645885414509 32.88591572226813, -87.20176595639896 32.88377323565189, -87.20484056497422 32.88291928394322, -87.20706345176815 32.88236220159904, -87.21236431652679 32.88067242034746, -87.2171489864354 32.87859720605191, -87.21767140876277 32.87847418672693, -87.2229661328569 32.87723664389679, -87.22827816471111 32.87460582630133, -87.22897032197187 32.87426954243271, -87.23358647497156 32.87224747867732, -87.23712897991336 32.86990839008801, -87.23891401366052 32.86831613063865, -87.24422266759683 32.86587544242866, -87.14650269865113 32.82473736538209))\n",
      "-87.18224253171091\n",
      "POLYGON ((-87.74302692893311 32.83995748540437, -87.74054689520469 32.84082647264969, -87.73630613878775 32.84140639525347, -87.74302692893311 32.83995748540437))\n",
      "-87.73995998764185\n",
      "POLYGON ((-87.72503217161567 32.85077025462511, -87.72993088587766 32.85002642591451, -87.73124493312385 32.85112278923191, -87.73192611462636 32.85555962914518, -87.73514796531846 32.85892970201235, -87.73551383684706 32.86001049997278, -87.74039320424068 32.86378501016927, -87.74132296867309 32.86447183436877, -87.7440825512044 32.86891848316366, -87.74562983391675 32.87002375432162, -87.75088660077932 32.87331443173031, -87.75617028928107 32.87253590545757, -87.75895859623324 32.87342162333026, -87.76143858388546 32.87412584708892, -87.76609775974293 32.87345440325699, -87.76672253283792 32.87329928542236, -87.76962994772994 32.87347044454507, -87.77199947493993 32.87356206011759, -87.77237612568162 32.87348287936892, -87.77728956778674 32.87173520248633, -87.78257477938934 32.87065772800467, -87.78785684074647 32.87007051805004, -87.79069392982348 32.86913055767409, -87.79314634629574 32.86825070526414, -87.72503217161567 32.85077025462511))\n",
      "-87.75896356916557\n",
      "POLYGON ((-88.22067889676072 32.87335852879227, -88.21540398305801 32.87055281970775, -88.21492992339337 32.87025151141264, -88.21013088521792 32.86665248408368, -88.20841980760434 32.86581178049768, -88.2048577423777 32.86311521856247, -88.19958187693 32.86171648662795, -88.19891243087633 32.86136861331879, -88.19430860064676 32.8587363173704, -88.1890329239211 32.85748915874174, -88.18660538052812 32.85692148949488, -88.18375696199556 32.85650510526452, -88.18078560910585 32.85691464247586, -88.17847815344787 32.85725324187909, -88.17752776996268 32.8569107007195, -88.17320446713538 32.85505869809333, -88.16792727696358 32.85492109829146, -88.16340916486028 32.85689257882072, -88.16264429302271 32.85785515031245, -88.15736060491807 32.86096888325765, -88.15717392494045 32.86131760315769, -88.15522642644268 32.86574838750609, -88.15278339313261 32.87017847222541, -88.15223158834937 32.87461122613837, -88.1520542016226 32.87518236400046, -88.15023670008614 32.87904186559921, -88.14817722333157 32.8834724033627, -88.14813638343921 32.88790587687471, -88.14873321111088 32.89234027786663, -88.14906978609959 32.89677430037719, -88.14865639376599 32.90120723202089, -88.14671292843556 32.90536886183532, -88.14656368180192 32.90563771441627, -88.14476244425931 32.91006855027825, -88.14256924302961 32.91449879894325, -88.14140864149473 32.9165179352555, -88.13875648828713 32.91892650861736, -88.13822929793633 32.92335922967307, -88.13945145747792 32.9277946455908, -88.14137959463926 32.93010751015478, -88.14310637206643 32.93223373543908, -88.14665023738641 32.935552453337, -88.14786433229148 32.93667434313303, -88.15192638520966 32.93859340107528, -88.15720667676031 32.93970477489206, -88.16248657962957 32.94109128551556, -88.16260357965328 32.94112863003026, -88.16776592986869 32.94286891377715, -88.17304566250868 32.9445516979557, -88.17832806009369 32.94480056918778, -88.18221413015551 32.94558700391069, -88.22067889676072 32.87335852879227))\n",
      "-88.17633344836422\n",
      "POLYGON ((-87.83005976682664 32.88255746593808, -87.83005229351033 32.88259658367549, -87.83016087785862 32.88703056307393, -87.83171489352327 32.8914703549651, -87.83292783497276 32.89590877575422, -87.8332585040412 32.90034364881826, -87.83509273804101 32.90478456891676, -87.83516652307259 32.90492641707409, -87.83825652941036 32.9092307036727, -87.84041123013418 32.91135548405396, -87.84401970415978 32.91368702485707, -87.84567528755416 32.91440991290116, -87.85094792333925 32.91595642964174, -87.8562227192317 32.9171363073872, -87.8595354626888 32.91818054332661, -87.86148683587965 32.92034695738199, -87.86675717280696 32.92244975983223, -87.86812509262087 32.92264645494408, -87.87203664712844 32.92284628840263, -87.87730620747507 32.92517750352729, -87.88044792303033 32.92449106550716, -87.83005976682664 32.88255746593808))\n",
      "-87.8499414047464\n",
      "POLYGON ((-87.89344306058982 32.90042889013652, -87.89610863212873 32.89614649946395, -87.89859473601027 32.89172168006163, -87.89859709578644 32.89171849816754, -87.90103972371678 32.88729661220633, -87.90333743344613 32.8828710333077, -87.90392546039421 32.88165980980619, -87.90638512472182 32.87844794773211, -87.90640556420968 32.87607473062466, -87.89344306058982 32.90042889013652))\n",
      "-87.9014539655917\n",
      "POLYGON ((-87.88971288155282 32.90794960580327, -87.89010656328166 32.90499238864849, -87.89327251144687 32.90076737091783, -87.88971288155282 32.90794960580327))\n",
      "-87.8910306520938\n",
      "POLYGON ((-88.58930952302948 32.98429639140587, -88.5853968369174 32.98530945587628, -88.58522566638199 32.98550236284176, -88.58509687033255 32.98531007968191, -88.57993519631663 32.983700679324, -88.57464269921776 32.98109897990521, -88.57446210972269 32.98089815800399, -88.56934630341507 32.9769622941972, -88.56883368501713 32.97647568601415, -88.58930952302948 32.98429639140587))\n",
      "-88.57960203902795\n",
      "POLYGON ((-88.53810042249621 33.00313280122178, -88.54298673848864 33.00328716793878, -88.54828060556824 33.00627587087785, -88.55356823552393 33.00669435936531, -88.55605042650664 33.00753459082564, -88.55886022601045 33.00875254254584, -88.56414976536041 33.00980621579078, -88.56943707226856 33.01000129719645, -88.57296851522311 33.01193583910433, -88.57473618072228 33.01438784562819, -88.57638957026079 33.01636256338337, -88.57806487574297 33.02079272571608, -88.57758658427619 33.0252272228934, -88.57722395734936 33.02966148731147, -88.5800771180412 33.03315394885674, -88.58090377340115 33.03408757859553, -88.58386998371235 33.03851501471777, -88.58538531184666 33.03997494245881, -88.5898182344383 33.04293612946682, -88.59068499098896 33.04367552669454, -88.59597787263928 33.04500422240799, -88.598835681222 33.04735026542517, -88.59941287587539 33.04891316515719, -88.53810042249621 33.00313280122178))\n",
      "-88.56392489775004\n",
      "POLYGON ((-88.54337850975783 33.16291357791224, -88.54262922134103 33.16273225030911, -88.538078592495 33.16130193528559, -88.53277970825158 33.1600476860208, -88.52748261018522 33.15950984857936, -88.52218679796684 33.15952182075924, -88.51689114342116 33.15960729645633, -88.51159830670784 33.16104107050672, -88.50630304968769 33.16137172823334, -88.50241856844191 33.16161402545112, -88.54337850975783 33.16291357791224))\n",
      "-88.52503103433354\n",
      "POLYGON ((-88.4904898131345 33.18051288673794, -88.49048048995796 33.18054563293225, -88.49045104178789 33.1806942130046, -88.48959722143508 33.18498030323893, -88.48843862214777 33.18941532662969, -88.48517621717495 33.19311559370769, -88.48442302573781 33.19385399403136, -88.47988625886809 33.19768113355391, -88.47867810942661 33.19829462086422, -88.47459224000488 33.20015253621195, -88.4692959875518 33.20142583251234, -88.46399612895138 33.20047619682258, -88.45882046939307 33.19831712290377, -88.4904898131345 33.18051288673794))\n",
      "-88.4764772027751\n",
      "POLYGON ((-88.42681948547234 33.21616948905598, -88.42692556571929 33.21619511003912, -88.42716226364013 33.21608049297264, -88.43221996476923 33.21208518241736, -88.43752103404864 33.21371691491388, -88.44078199944002 33.21606888461384, -88.44282544321678 33.21774624117701, -88.44812746442443 33.21967501064078, -88.45342306128384 33.21695285655034, -88.45872206862255 33.21673725344755, -88.4640246483589 33.21883751441074, -88.46666198621219 33.22047626162396, -88.46932870484348 33.22169338840031, -88.47463242222496 33.2241387345627, -88.47603758664535 33.22489898779202, -88.47993690899003 33.2268548742065, -88.48523951283943 33.2283393816372, -88.49036909107072 33.22931463592728, -88.48769214486127 33.23375161428946, -88.48615872698663 33.23818712090927, -88.48525966566123 33.23954379228002, -88.48234286821975 33.24262544047657, -88.4799670356884 33.2441911123856, -88.47727126504043 33.24706516853429, -88.4746753005868 33.2497085867119, -88.47276217119384 33.25150403687483, -88.46999133348184 33.25594077990423, -88.46938871066548 33.25882229156725, -88.46898863795121 33.26037545706862, -88.4693950489452 33.26274086595519, -88.46968138142725 33.26480820421161, -88.46940322893724 33.26779716221954, -88.4692264164819 33.26924225633788, -88.46940596844529 33.2694903645241, -88.4730506941935 33.27367136473354, -88.47471901120463 33.27574753028455, -88.47690111695375 33.27810034130529, -88.47645463227715 33.28253441051508, -88.47631905877556 33.28696810606904, -88.47628491504518 33.29140167971893, -88.48005494028446 33.29470578950417, -88.48277618436643 33.29582729283496, -88.48536662955732 33.29892781042292, -88.48750662283371 33.30025484422855, -88.48537095461357 33.30132599310932, -88.48007230244606 33.30467086205861, -88.47995018645936 33.30469787142056, -88.48007284470332 33.30498200384947, -88.48067786966622 33.30913050339528, -88.48008376202483 33.31124577626171, -88.47962206610605 33.31356533078727, -88.48008994612843 33.31479323914517, -88.48139657819601 33.31799667394202, -88.48254528151475 33.32242877622305, -88.48541498308447 33.32572578396167, -88.49072046674638 33.32546331779994, -88.49539913538283 33.32241213975128, -88.49602023471728 33.32225218214717, -88.50132609295844 33.3223006603739, -88.50663179349851 33.3222690843117, -88.50710966000322 33.32239580642909, -88.51193922479949 33.3230589240338, -88.51724572150466 33.32337381870924, -88.51984495771811 33.32237676221068, -88.52151985149996 33.31794063989894, -88.52009298346717 33.31350931223553, -88.51966869131587 33.30907643476488, -88.52251325760352 33.30622177057748, -88.42681948547234 33.21616948905598))\n",
      "-88.57658112699441\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POLYGON ((-88.09297598034306 33.25113103845684, -88.09306849802455 33.251362761465, -88.09466160264948 33.25579930366005, -88.09782969654424 33.26023882019854, -88.09825200268629 33.26069934043359, -88.1035439340432 33.26458269344208, -88.10364788675865 33.26468310515776, -88.10401862252995 33.26911730722803, -88.10352888139766 33.2703836938766, -88.10266735443197 33.27354836845205, -88.10036086525729 33.2779776455526, -88.09819794627087 33.2810501449646, -88.09695689196344 33.28240484350101, -88.09673742279992 33.28683796327235, -88.09724214437723 33.29127245089927, -88.0981683266091 33.29219408936279, -88.10346607046412 33.29457541552168, -88.10822242613163 33.29572608974299, -88.10876698771987 33.2958201217169, -88.10943521276408 33.29572824779169, -88.11407301738535 33.29506073053361, -88.11937936887304 33.2941293338919, -88.12468642105995 33.29285153882661, -88.12999365011534 33.29142819199448, -88.13076803897023 33.29133076211241, -88.13529954124846 33.29052621537562, -88.14060794555465 33.28841466575672, -88.14591269992306 33.28791176516663, -88.14946311114984 33.28692578200997, -88.151219759067 33.28625261674384, -88.15652862659846 33.28358021322299, -88.15790784497284 33.28250422766939, -88.16184251676692 33.2781273377196, -88.1618747789324 33.27807614222676, -88.16429511920001 33.2736458299901, -88.16704811205095 33.26921596010737, -88.16716240146775 33.2690965274211, -88.17047272000976 33.26478684177773, -88.17247738832849 33.26216959332362, -88.17376432439933 33.26035749188011, -88.1769520611503 33.25592792656575, -88.17779233307127 33.25477849510931, -88.1804701364687 33.25149865738016, -88.18310300004873 33.24939694420473, -88.09297598034306 33.25113103845684))\n",
      "-88.13252791592943\n",
      "POLYGON ((-88.51900773695624 33.337125741002, -88.51727752230798 33.33800841298723, -88.51335877123242 33.34012075506662, -88.51197598680217 33.3404663143883, -88.50667188215922 33.34181801702341, -88.5013916507724 33.34457159518824, -88.50137049259023 33.34461681650389, -88.49737381206164 33.34901065097079, -88.49607590104938 33.35111556874497, -88.4950157960116 33.35344737920232, -88.49608638673411 33.35654872171254, -88.49634778471277 33.35787912722316, -88.49703240527701 33.36231171922042, -88.49610092694454 33.36408071063759, -88.49336163075201 33.36675018021118, -88.49080074372512 33.36843445083228, -88.48549661356053 33.37090019853073, -88.4801878518338 33.37088792418096, -88.4794637298526 33.37120144335178, -88.47488314140635 33.37327717854336, -88.47166541886224 33.37564421781718, -88.46958073094827 33.37726161905501, -88.46688269546536 33.3800831819925, -88.46459231428695 33.38451927286368, -88.46428537868455 33.38606981731714, -88.46385328337865 33.38895361218694, -88.46429122952648 33.3898095903678, -88.4667571372723 33.3933839195293, -88.46961355461966 33.39745070192306, -88.47002016266418 33.39781378916263, -88.47492785856562 33.39978231178705, -88.48024094497325 33.40125405272585, -88.48555115798898 33.4010389893911, -88.4908628725301 33.40163982496486, -88.49616708863965 33.39832426514322, -88.49731960133077 33.3977795884582, -88.50147430621841 33.39671685964579, -88.50677858696234 33.39377665883425, -88.50869724869425 33.39333001892972, -88.51208694048701 33.39293015066296, -88.51739512671377 33.39205423429659, -88.51843197061362 33.39331548086577, -88.51884513697914 33.39774837517624, -88.51741103667712 33.39935664815081, -88.51432159891988 33.40218877305992, -88.51210878091291 33.40324409266768, -88.50680397933395 33.406125035507, -88.50580171267163 33.40663480575471, -88.50149819728051 33.40869127265049, -88.49927745874533 33.41107749516289, -88.49653135744724 33.41551480448241, -88.49672098500926 33.41994807675383, -88.49793609946809 33.42437993865855, -88.50153662628698 33.42794011045563, -88.50248013494461 33.42880718153021, -88.50685840268036 33.43257042039722, -88.50746461497666 33.43323361875275, -88.51900773695624 33.337125741002))\n",
      "-88.49575088435688\n",
      "POLYGON ((-88.53323121325434 33.35221747328099, -88.53103421012699 33.34896016549821, -88.52951245825831 33.34586853881778, -88.53323121325434 33.35221747328099))\n",
      "-88.53125929387988\n",
      "POLYGON ((-88.54385700400894 33.37550032963927, -88.54383281756897 33.3755397950253, -88.5390632525464 33.37998152490373, -88.53860686763394 33.38032381616269, -88.53782953910203 33.37998361162855, -88.53438180672603 33.37555585455105, -88.53367730322026 33.37112350111502, -88.53477717411252 33.36668812471739, -88.53500599310361 33.36225420768329, -88.53601630783137 33.35781898154028, -88.53457538669463 33.3545060921959, -88.54385700400894 33.37550032963927))\n",
      "-88.53797936123452\n",
      "POLYGON ((-88.16968337545705 33.44861071532345, -88.16764746161203 33.45099157497486, -88.16682003813146 33.45225623149663, -88.16481809683603 33.45542137283864, -88.16183495413242 33.45985093000178, -88.1614907579153 33.46027635239268, -88.15923134118687 33.46428089588476, -88.1587556277178 33.46871377435124, -88.16146604999165 33.47302024813709, -88.1667856393818 33.47058403525486, -88.17210243648888 33.46951999886416, -88.17443796977938 33.46873436395283, -88.17742036986125 33.46773323436391, -88.18170416797615 33.46430970683347, -88.18274265135017 33.46324704692108, -88.18582219133486 33.4598809991102, -88.18806642909126 33.45750343183771, -88.18923222335033 33.45545137469738, -88.19029976464442 33.45101903550598, -88.19148932889999 33.4465868326018, -88.19208904833442 33.44215397038543, -88.18809461096791 33.44021268490674, -88.18277920877547 33.44163635125626, -88.18165483535184 33.44214198413658, -88.17746304444981 33.44339460872027, -88.17400953158624 33.44501190789374, -88.16968337545705 33.44861071532345))\n",
      "-88.17576365080502\n",
      "POLYGON ((-88.30914148122388 33.5087310686432, -88.30600087136142 33.51316405979719, -88.30495943478822 33.51599370938972, -88.30396990219141 33.51759719928111, -88.30495882595426 33.51818743522649, -88.31027626921511 33.51980201159718, -88.3149546696403 33.51759893055549, -88.31559454857192 33.51728592962499, -88.32091238877248 33.51616104444523, -88.32623005962392 33.51602792486798, -88.32893482283623 33.51759969931249, -88.33154786251876 33.52070912006698, -88.33396575203325 33.52203312675504, -88.3351450181955 33.52646660747443, -88.33154807389619 33.52756166309928, -88.32622963594048 33.52976380305951, -88.32091116749359 33.52935911143227, -88.31814074342648 33.53089982651161, -88.3209107912872 33.53342325121298, -88.32622953316744 33.53309453700936, -88.33154824198522 33.53300946594151, -88.3368671372879 33.53505219625323, -88.33784628370509 33.53089998057973, -88.33756331801263 33.52646647346747, -88.34218440431694 33.52490586733458, -88.34750264856839 33.52508199540384, -88.30914148122388 33.5087310686432))\n",
      "-88.32384445185392\n",
      "POLYGON ((-88.35146464875379 33.51872682444157, -88.35150147297774 33.51759763864366, -88.35033895498047 33.51316430718882, -88.35044051861358 33.50873075824869, -88.34922283676039 33.50429743626135, -88.34749763823203 33.50186800874953, -88.34218065613749 33.50059260846759, -88.33686401794222 33.50133585460942, -88.33154731116907 33.50282621621425, -88.32623044674865 33.50347062532251, -88.32239329263967 33.50306503342659, -88.35146464875379 33.51872682444157))\n",
      "-88.34106638439488\n",
      "POLYGON ((-88.40086304486091 33.6416503083897, -88.40091003425589 33.64171763562762, -88.40078230213666 33.64183991309057, -88.3954568116277 33.64177317104365, -88.39539871247966 33.64176908731122, -88.40086304486091 33.6416503083897))\n",
      "-88.39900916160576\n",
      "POLYGON ((-88.39013790374906 33.65092677912148, -88.39499251868622 33.6550215284398, -88.3954673771093 33.65540466602292, -88.40079391245257 33.65570972924698, -88.40172512576049 33.65501773960852, -88.40611634720761 33.65141263234055, -88.40674376785465 33.65058115516474, -88.40831614836225 33.64614660771009, -88.4114360866848 33.64486451469448, -88.41676153334005 33.64468829164525, -88.42208746851529 33.64498228147637, -88.39013790374906 33.65092677912148))\n",
      "-88.3951944422728\n",
      "POLYGON ((-88.14234999467615 33.67491315264793, -88.13976489198583 33.67401722588798, -88.13443365653634 33.67569121733142, -88.13361722094008 33.67637453989949, -88.14234999467615 33.67491315264793))\n",
      "-88.1381445611236\n",
      "POLYGON ((-88.72238987769157 33.76843652384681, -88.72234614326992 33.76968537776953, -88.72624242119834 33.77302214245215, -88.73025700353541 33.76965996887477, -88.72838468868363 33.76523248200903, -88.72855543004104 33.76327135422841, -88.72238987769157 33.76843652384681))\n",
      "-88.72644205569841\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POLYGON ((-88.54143296330368 33.82333083902631, -88.54503057723434 33.82490284498289, -88.54616983683933 33.82332261821379, -88.54604412807129 33.82247568232404, -88.54143296330368 33.82333083902631))\n",
      "-88.54432913140134\n",
      "POLYGON ((-88.24642928224151 33.88112183072212, -88.2470939186511 33.87668872913978, -88.24611846153718 33.87607396403728, -88.24435668223161 33.87668688159805, -88.24077543415179 33.87897755308873, -88.23824169671789 33.88111608258338, -88.23606131227778 33.8855480150527, -88.23933144607707 33.88674242436772, -88.24642928224151 33.88112183072212))\n",
      "-88.24189931341498\n",
      "POLYGON ((-88.22590956316562 33.88997350915045, -88.22474305728649 33.88859755475809, -88.22445029806056 33.88997230015734, -88.22464390161267 33.89005329766938, -88.22590956316562 33.88997350915045))\n",
      "-88.22503247090957\n",
      "POLYGON ((-88.76957192675654 33.90113500799949, -88.76721520521798 33.90254080427768, -88.76765417881073 33.90414377117036, -88.76957192675654 33.90113500799949))\n",
      "-88.76814710359508\n",
      "POLYGON ((-88.19838723828222 33.90355623848879, -88.20114063458189 33.90325101113384, -88.19801727912269 33.90274265715955, -88.19838723828222 33.90355623848879))\n",
      "-88.19918171732891\n",
      "POLYGON ((-88.9210083701212 34.00736673019391, -88.9198678671503 34.00626551764785, -88.91804556470055 34.00388642540848, -88.91447563083855 33.99988524237018, -88.91138076252793 34.00391795814987, -88.90916771666409 34.00574061925931, -88.90383536688772 34.00809118228418, -88.9027461282916 34.00839183912124, -88.89850372550674 34.01059097605823, -88.89668897238946 34.01285330713612, -88.89605148394126 34.01728977557322, -88.89633025835549 34.02172204832453, -88.89519605073014 34.02616079184267, -88.89519419546187 34.03059434989319, -88.89555273124626 34.03502625739992, -88.89869113709194 34.03853876029966, -88.90212786724578 34.03499599429459, -88.90400603247622 34.03330883254722, -88.90700185542059 34.03053978222037, -88.90932747784208 34.02913263544973, -88.91466198361807 34.0269259762302, -88.9183336953363 34.02605280920159, -88.91999396844423 34.02439645432708, -88.9210083701212 34.00736673019391))\n",
      "-88.90779002622331\n",
      "POLYGON ((-89.11051314398728 34.38069759895857, -89.10961284039809 34.38411335582868, -89.10559587201102 34.38857216590757, -89.10521852052747 34.38908736091872, -89.10238183110005 34.39302580746229, -89.10003829002818 34.397474000918, -89.10054796182362 34.40190438690578, -89.10116342149574 34.40633411237214, -89.10077246379873 34.41077011971458, -89.10024911861477 34.41520695368123, -89.10009284446852 34.41599861781231, -89.09982992134563 34.41964312493307, -89.10015188301449 34.42242789545733, -89.10045629558056 34.42407279482565, -89.10396453737302 34.42848446000806, -89.10559741927439 34.43008006473692, -89.1077540787675 34.42846070762762, -89.11094110526901 34.42669121129527, -89.11194399370818 34.42400075735134, -89.11336350350734 34.41955820487021, -89.11404106938724 34.41512034888616, -89.1161942086011 34.41368521799622, -89.11051314398728 34.38069759895857))\n",
      "-89.10741693258166\n",
      "POLYGON ((-88.50128206668207 33.30015220474133, -88.5014186360759 33.30023622552528, -88.50128268875348 33.30046527895349, -88.50100418489893 33.30023680075912, -88.50128206668207 33.30015220474133))\n",
      "-88.50123511413493\n",
      "POLYGON ((-88.06628129406114 33.31315649954487, -88.06574600609778 33.313377065186, -88.06103085877017 33.31780047736504, -88.06520786520167 33.32224298018492, -88.06625275660132 33.32257986923486, -88.06678983728175 33.32224635365389, -88.06912175154956 33.31781772523922, -88.06810301557961 33.31338204703787, -88.06628129406114 33.31315649954487))\n",
      "-88.06564917842272\n",
      "POLYGON ((-88.02384097718418 33.31275020225528, -88.01853543210017 33.31284753483066, -88.01779825528946 33.31326564487983, -88.01321328567738 33.31750246730191, -88.01292826505427 33.31768683184392, -88.00912350712335 33.32211056610593, -88.0121980342107 33.32655202406363, -88.01317894775886 33.32693384443223, -88.01847706997788 33.32915227647502, -88.02377772259339 33.33072864659939, -88.02552834727041 33.33101899837438, -88.02908060058328 33.33171694723806, -88.03106208889093 33.33103245868332, -88.03439237939963 33.33012064706462, -88.03969816106664 33.33026324794213, -88.04500787345111 33.32921056200384, -88.04949129347298 33.32664198001757, -88.05032475034479 33.32588779185695, -88.0553591692296 33.32222158293101, -88.05034357083986 33.32002935837687, -88.04746180115382 33.3177702902066, -88.04504803664926 33.31694245447711, -88.03974062106005 33.31753197377119, -88.03444003357106 33.316088877413, -88.02914056113661 33.31437402432883, -88.02635317086281 33.31328687303194, -88.02384097718418 33.31275020225528))\n",
      "-88.0303170543531\n",
      "POLYGON ((-88.48541908753809 33.3279991687799, -88.48011228516172 33.32760355540422, -88.47480565969398 33.32728403753004, -88.46972134636776 33.33131114674265, -88.46950792034055 33.33241853125058, -88.46884825492712 33.3357456631393, -88.46704565353349 33.34018120826954, -88.4675541606084 33.34461417285075, -88.4695338749533 33.34841302409973, -88.47016068246388 33.34904476766791, -88.47484806082487 33.35246397714479, -88.47691719765616 33.35347037598627, -88.48015934713106 33.3545694029363, -88.4854668434395 33.35443493776303, -88.48682001188529 33.35345812208328, -88.49076819021082 33.35101803437352, -88.49492421452511 33.34901396839457, -88.49095945156542 33.34458571430804, -88.49075569595222 33.34433030574398, -88.48952725508147 33.3401540356807, -88.4895659114685 33.33572045321448, -88.48794872119095 33.33128900397461, -88.48541908753809 33.3279991687799))\n",
      "-88.47978103633328\n",
      "POLYGON ((-88.27303286738855 33.54782537363606, -88.26771493142994 33.54545148522166, -88.26239366468323 33.54777246342791, -88.26180226469187 33.54861623672446, -88.26180716408258 33.5530497712052, -88.262389329524 33.55338802690379, -88.26770927056785 33.55342231302919, -88.26870143663579 33.55305331318554, -88.27303183963353 33.54941037932836, -88.27366329775538 33.54862210273624, -88.27303286738855 33.54782537363606))\n",
      "-88.26695973459782\n",
      "POLYGON ((-88.13449152012116 33.65013937517555, -88.13378492453312 33.65045177113404, -88.13248277827466 33.6548832299585, -88.13432376649743 33.65931969700614, -88.13447007963826 33.65961034826805, -88.1379465205447 33.66375885339323, -88.13978272765328 33.66592164610492, -88.14129097121459 33.66376397397471, -88.14416509671673 33.65933477292602, -88.1444663489097 33.65490169359445, -88.13981233823364 33.6524755861928, -88.13525203337821 33.65045407591669, -88.13449152012116 33.65013937517555))\n",
      "-88.13852308787816\n",
      "POLYGON ((-88.80038671904454 33.67735130941325, -88.79971401332894 33.68074429040875, -88.80040936212518 33.6814723483364, -88.80115606603201 33.68073879192978, -88.80038671904454 33.67735130941325))\n",
      "-88.80042027106971\n",
      "POLYGON ((-88.73197119633548 33.85685728693437, -88.72662425877377 33.85512457268636, -88.72129714166854 33.85765684617674, -88.72036654348982 33.85836245899573, -88.72131387195661 33.86129635117435, -88.72208801606325 33.86279054387997, -88.72667181831616 33.86533156347846, -88.73200780960025 33.86461080863383, -88.73310624662176 33.86275504214747, -88.73258294418548 33.85832321724106, -88.73197119633548 33.85685728693437))\n",
      "-88.72719637420853\n",
      "POLYGON ((-88.59327555937895 33.8938354161389, -88.59151488783965 33.89417124661779, -88.59328679219658 33.89745843729338, -88.59369434928948 33.89416662005761, -88.59327555937895 33.8938354161389))\n",
      "-88.59283166657045\n",
      "POLYGON ((-88.57724976508239 33.89303614002128, -88.57191098125988 33.89385047287551, -88.57132195091351 33.89421246035428, -88.57191305289091 33.89457744833653, -88.57725773753138 33.89577362700559, -88.57882452142265 33.89419752490284, -88.57724976508239 33.89303614002128))\n",
      "-88.57553281831515\n",
      "POLYGON ((-88.74833379633279 33.92723839600181, -88.74399842373296 33.92922210220838, -88.74300447594017 33.93010963086803, -88.73766490477348 33.93091908610496, -88.736098792655 33.93368187567562, -88.73768461458448 33.93502633228505, -88.74303458642898 33.93630316089892, -88.74837018075245 33.93462784510459, -88.74973805887309 33.93363625798364, -88.75313398080294 33.92919113579921, -88.74833379633279 33.92723839600181))\n",
      "-88.74470324521567\n",
      "POLYGON ((-88.91978781800553 33.99474910816089, -88.91947030902371 33.99501255518623, -88.91980348752379 33.99700379497606, -88.92006891224115 33.99500969909545, -88.91978781800553 33.99474910816089))\n",
      "-88.91978029616179\n",
      "POLYGON ((-88.80967074238167 34.37927160194319, -88.80634863661975 34.38121882656681, -88.80432233068262 34.38342628089718, -88.80331021311984 34.38566404641599, -88.80124366605544 34.39010547221526, -88.80056908925251 34.39454158965241, -88.8026394353019 34.39896723775779, -88.80442249001946 34.40113232064056, -88.80980643094559 34.40299161551047, -88.81518030903308 34.40307064156537, -88.82053237392032 34.39941938453097, -88.82084240671985 34.39889635736213, -88.82216080683116 34.39445755429594, -88.82209079500991 34.39002428754071, -88.82045875564921 34.38683234723709, -88.81959499721243 34.3856006554154, -88.81505448586044 34.38131929241265, -88.81462584694582 34.38118668793344, -88.80967074238167 34.37927160194319))\n",
      "-88.81145598295235\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n",
      "[]\n",
      "[]\n",
      "[]\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "               pyproj.Proj(init='epsg:3857'))\n",
    "zdr_areas = []\n",
    "zdr_centroid_lon = []\n",
    "zdr_centroid_lat = []\n",
    "zdr_mean = []\n",
    "zdr_max = []\n",
    "\n",
    "for col in zdrc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "                print(new_shape)\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        boundary = np.asarray(poly.boundary.xy)\n",
    "        polypath = Path(boundary.transpose())\n",
    "        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "        mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "        mean = np.mean(ZDRmasked[mask])\n",
    "        if projected_area > 20 * units('km^2') and mean > 3.0:\n",
    "            zdr_areas.append((projected_area))\n",
    "            zdr_centroid_lon.append((poly.centroid.x))\n",
    "            zdr_centroid_lat.append((poly.centroid.y))\n",
    "            zdr_mean.append((mean))\n",
    "            zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "        #print(s_new)\n",
    "print(zdr_areas)\n",
    "print(zdr_centroid_lon)\n",
    "print(zdr_centroid_lat)\n",
    "print(zdr_mean)\n",
    "print(zdr_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.003594768552302712\n"
     ]
    }
   ],
   "source": [
    "for i in range(len(zdrc.collections)):\n",
    "    p = zdrc.collections[i].get_paths()[0]\n",
    "    v = p.vertices\n",
    "    x = v[:,0]\n",
    "    y = v[:,1]\n",
    "    if len(x)>2:\n",
    "        poly = sp.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "        print(i, poly.area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "print(len(zdrc.collections))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<enumerate object at 0x000001BAC8480B88>\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "print(enumerate(contour_path.to_polygons()))\n",
    "for cp in enumerate(contour_path.to_polygons()):\n",
    "    print(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "-88.2279510443226\n",
      "1\n",
      "-88.92762973301487\n",
      "1\n",
      "-88.90014267895808\n",
      "1\n",
      "-88.86370845890099\n",
      "1\n",
      "-89.2454798303056\n",
      "1\n",
      "-89.27887977518918\n",
      "[<Quantity(282.68846704195624, 'kilometer ** 2')>, <Quantity(50.18036435956034, 'kilometer ** 2')>, <Quantity(192.25203767886327, 'kilometer ** 2')>, <Quantity(39.52924508427346, 'kilometer ** 2')>]\n",
      "[-88.92762973301487, -88.86370845890099, -89.2454798303056, -89.27887977518918]\n",
      "[34.25133703035296, 34.45675677315477, 34.49475900749556, 34.62730337540038]\n"
     ]
    }
   ],
   "source": [
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "               pyproj.Proj(init='epsg:3857'))\n",
    "kdp_areas = []\n",
    "kdp_centroid_lon = []\n",
    "kdp_centroid_lat = []\n",
    "\n",
    "from shapely import geometry\n",
    "for col in kdpc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print(1)\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        if projected_area > 20 * units('km^2'):\n",
    "            kdp_areas.append((projected_area))\n",
    "            kdp_centroid_lon.append((poly.centroid.x))\n",
    "            kdp_centroid_lat.append((poly.centroid.y))\n",
    "\n",
    "        #print(s_new)\n",
    "print(kdp_areas)\n",
    "print(kdp_centroid_lon)\n",
    "print(kdp_centroid_lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "hi\n",
      "[]\n",
      "[]\n",
      "[]\n",
      "[]\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),\n",
    "               pyproj.Proj(init='epsg:3857'))\n",
    "zdr_areas = []\n",
    "zdr_centroid_lon = []\n",
    "zdr_centroid_lat = []\n",
    "zdr_mean = []\n",
    "zdr_max = []\n",
    "\n",
    "for col in zdrc.collections:\n",
    "    # Loop through all polygons that have the same intensity level\n",
    "    print('hi')\n",
    "    for contour_path in col.get_paths(): \n",
    "        # Create the polygon for this intensity level\n",
    "        # The first polygon in the path is the main one, the following ones are \"holes\"\n",
    "        print('hi')\n",
    "        for ncp,cp in enumerate(contour_path.to_polygons()):\n",
    "            print('hi')\n",
    "            cpa = np.asarray(cp[:])\n",
    "            x = cpa[:,0]\n",
    "            y = cpa[:,1]\n",
    "            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])\n",
    "            if ncp == 0:\n",
    "                poly = new_shape\n",
    "                print('hi')\n",
    "            else:\n",
    "                # Remove the holes if there are any\n",
    "                poly = poly.difference(new_shape)\n",
    "                print('hi')\n",
    "\n",
    "        # do something with polygon\n",
    "        #print(poly.area) \n",
    "        #print(poly.centroid.x)\n",
    "        s_new = transform(proj, poly)\n",
    "        projected_area = (transform(proj, poly).area * units('m^2')).to('km^2')\n",
    "        boundary = np.asarray(poly.boundary.xy)\n",
    "        polypath = Path(boundary.transpose())\n",
    "        coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T # create an Mx2 array listing all the coordinates in field\n",
    "        mask = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)\n",
    "        mean = np.mean(ZDRmasked[mask])\n",
    "        if projected_area > 20 * units('km^2') and mean > 3.0:\n",
    "            zdr_areas.append((projected_area))\n",
    "            zdr_centroid_lon.append((poly.centroid.x))\n",
    "            zdr_centroid_lat.append((poly.centroid.y))\n",
    "            zdr_mean.append((mean))\n",
    "            zdr_max.append((np.max(ZDRmasked[mask])))\n",
    "        #print(s_new)\n",
    "print(zdr_areas)\n",
    "print(zdr_centroid_lon)\n",
    "print(zdr_centroid_lat)\n",
    "print(zdr_mean)\n",
    "print(zdr_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}