<jupyter_start><jupyter_text>790 Project (ETC) Plotting ScriptsThis notebook contains all of the scripts used to generate figures for the ETC climatology study performed for GEOG 790 (SP 19)Required modules: Numpy, Matplotlib, Cartopy, Scipy, Shapely, GeopandasYou will also need descartes to plot the polygons from geopandas (conda install -c conda-forge descartes)<jupyter_code>import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.path as mplPath
from matplotlib.colors import Normalize
import cartopy
import cartopy.crs as ccrs
import scipy
from scipy import spatial
from scipy import stats
from shapely import geometry
from shapely.geometry import shape, Point, Polygon, LineString, MultiLineString
import pickle
import geopandas as gpd<jupyter_output><empty_output><jupyter_text>SettingsYou can adjust the naming convention, and the plot extents using these settings<jupyter_code>#plotExtent = [-125, -70, 23, 60] # Where to show the plot (LonMin, LonMax, LatMin, LatMax)
plotExtent = [-175, -5, 5, 85] #[-175, -10, 5, 85]

censusGrid = [-160, -40, 20, 70] # The grid used for storm_census.py

zoneColors = {
    "Clipper": "#3687FF",
    "Northwest": "#6C179E",
    "Colorado": "#F3A200",
    "GreatBasin": "#A45E45",
    "GulfOfMexico": "#4B9F0B",
    "EastCoast": "#DE55D9",
    "Other": "#666666",
}
typeStr = {
    "Clipper": "Alberta Clipper",
    "Northwest": "Northwestern",
    "Colorado": "Colorado Low",
    "GreatBasin": "Great Basin Low",
    "GulfOfMexico": "Gulf of Mexico Low",
    "EastCoast": "East Coast Low",
    "Other": "Other Lows",
}

clipper_zone = [[-115,50],
                [-105,50],
                [-105,55],
                [-110,55],
                [-110,60],
                [-125,60],
                [-125,55],
                [-115,55]]
northwst_zone = [[-125, 60],
                [-115,60],
                [-115,65],
                [-125,65]]
colorado_zone = [[-105, 40],
                [-100,40],
                [-100,35],
                [-105,35]]
basin_zone =    [[-120, 45],
                [-115,45],
                [-115,35],
                [-120,35]]
gom_zone =      [[-100, 25],
                [-90,25],
                [-90,30],
                [-100,30]]
eastcst_zone =  [[-80, 30],
                [-75,30],
                [-75,35],
                [-65,35],
                [-65,45],
                [-70,45],
                [-70,40],
                [-80,40]]<jupyter_output><empty_output><jupyter_text>Data LoadingRun the below cell block one time to load the data<jupyter_code>det_storms = np.load('storm_det_slp.npz', encoding='latin1', allow_pickle=True)
stormPosn = det_storms['storms']

tracked_storms = np.load('storm_track_slp.npz', encoding='latin1', allow_pickle=True)
storms = tracked_storms['storms']

census_storms = gpd.read_file("storm_census.geojson")

with open('narrdomain.data', 'rb') as filehandle:
    narr_list = pickle.load(filehandle)
with open('narrdomain_analysis.data', 'rb') as filehandle:
    narr_list_analysis = pickle.load(filehandle)   
    
states = cartopy.feature.NaturalEarthFeature(category='cultural',
                                            name='admin_1_states_provinces_lakes',
                                            scale='110m',
                                            facecolor='none')<jupyter_output><empty_output><jupyter_text>Helper Functions Run this code block to get all of the required helper functions needed, this is part of the original stormTracking git repository seen here: https://github.com/Phantom139/stormTracking/blob/master/storm_functions.py<jupyter_code>def calculate_bergeron(stormObject):
    bergeronList = []
    for i in range(len(stormObject['amp'])):
        pressures = stormObject['amp'][i:i+8]
        latitudes = stormObject['lat'][i:i+8]

        angFactor = np.sin(np.radians(60)) / np.sin(np.radians(latitudes))

        diffs = np.diff(pressures) / 100 #Convert Pa to mb
        negatives = np.sum(val for val in diffs if val < 0)	

        bFactor = (-1 * np.sum(negatives)) / (24)

        final = bFactor * angFactor

        bergeronList.append(bFactor)
    return bergeronList

def get_projection_object(type="Lambert"):
    # Cartopy has a "globe" object to define more projection standards, we create this first
    globe = ccrs.Globe(ellipse=None,
                       semimajor_axis=6370000,
                       semiminor_axis=6370000,
                       nadgrids="@null")    
    # Now we can create the projection object
    cen_lat  = float(50.0)
    cen_lon  = float(-107.0)
    std_pll  = [50.0]
    cutoff   = -30.0
    if type == "Lambert":
        projObj = ccrs.LambertConformal(central_latitude=cen_lat, 
                                        central_longitude=cen_lon,
                                        standard_parallels=std_pll,
                                        globe=globe,
                                        cutoff=cutoff)
    elif type == "EqualArea":
        projObj = ccrs.AlbersEqualArea(central_latitude=cen_lat,
                                       central_longitude=cen_lon,
                                       standard_parallels=std_pll,
                                       globe=globe)
    else:
        print("get_projection_object(): Invalid type " + type + " sent to function")
        return None
    return projObj<jupyter_output><empty_output><jupyter_text>Additional CalculationsThis code block handles additional pre-plotting calculations that may be required by multiple code blocks<jupyter_code>poly_clipper = Polygon(clipper_zone)
poly_northwst = Polygon(northwst_zone)
poly_colorado = Polygon(colorado_zone)
poly_basin = Polygon(basin_zone)
poly_gom = Polygon(gom_zone)
poly_eastcst = Polygon(eastcst_zone)

gdf_clipper = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[poly_clipper])
gdf_northwst = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[poly_northwst])
gdf_colorado = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[poly_colorado])
gdf_basin = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[poly_basin])
gdf_gom = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[poly_gom])
gdf_eastcst = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[poly_eastcst])

# Units of km ^2
clipper_area = gdf_clipper['geometry'].to_crs('epsg:3395').map(lambda p: p.area / 10**6)[0]
northwst_area = gdf_northwst['geometry'].to_crs('epsg:3395').map(lambda p: p.area / 10**6)[0]
colorado_area = gdf_colorado['geometry'].to_crs('epsg:3395').map(lambda p: p.area / 10**6)[0]
basin_area = gdf_basin['geometry'].to_crs('epsg:3395').map(lambda p: p.area / 10**6)[0]
gom_area = gdf_gom['geometry'].to_crs('epsg:3395').map(lambda p: p.area / 10**6)[0]
eastcst_area = gdf_eastcst['geometry'].to_crs('epsg:3395').map(lambda p: p.area / 10**6)[0]

print(f"Clipper Area {clipper_area} km^2")
print(f"Northwest Area {northwst_area} km^2")
print(f"Colorado Area {colorado_area} km^2")
print(f"Basin Area {basin_area} km^2")
print(f"Gulf Area {gom_area} km^2")
print(f"East Coast Area {eastcst_area} km^2")

class_area = {
    "Clipper": clipper_area,
    "Northwest": northwst_area,
    "Colorado": colorado_area,
    "GreatBasin": basin_area,
    "GulfOfMexico": gom_area,
    "EastCoast": eastcst_area
}<jupyter_output><empty_output><jupyter_text>Plotting Functions  Classification ZonesThis will plot a figure showing all of the classification zones, an optional parameter is provided to include the storm markers for all non-classified storms (Other Lows)<jupyter_code>def plot_zones(stormObj = None):
    print("plot_zones()")
      
    narr_domain = Polygon(narr_list)
    narr_domain_inner = Polygon(narr_list_analysis)
            
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    #ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent(plotExtent, crs=ccrs.PlateCarree())  
    #ax.set_xlim([-125, -70])
    #ax.set_ylim([23, 49])

    # Draw our plot, coastlines first, then the contours.
    ax.add_feature(states, edgecolor='k')    
    ax.coastlines()

    clipperPoly = mplPath.Path(clipper_zone)
    northwestPoly = mplPath.Path(northwst_zone)
    coloradoPoly = mplPath.Path(colorado_zone)
    basinPoly = mplPath.Path(basin_zone)
    gomPoly = mplPath.Path(gom_zone)
    eastCoastPoly = mplPath.Path(eastcst_zone)

    patch1 = matplotlib.patches.PathPatch(clipperPoly, facecolor=zoneColors["Clipper"], transform=ccrs.PlateCarree(), label="Alberta Clipper")
    patch2 = matplotlib.patches.PathPatch(northwestPoly, facecolor=zoneColors["Northwest"], transform=ccrs.PlateCarree(), label="Northwestern")
    patch3 = matplotlib.patches.PathPatch(coloradoPoly, facecolor=zoneColors["Colorado"], transform=ccrs.PlateCarree(), label="Colorado Low")
    patch4 = matplotlib.patches.PathPatch(basinPoly, facecolor=zoneColors["GreatBasin"], transform=ccrs.PlateCarree(), label="Great Basin Low")
    patch5 = matplotlib.patches.PathPatch(gomPoly, facecolor=zoneColors["GulfOfMexico"], transform=ccrs.PlateCarree(), label="Gulf of Mexico Low")
    patch6 = matplotlib.patches.PathPatch(eastCoastPoly, facecolor=zoneColors["EastCoast"], transform=ccrs.PlateCarree(), label="East Coast Low")
    ax.add_patch(patch1)
    ax.add_patch(patch2)
    ax.add_patch(patch3)
    ax.add_patch(patch4)
    ax.add_patch(patch5)
    ax.add_patch(patch6)
    
    #ax.add_geometries([narr_domain], crs=ccrs.PlateCarree(), facecolor = 'none', edgecolor='black', alpha=0.5, label="NARR Domain")
    # RF Note: Comment this line out to remove the NARR inner domain bounding box
    ax.add_geometries([narr_domain_inner], crs=ccrs.PlateCarree(), facecolor = 'none', edgecolor='red', alpha=0.5, label="Analysis Domain")

    if stormObj is not None:
        print("plot_zones(): stormObj is not None")
        print(" Entering loop " + str(len(stormObj)) + " points to process...")
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic'):
                lon, lat = stormObj[ed]['lon'][0], stormObj[ed]['lat'][0]
                if(lon < 0):
                    lon += 360
                stormPoint = Point(lon, lat)
                
                if (stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12]):
                    
                    if(narr_domain_inner.contains(stormPoint)):                
                        if (stormObj[ed]['classification'] == 'Other'):
                            plt.plot(lon, lat, 'bo', alpha=0.25, markeredgewidth=0, markersize=2, transform=ccrs.PlateCarree())
                        elif (stormObj[ed]['classification'] != 'Other'):
                            plt.plot(lon, lat, 'ro', alpha=0.25, markeredgewidth=0, markersize=2, transform=ccrs.PlateCarree())
                    else:
                        plt.plot(lon, lat, 'ko', alpha=0.25, markeredgewidth=0, markersize=2, transform=ccrs.PlateCarree())
                        
            if(ed % 500 == 0):
                print("%d / %d points processed (%.2f Percent)" % (ed, len(stormObj), ((ed/len(stormObj))*100)))
        plt.legend(framealpha = 1.0)
        plt.title('Storm Classification Zones (Other Marked)')
        plt.savefig('figures/storm_classes_with_other', bbox_inches='tight', pad_inches=0.05, dpi=300)
    else:
        print("plot_zones(): stormObj is None")
        plt.legend(framealpha = 1.0)
        plt.title('Storm Classification Zones')
        plt.savefig('figures/storm_classes', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_zones(): Done")
plot_zones()
plot_zones(storms)<jupyter_output><empty_output><jupyter_text>Track LinesPlots a "spaghetti-plot" style plot for all storm tracks, can be split for Bomb cyclones only with a single argument<jupyter_code>def plot_tracks(stormObj, bombOnly = False):
    print("plot_tracks(bombOnly = " + str(bombOnly) + ")")
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    #ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent(plotExtent, crs=ccrs.PlateCarree())   

    # Draw our plot, coastlines first, then the contours.
    ax.add_feature(states, edgecolor='k')    
    ax.coastlines()
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic'):
            if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12]):
                lon, lat = stormObj[ed]['lon'], stormObj[ed]['lat']
                pLon = []
                pLat = []
                
                prvLon = lon[0]
                for i, iL in enumerate(lon):
                    if(np.abs(iL - prvLon) > 180):
                        break
                    pLon.append(iL)
                    pLat.append(lat[i])
                
                bergeron = calculate_bergeron(stormObj[ed])
                if(bombOnly):
                    if(any(b >= 1 for b in bergeron)):                
                        plt.plot(pLon, pLat, 'r-', linewidth=0.5, alpha=0.35, transform=ccrs.PlateCarree())
                else:
                    plt.plot(pLon, pLat, 'r-', linewidth=0.5, alpha=0.35, transform=ccrs.PlateCarree())                    
    # Show the plot.
    plt.legend(framealpha = 1.0)
    plt.title("Storm Tracks (Bomb Cyclones)" if bombOnly == True else "Storm Tracks")
    plt.savefig("figures/bomb_tracks" if bombOnly == True else "figures/storm_tracks", bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_tracks(): Done")
#plot_tracks(storms, bombOnly=False) # RF Note: This takes a LONG time to run, so be mindful of that.
plot_tracks(storms, bombOnly=True)<jupyter_output><empty_output><jupyter_text>Mean Cyclone TrackGenerates a plot showing the mean line track for cyclones, can be split into bomb cyclones only<jupyter_code>def plot_mean_track(stormObj, bombOnly = False):
    print("plot_mean_track(bombOnly = " + str(bombOnly) + ")")
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    #ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent(plotExtent, crs=ccrs.PlateCarree())   

    # Draw our plot, coastlines first, then the contours.
    ax.add_feature(states, edgecolor='k')    
    ax.coastlines()

    # Process Variables
    monthList = [1, 2, 3, 4, 10, 11, 12]
    yearList = np.arange(1979, 2020, 1)
    cycloneList = ['Clipper', 'Northwest', 'Colorado', 'GreatBasin', 'GulfOfMexico', 'EastCoast']

    for classification in cycloneList:
        # Find the longest cyclone.
        longest = 0
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in monthList
                and stormObj[ed]['year'][0] in yearList and stormObj[ed]['classification'] == classification):
                count = len(stormObj[ed]['amp'])
                if(count > longest):
                    longest = count
        # Test classification
        meanLat = np.zeros((longest))
        meanLon = np.zeros((longest))
        count = np.zeros((longest))
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic'):
                # Test months
                if(storms[ed]['month'][0] in monthList):
                    # Test years
                    if(storms[ed]['year'][0] in yearList):
                        if(storms[ed]['classification'] == classification):
                            bergeron = calculate_bergeron(stormObj[ed])
                            if(bombOnly == True):
                                if(any(b >= 1 for b in bergeron)):	
                                    for i in range(len(stormObj[ed]['lat'])):
                                        meanLat[i] += stormObj[ed]['lat'][i] 
                                        meanLon[i] += stormObj[ed]['lon'][i]
                                        count[i] += 1
                            else:
                                for i in range(len(stormObj[ed]['lat'])):
                                    meanLat[i] += stormObj[ed]['lat'][i] 
                                    meanLon[i] += stormObj[ed]['lon'][i]
                                    count[i] += 1
        meanLat[:] = meanLat[:] / count[:]
        meanLon[:] = meanLon[:] / count[:]
        
        pLon = []
        pLat = []

        prvLon = meanLon[0]
        for i, iL in enumerate(meanLon):
            if(np.abs(iL - prvLon) > 180):
                break
            pLon.append(iL)
            pLat.append(meanLat[i])        

        print(classification + ": stdLat: %.3f" % np.std(meanLat))

        #print(meanLat)
        #print(meanLon)
        #print(count)
        # Draw the lines.
        plt.plot(pLon[0:35], pLat[0:35], linestyle='-', color=zoneColors[classification], alpha=1, markeredgewidth=0, transform=ccrs.PlateCarree())

    legendVals = list(typeStr.values())
    legendVals.remove("Other Lows")

    plt.legend(legendVals, loc=0, framealpha = 1.0)
    if(bombOnly == True):
        plt.title("Mean Bomb Cyclone Track (" + str(yearList[0]) + " - " + str(yearList[-1]) + ")")
        plt.savefig('figures/mean_track_bomb', bbox_inches='tight', pad_inches=0.05, dpi=300)
    else:
        plt.title("Mean Cyclone Track (" + str(yearList[0]) + " - " + str(yearList[-1]) + ")")
        plt.savefig('figures/mean_track', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_mean_track(): Done.")
plot_mean_track(storms, bombOnly = False)
plot_mean_track(storms, bombOnly = True)<jupyter_output><empty_output><jupyter_text>Four Panel<jupyter_code>def make_four_panel_tracks(stormObj):
    print("make_four_panel_tracks()")
    
    plotExtent4P = [-125, -70, 23, 60] #plotExtent
    
    
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    gs = fig.add_gridspec(2, 2)
    ax1 = fig.add_subplot(gs[0, 0], projection=projObj)
    ax2 = fig.add_subplot(gs[0, 1], projection=projObj)
    ax3 = fig.add_subplot(gs[1, 0], projection=projObj)
    ax4 = fig.add_subplot(gs[1, 1], projection=projObj)
    #ax = plt.axes(projection=ccrs.PlateCarree())
    
    monthList = [1, 2, 3, 4, 10, 11, 12]
    yearList = np.arange(1979, 2020, 1)
    cycloneList = ['Clipper', 'Northwest', 'Colorado', 'GreatBasin', 'GulfOfMexico', 'EastCoast']    
       
    
    print("Begin Plot 1 - All Cyclone Tracks")
    ax1.set_extent(plotExtent4P, crs=ccrs.PlateCarree())   
    ax1.add_feature(states, edgecolor='k', linewidth=0.5)    
    ax1.coastlines()
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic'):
            if(stormObj[ed]['month'][0] in monthList and stormObj[ed]['year'][0] in yearList):
                if(storms[ed]['classification'] in cycloneList):
                    lon, lat = stormObj[ed]['lon'], stormObj[ed]['lat']
                    #lon[lon <  0] = 360 + lon[lon < 0]
                    ax1.plot(lon, lat, color=zoneColors[stormObj[ed]['classification']], linewidth=0.2, alpha=0.75, transform=ccrs.PlateCarree())
                
    print("Plot 1 Done")
    print("Begin Plot 2 - Bomb Cyclone Tracks")
    ax2.set_extent(plotExtent4P, crs=ccrs.PlateCarree())   
    ax2.add_feature(states, edgecolor='k', linewidth=0.5)    
    ax2.coastlines()
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic'):
            if(stormObj[ed]['month'][0] in monthList and stormObj[ed]['year'][0] in yearList):
                if(storms[ed]['classification'] in cycloneList):
                    bergeron = calculate_bergeron(stormObj[ed])
                    if(any(b >= 1 for b in bergeron)):                 
                        lon, lat = stormObj[ed]['lon'], stormObj[ed]['lat']
                        #lon[lon <  0] = 360 + lon[lon < 0]
                        ax2.plot(lon, lat, color=zoneColors[stormObj[ed]['classification']], linewidth=0.2, alpha=0.75, transform=ccrs.PlateCarree())
                    
    print("Plot 2 Done")
    print("Begin Plot 3 - Mean Cyclone Track")

    ax3.set_extent(plotExtent4P, crs=ccrs.PlateCarree())   
    ax3.add_feature(states, edgecolor='k', linewidth=0.5)    
    ax3.coastlines()    
    
    for classification in cycloneList:
        # Find the longest cyclone.
        longest = 0
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in monthList
                and stormObj[ed]['year'][0] in yearList and stormObj[ed]['classification'] == classification):
                count = len(stormObj[ed]['amp'])
                if(count > longest):
                    longest = count
        # Test classification
        meanLat = np.zeros((longest))
        meanLon = np.zeros((longest))
        count = np.zeros((longest))
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic'):
                # Test months
                if(storms[ed]['month'][0] in monthList):
                    # Test years
                    if(storms[ed]['year'][0] in yearList):
                        if(storms[ed]['classification'] == classification):
                            for i in range(len(stormObj[ed]['lat'])):
                                meanLat[i] += stormObj[ed]['lat'][i] 
                                meanLon[i] += stormObj[ed]['lon'][i]
                                count[i] += 1
        meanLat[:] = meanLat[:] / count[:]
        meanLon[:] = meanLon[:] / count[:]

        ax3.plot(meanLon[::2], meanLat[::2], 
                 linestyle='-', color=zoneColors[classification], 
                 alpha=1, markeredgewidth=0, transform=ccrs.PlateCarree())    
    
    print("Plot 3 Done")
    print("Begin Plot 4 - Mean Bomb Cyclone Track")
    
    ax4.set_extent(plotExtent4P, crs=ccrs.PlateCarree())   
    ax4.add_feature(states, edgecolor='k', linewidth=0.5)    
    ax4.coastlines()    
    
    for classification in cycloneList:
        # Find the longest cyclone.
        longest = 0
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in monthList
                and stormObj[ed]['year'][0] in yearList and stormObj[ed]['classification'] == classification):
                count = len(stormObj[ed]['amp'])
                if(count > longest):
                    longest = count
        # Test classification
        meanLat = np.zeros((longest))
        meanLon = np.zeros((longest))
        count = np.zeros((longest))
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic'):
                # Test months
                if(storms[ed]['month'][0] in monthList):
                    # Test years
                    if(storms[ed]['year'][0] in yearList):
                        if(storms[ed]['classification'] == classification):
                            bergeron = calculate_bergeron(stormObj[ed])
                            if(any(b >= 1 for b in bergeron)):                         
                                for i in range(len(stormObj[ed]['lat'])):
                                    meanLat[i] += stormObj[ed]['lat'][i] 
                                    meanLon[i] += stormObj[ed]['lon'][i]
                                    count[i] += 1
        meanLat[:] = meanLat[:] / count[:]
        meanLon[:] = meanLon[:] / count[:]

        ax4.plot(meanLon, meanLat, linestyle='-', 
                 color=zoneColors[classification], alpha=1, 
                 markeredgewidth=0, transform=ccrs.PlateCarree())
        
        legendVals = list(typeStr.values())
        legendVals.remove("Other Lows")
        ax4.legend(legendVals, loc=0, framealpha = 1.0, prop={'size': 6})        
        
    
    print("Plot 4 Done")
    print("Done.")
make_four_panel_tracks(storms)<jupyter_output><empty_output><jupyter_text>Gridded Frequency PlotsPlots data using the storm_census counts<jupyter_code>print(census_storms)
def plot_census_frequencies(stormCensus):
    print("plot_census_frequencies()")

    projObj = get_projection_object(type="EqualArea")

    fig = plt.figure(figsize=(12,9))
    gs = fig.add_gridspec(2, 3)
    ax1 = fig.add_subplot(gs[0, 0], projection=projObj)
    ax2 = fig.add_subplot(gs[0, 1], projection=projObj)
    ax3 = fig.add_subplot(gs[0, 2], projection=projObj)    
    ax4 = fig.add_subplot(gs[1, 0], projection=projObj)
    ax5 = fig.add_subplot(gs[1, 1], projection=projObj)
    ax6 = fig.add_subplot(gs[1, 2], projection=projObj)      
    
    range1 = np.arange(1979, 2000, 1)
    range2 = np.arange(2000, 2020, 1)
    
    print("plot_census_frequencies(): Running pre-plot calculations")
    stormCensus["total_1"] = 0
    stormCensus["total_2"] = 0
    stormCensus["total_3"] = 0
    stormCensus["total_4"] = 0
    stormCensus["total_5"] = 0
    stormCensus["total_6"] = 0
    for t in range1:
        stormCensus["total_1"] += stormCensus["count_"+str(t)]
    for t in range2:
        stormCensus["total_2"] += stormCensus["count_"+str(t)]
    stormCensus["total_3"] = stormCensus["total_2"] - stormCensus["total_1"]
    for t in range1:
        stormCensus["total_4"] += stormCensus["count_bomb_"+str(t)]
    for t in range2:
        stormCensus["total_5"] += stormCensus["count_bomb_"+str(t)]
    stormCensus["total_6"] = stormCensus["total_5"] - stormCensus["total_4"] 
    
    crs = projObj
    crs_proj4 = crs.proj4_init
    censusProject = stormCensus.to_crs(crs_proj4)
    
    # Distribution of cyclone tracks and intensity
    print("plot_census_frequencies(): storm_track_distribution 1")
    #fig = plt.figure(figsize=(12,9))
    ax1.set_extent(plotExtent, crs=ccrs.PlateCarree())   
    ax1.add_feature(states, edgecolor='k')    
    ax1.coastlines()
    ax1.set(title = "1979 - 1999")   
    
    cmap1 = plt.cm.viridis
    cmap1.set_under(color='white')
    norm1 = Normalize(vmin=1, vmax=200)
    
    PCM1 = censusProject.plot(column='total_1', cmap=cmap1, norm=norm1, ax=ax1)
    

    print("plot_census_frequencies(): storm_track_distribution 2")
    ax2.set_extent(plotExtent, crs=ccrs.PlateCarree())   
    ax2.add_feature(states, edgecolor='k')    
    ax2.coastlines()  
    ax2.set(title = "2000 - 2019")
    
    censusProject.plot(column='total_2', cmap=cmap1, norm=norm1, ax=ax2)
    
    sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
    sm1.set_array([])
    plt.colorbar(sm1, ax=[ax1, ax2], orientation='horizontal', 
                 ticks=np.linspace(1,201,10), boundaries=np.arange(1,201,10))
    
    print("plot_census_frequencies(): storm_track difference")
 
    ax3.set_extent(plotExtent, crs=ccrs.PlateCarree())   
    ax3.add_feature(states, edgecolor='k')    
    ax3.coastlines()
    ax3.set(title = "20-Year Difference")
    
    cmap2 = plt.cm.RdBu
    norm2 = Normalize(vmin=-50, vmax=50)  
    
    censusProject.plot(column='total_3', cmap=cmap2, norm=norm2, ax=ax3)   
    
    #####
    #####
    #####
    
    print("plot_census_frequencies(): storm_track_distribution (bomb 1)")
    ax4.set_extent(plotExtent, crs=ccrs.PlateCarree())   
    ax4.add_feature(states, edgecolor='k')    
    ax4.coastlines()
    ax4.set(title = "1979 - 1999")
    
    cmap3 = plt.cm.viridis
    cmap3.set_under(color='white')
    norm3 = Normalize(vmin=1, vmax=30)
    
    censusProject.plot(column='total_4', cmap=cmap3, norm=norm3, ax=ax4) 

    print("plot_census_frequencies(): storm_track_distribution (bomb 2)")
    
    ax5.set_extent(plotExtent, crs=ccrs.PlateCarree())   
    ax5.add_feature(states, edgecolor='k')    
    ax5.coastlines()  
    ax5.set(title = "2000 - 2019")
    
    censusProject.plot(column='total_5', cmap=cmap3, norm=norm3, ax=ax5)   
        
    print("plot_census_frequencies(): storm_track difference bomb")
   
    ax6.set_extent(plotExtent, crs=ccrs.PlateCarree())   
    ax6.add_feature(states, edgecolor='k')    
    ax6.coastlines()
    ax6.set(title = "20-Year Difference")
    
    cmap4 = plt.cm.RdBu
    norm4 = Normalize(vmin=-10, vmax=10)  
    
    censusProject.plot(column='total_6', cmap=cmap4, norm=norm4, ax=ax6) 

    plt.suptitle('Distribution of Cyclone Tracks')
    plt.savefig('figures/panel_test.png', bbox_inches='tight', pad_inches=0.05, dpi=300)
plot_census_frequencies(census_storms)
print(census_storms['count_1979'])
def plot_theil_slopes(stormCensus):
    years = np.arange(1979, 2020, 1)

    census_vals = np.zeros((len(years), len(stormCensus['count_1979'])))

    slope = np.zeros(len(stormCensus['count_1979']))
    intercept = np.zeros(len(stormCensus['count_1979'])) 
    lo = np.zeros(len(stormCensus['count_1979']))
    hi = np.zeros(len(stormCensus['count_1979']))
    tau = np.zeros(len(stormCensus['count_1979']))
    pval = np.zeros(len(stormCensus['count_1979']))

    for i, y in enumerate(years):
        census_vals[i] = stormCensus['count_' + str(y)]

    for x in range(len(stormCensus['count_1979'])):
        slope[x], intercept[x], lo[x], hi[x] = scipy.stats.mstats.theilslopes(census_vals[:, x], years)
        tau[x], pval[x] = scipy.stats.kendalltau(years, census_vals[:, x])
        
    stormCensus['theil_slope'] = slope
    stormCensus['theil_inter'] = intercept
    stormCensus['theil_lo'] = lo
    stormCensus['theil_up'] = hi
    stormCensus['kendall_tau'] = tau
    stormCensus['kendall_p'] = pval

    stormCensus['isSig'] = pval < 0.1
    stormCensus['isSig'] = stormCensus['isSig'].astype(int)
    stormCensus['isSig'] = stormCensus['isSig'].replace({0:np.nan})
    
    print(stormCensus['isSig'].value_counts())
    print(len(stormCensus['isSig']))
    
    projObj = get_projection_object()     
    crs = projObj
    crs_proj4 = crs.proj4_init
    censusProject = stormCensus.to_crs(crs_proj4)    
   
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    #ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent(plotExtent, crs=ccrs.PlateCarree())   

    # Draw our plot, coastlines first, then the contours.
    ax.add_feature(states, edgecolor='k')    
    ax.coastlines()    
    
    cmap = plt.cm.RdBu
    norm = Normalize(vmin=-0.25, vmax=0.25)  
    
    censusProject.plot(column='theil_slope', cmap=cmap, norm=norm, ax=ax, legend=True)
    censusProject.plot(column='isSig', vmin=0.01, hatch='...', 
                       facecolor=None, edgecolor=None, ax=ax, alpha=0.01)
plot_theil_slopes(census_storms)<jupyter_output><empty_output><jupyter_text>Bar Plots Bomb Frequency by Type<jupyter_code>def plot_bomb_frequency_typed(stormObj):
    print("plot_bomb_frequency_typed()")
    fig = plt.figure(figsize=(12,9))
    yrs = np.arange(1979, 2020, 1)
    totals_year = np.zeros((len(yrs)))
    
    classList = list(typeStr.keys())
    classList.remove("Other")
    
    for className in classList:
        count_list_ind = np.zeros((len(yrs)))
        print(className)
        for y in yrs:
            indCount = 0
            for ed in range(len(stormObj)):
                
                if (stormObj[ed]['type'] == 'cyclonic'):
                    if(stormObj[ed]['classification'] == className):
                        if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12] and stormObj[ed]['year'][0] == y):
                            bergeron = calculate_bergeron(stormObj[ed])
                            if(any(b >= 1 for b in bergeron)):                            
                                indCount += 1
            count_list_ind[y - 1979] = indCount
            print(str(y) + ": " + str(indCount))
        #Draw the plot
        plt.bar(yrs, count_list_ind, bottom=totals_year, color=zoneColors[className])
        totals_year += count_list_ind
    plt.ylim(0, 20)
    plt.xlabel("Year")
    plt.ylabel("Number of Cyclones")
    plt.legend(typeStr.values(), loc=0, framealpha = 1.0)
    plt.title("Number of Bomb Cyclones by Cyclone Type")
    plt.savefig('figures/bomb_frequency_typed', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_bomb_frequency_typed(): Done")
plot_bomb_frequency_typed(storms)<jupyter_output><empty_output><jupyter_text>Frequencies of all cyclones by cyclone typeThis plot is for "all" cyclones per year with "Other Lows" removed due to the large magnitude of difference<jupyter_code>def plot_frequency_typed(stormObj):
    print("plot_frequency_typed()")
    fig = plt.figure(figsize=(12,9))
    yrs = np.arange(1979, 2020, 1)
    totals_year = np.zeros((len(yrs)))
    totals_class = {}
    for className in typeStr.keys():
        count_list_ind = np.zeros((len(yrs)))
        for y in yrs:
            indCount = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == className and className != "Other"):
                    if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12] and stormObj[ed]['year'][0] == y):
                        indCount += 1
            count_list_ind[y - 1979] = indCount
        totals_class[className] = count_list_ind
        #Draw the plot
        plt.bar(yrs, count_list_ind, bottom=totals_year, color=zoneColors[className])
        totals_year += count_list_ind

    legendList = list(typeStr.values())
    legendList.remove("Other Lows")

    plt.legend(legendList, loc=0, framealpha = 1.0)
    plt.xlabel("Year")
    plt.ylabel("Number of Cyclones")
    plt.title("Number of Cyclones by Cyclone Type")
    plt.savefig('figures/frequency_typed', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_frequency_typed(): Done")
    
    print("\nSummary Statistics:")
    print("Totals By Year:")
    for i, t in enumerate(totals_year):
        print(str(i+1979) + ": " + str(t))
    print("Mean: " + str(np.average(totals_year)))
    print("Max: " + str(np.max(totals_year)))
    print("Min: " + str(np.min(totals_year)))
    print("\nStatistics By Class:")
    for className in typeStr.keys():
        print(className)
        for i, t in enumerate(totals_class[className]):
            print(str(i+1979) + ": " + str(t))
        print("Mean: " + str(np.average(totals_class[className])))
        print("Max: " + str(np.max(totals_class[className])))
        print("Min: " + str(np.min(totals_class[className])))
        print("\n")
plot_frequency_typed(storms)<jupyter_output><empty_output><jupyter_text>Cyclone DensityPlots a scatter plot of cyclone count density for each year.<jupyter_code>def plot_count_density(stormObj):
    print("plot_count_density()")
    
    yrs = np.arange(1979, 2020, 1)
    totals_year = {}
    
    classes = ["Clipper", "Colorado", "Northwest", "GreatBasin", "GulfOfMexico", "EastCoast"]
    for className in classes:
        count_list_ind = np.zeros((len(yrs)))
        for y in yrs:
            indCount = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == className and className != "Other"):
                    if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12] and stormObj[ed]['year'][0] == y):
                        indCount += 1
            count_list_ind[y - 1979] = (indCount / class_area[className]) * 10e6
        totals_year[className] = count_list_ind
        
    fig = plt.figure(figsize=(12, 9))
    #for className in classes:
    #    plt.scatter(yrs, totals_year[className], color=zoneColors[className])
    
    #pct_year = {}
    #for className in classes:
    #    pct_year[className] = np.zeros((len(yrs)))
    for i, y in enumerate(yrs):
        pct_year = 0
        total = 0
        # Fetch the total # of cyclones for the year
        for className in classes:
            total += totals_year[className][i]
        # Calculate the percentage for each zone
        for className in classes:
            pct_class = (totals_year[className][i] / total) * 100
            plt.bar(y, pct_class, bottom=pct_year, color=zoneColors[className])
            pct_year += pct_class
    
    print(pct_year)
    plt.legend(classes, loc=0, framealpha = 1.0)
    plt.xlabel("Year")
    plt.ylabel("Percent")
    #plt.ylabel("Cyclone Count Density ((count / $km^2$) * $10^-6$)")
    plt.title("Percent of Classified Cyclones by Zone (Normalized by Area)")
    
    print("plot_count_density(): Done")
plot_count_density(storms)<jupyter_output><empty_output><jupyter_text>Three-Panel BarsCombines the three above plots into a single plot<jupyter_code>def three_panel_bar_stats(stormObj):
    print("three_Panel_bar_stats()")
    fig = plt.figure(figsize=(12, 9))
    
    gs = fig.add_gridspec(3, 1)
    ax1 = fig.add_subplot(gs[0, 0])
    ax2 = fig.add_subplot(gs[1, 0]) 
    ax3 = fig.add_subplot(gs[2, 0])
    
    print("three_Panel_bar_stats(): Start Panel A")
    
    classes = ["Clipper", "Colorado", "Northwest", "GreatBasin", "GulfOfMexico", "EastCoast"]
    
    yrs = np.arange(1979, 2020, 1)
    totals_year = np.zeros((len(yrs)))
    for className in classes:
        count_list_ind = np.zeros((len(yrs)))
        for y in yrs:
            indCount = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == className):
                    if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12] and stormObj[ed]['year'][0] == y):
                        indCount += 1
            count_list_ind[y - 1979] = indCount
        #Draw the plot
        ax1.bar(yrs, count_list_ind, bottom=totals_year, color=zoneColors[className])
        totals_year += count_list_ind

    legendList = list(typeStr.values())
    legendList.remove("Other Lows")

    ax1.legend(legendList, loc=0, framealpha = 1.0)
    ax1.set(xlabel="Year", ylabel="Number of Cyclones", title="Cyclogenesis Frequency by Region")  
    
    print("three_Panel_bar_stats(): Start Panel B")
    totals_year = np.zeros((len(yrs)))
    for className in classes:
        count_list_ind = np.zeros((len(yrs)))
        for y in yrs:
            indCount = 0
            for ed in range(len(stormObj)):
                
                if (stormObj[ed]['type'] == 'cyclonic'):
                    if(stormObj[ed]['classification'] == className):
                        if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12] and stormObj[ed]['year'][0] == y):
                            bergeron = calculate_bergeron(stormObj[ed])
                            if(any(b >= 1 for b in bergeron)):                            
                                indCount += 1
            count_list_ind[y - 1979] = indCount
        #Draw the plot
        ax2.bar(yrs, count_list_ind, bottom=totals_year, color=zoneColors[className])
        totals_year += count_list_ind
    ax2.set_ylim(0, 14)  
    ax2.set(xlabel="Year", ylabel="Number of Bomb Cyclones", title="Bomb Cyclone Frequency by Region") 
    
    print("three_Panel_bar_stats(): Start Panel C")
    
    totals_year = {}
    for className in classes:
        count_list_ind = np.zeros((len(yrs)))
        for y in yrs:
            indCount = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == className):
                    if(stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12] and stormObj[ed]['year'][0] == y):
                        indCount += 1
            count_list_ind[y - 1979] = (indCount / class_area[className]) * 10e6
        totals_year[className] = count_list_ind  
        
    for i, y in enumerate(yrs):
        pct_year = 0
        total = 0
        # Fetch the total # of cyclones for the year
        for className in classes:
            total += totals_year[className][i]
        # Calculate the percentage for each zone
        for className in classes:
            pct_class = (totals_year[className][i] / total) * 100
            ax3.bar(y, pct_class, bottom=pct_year, color=zoneColors[className])
            pct_year += pct_class     
            
    ax3.set(xlabel="Year", ylabel="Percent", title="Percent of Classified Cyclones by Zone (Normalized by Area)")
    
    print("three_Panel_bar_stats(): Done")
three_panel_bar_stats(storms)<jupyter_output><empty_output><jupyter_text>Comparative Bar ChartCompares the number of classified cyclones to the number of "other" cyclones, non-bomb and bombs sorted into two sub-plots<jupyter_code>def plot_comparative_counts(stormObj):
    print("plot_comparative_counts()")
    
    yrs = np.arange(1979, 2020, 1)
    
    total_classified_year = np.zeros((len(yrs)))
    total_classified_bomb_year = np.zeros((len(yrs)))
    total_nonclassified_year = np.zeros((len(yrs)))
    total_nonclassified_bomb_year = np.zeros((len(yrs)))
    
    for ed in range(len(stormObj)):
        if(stormObj[ed]['year'][0] in yrs 
           and stormObj[ed]['type'] == 'cyclonic' 
           and stormObj[ed]['month'][0] in [1, 2, 3, 4, 10, 11, 12]):
            isBomb = False
            isClassified = False

            yI = int(stormObj[ed]['year'][0] - 1979)           
            bergeron = calculate_bergeron(stormObj[ed])
            if(any(b >= 1 for b in bergeron)):            
                isBomb = True
            if(stormObj[ed]['classification'] != 'Other'):
                isClassified = True

            if(isClassified):
                if(isBomb):
                    total_classified_bomb_year[yI] += 1
                total_classified_year[yI] += 1
            else:
                if(isBomb):
                    total_nonclassified_bomb_year[yI] += 1
                total_nonclassified_year[yI] += 1    
    
    fig = plt.figure(figsize=(12, 9))
    
    gs = fig.add_gridspec(1, 2)
    ax1 = fig.add_subplot(gs[0, 0])
    ax2 = fig.add_subplot(gs[0, 1])    
    
    ax1.bar(yrs, total_classified_year, color=zoneColors['Clipper'], label="Classified Lows")
    ax1.bar(yrs, total_nonclassified_year, bottom=total_classified_year, color=zoneColors['Other'], label="Other Lows")
    ax1.set(xlabel='Year', ylabel='Count', title='All Cyclones')
    
    ax1.legend()
    
    ax2.bar(yrs, total_classified_bomb_year, color=zoneColors['Clipper'])
    ax2.bar(yrs, total_nonclassified_bomb_year, bottom=total_classified_bomb_year, color=zoneColors['Other'])
    ax2.set(xlabel='Year', ylabel='Count', title='Bomb Cyclones')
    
    plt.suptitle("Comparison of Classified vs. Other Cyclones")
    
    print("Storm Counts")
    for i in range(len(total_classified_year)):
        print(str(i + 1979) + ": C: " + str(total_classified_year[i]) + ", NC: " + str(total_nonclassified_year[i]))
    print("Bomb Counts")
    for i in range(len(total_classified_year)):
        print(str(i + 1979) + ": C: " + str(total_classified_bomb_year[i]) + ", NC: " + str(total_nonclassified_bomb_year[i]))        
    
    print("Stats")
    print("Other Max: " + str(np.max(total_nonclassified_year)))
    print("Other Min: " + str(np.min(total_nonclassified_year)))
    print("Other Avg: " + str(np.average(total_nonclassified_year)))
    print("Other Bomb Max: " + str(np.max(total_nonclassified_bomb_year)))
    print("Other Bomb Min: " + str(np.min(total_nonclassified_bomb_year)))
    print("Other Bomb Avg: " + str(np.average(total_nonclassified_bomb_year)))
    
    print("plot_comparative_counts(): Done")
plot_comparative_counts(storms)<jupyter_output><empty_output><jupyter_text>Bar Chart of Time of Maximum Pressure FallPlots a bar chart of how many hours after cyclogenesis the maximum pressure fall occurs.<jupyter_code>def plot_strength_bars(stormObj):
    print("plot_strength_bars()")
    fig = plt.figure(figsize=(12,9))

    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    dt = 3

    # Find the longest cyclone.
    longest = 0
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            count = len(stormObj[ed]['amp'])
            if(count > longest):
                longest = count

    #totals_column = np.zeros((longest))
    for className in typeStr.keys():
        fall_array = np.zeros((longest))
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years and stormObj[ed]['classification'] == className):
                dPressure = np.diff(stormObj[ed]['amp'])
                largestFall = np.nanmin(dPressure)
                index_of_largest_fall = np.argmin(dPressure)
                time_of_fall = index_of_largest_fall * dt
                fall_array[index_of_largest_fall] += 1
        #idx_of_last_non_zero = np.max(np.nonzero(fall_array))
        #splitArray = fall_array[0:idx_of_last_non_zero+1]
        times = np.arange(0, (len(fall_array)*3), 3)
        
        plt.plot(times+1, fall_array, 'o', color=zoneColors[className])
        
        finiteYmask = np.isfinite(np.array(fall_array))
        Yclean = np.array(fall_array)[finiteYmask]
        Xclean = np.array(times)[finiteYmask]

        slope, intercept, r_value, p_value, std_err = stats.linregress(Xclean, Yclean)
        line = slope*(times)+intercept
        plt.plot((times), line, color=zoneColors[className], linestyle='--') 
        
        print("Statistics " + className + ": m: %.3f, b: %.3f, r: %.3f, p: %.3f, stdErr: %.3f, stdDev: %.3f" % (slope, intercept, r_value, p_value, std_err, np.std(Yclean)))

        #plt.bar(times+1, fall_array, bottom=totals_column, width=2, color=zoneColors[className])
        #totals_column += fall_array
        
    for type in typeStr.keys():
        plt.plot(0, 0, 'o', alpha=1, color=zoneColors[type], label=typeStr[type])        
        
    idx_of_last_non_zero = np.max(np.nonzero(fall_array))
    plt.xlim(0, idx_of_last_non_zero*3)
    plt.legend(loc=1, framealpha = 1.0)
    plt.xlabel("Time after cyclogenesis (hours)")
    plt.ylabel("Count")
    plt.yscale('log')
    plt.ylim(-0.001, 3000)
    plt.title("Time of Maximum Pressure Fall")
    plt.savefig('figures/pressure_fall_occurance', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_strength_bars(): Done")
plot_strength_bars(storms)<jupyter_output><empty_output><jupyter_text>Bar chart of maximum pressure fall (Bombs only)Same as above, but split for bomb cyclones only<jupyter_code>def plot_strength_bars_bombs(stormObj):
    print("plot_strength_bars_bombs()")
    fig = plt.figure(figsize=(12,9))

    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    dt = 3

    # Find the longest cyclone.
    longest = 0
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            count = len(stormObj[ed]['amp'])
            if(count > longest):
                longest = count

    totals_column = np.zeros((longest))
    for className in typeStr.keys():
        fall_array = np.zeros((longest))
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years and stormObj[ed]['classification'] == className):
                bergeron = calculate_bergeron(stormObj[ed])
                if(any(b >= 1 for b in bergeron)):
                    dPressure = np.diff(stormObj[ed]['amp'])
                    largestFall = np.nanmin(dPressure)
                    index_of_largest_fall = np.argmin(dPressure)
                    time_of_fall = index_of_largest_fall * dt
                    fall_array[index_of_largest_fall] += 1
        times = np.arange(0, (len(fall_array)*3), 3)
        plt.bar(times+1, fall_array, bottom=totals_column, width=2, color=zoneColors[className])
        totals_column += fall_array
    idx_of_last_non_zero = np.max(np.nonzero(fall_array))
    plt.xlim(0, idx_of_last_non_zero*3)
    plt.legend(typeStr.values(), loc=0, framealpha = 1.0)
    plt.xlabel("Time after cyclogenesis (hours)")
    plt.ylabel("Count")
    plt.title("Time of Maximum 3-Hour Pressure Fall (Bomb Cyclones Only)")
    plt.savefig('figures/pressure_fall_bombs_occurance', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_strength_bars_bombs(): Done")
plot_strength_bars_bombs(storms)<jupyter_output><empty_output><jupyter_text>Mean Plots w/ Statistics Mean Cyclogenesis LatitudeGenerates a plot of mean cyclogenesis latitude, prints output of basic statistics<jupyter_code>def plot_mean_lat(stormObj):
    print("plot_mean_lat()")
    fig = plt.figure(figsize=(12,9))

    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    cyclGroups = ["Clipper", "Colorado", "Northwest", "GreatBasin", "GulfOfMexico", "EastCoast"]

    for year in years:
        for type in cyclGroups:
            meanLat = 0
            count = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == type and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] == year):
                    lat = stormObj[ed]['lat'][0]
                    meanLat += lat
                    count += 1
            if(count != 0):
                mLat = meanLat / count
                plt.plot(year, mLat, '-o', color=zoneColors[type])
    # Regression Lines
    for type in cyclGroups:
        latList = []
        for year in years:
            lat = 0
            meanLat = 0
            mLat = 0
            count = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == type and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] == year):
                    lat = stormObj[ed]['lat'][0]
                    meanLat += lat
                    count += 1
            if count != 0:
                mLat = meanLat / count
            else:
                mLat = np.nan

            latList.append(mLat)
            print(type + " (" + str(year) + "): " + str(mLat))

        finiteYmask = np.isfinite(np.array(latList))
        Yclean = np.array(latList)[finiteYmask]
        Xclean = np.array(years)[finiteYmask]

        slope, intercept, r_value, p_value, std_err = stats.linregress(Xclean, Yclean)
        line = slope*years+intercept
        plt.plot(years, line, color=zoneColors[type], linestyle='--')

        print("Statistics " + type + ": m: %.3f, b: %.3f, r: %.3f, p: %.3f, stdErr: %.3f, stdDev: %.3f" % (slope, intercept, r_value, p_value, std_err, np.std(Yclean)))
    # Add "Fake" Lines for Legend
    for type in cyclGroups:	
        plt.plot(0, 0, linestyle='-', alpha=1, color=zoneColors[type], label=typeStr[type])
    plt.xlim(1979, 2020)
    plt.xlabel("Year")
    plt.ylim(20, 70)
    plt.ylabel("Latitude")
    plt.legend(framealpha = 1.0)
    plt.title('Cyclogenesis Mean Latitude')
    plt.savefig('figures/mean_latitude', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_mean_lat(): Done")
plot_mean_lat(storms)<jupyter_output><empty_output><jupyter_text>Mean PressureSame as above, but for the mean cyclone lowest pressure for each year<jupyter_code>def plot_mean_pressure(stormObj):
    print("plot_mean_pressure()")
    fig = plt.figure(figsize=(12,9))

    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    cyclGroups = ["Clipper", "Colorado", "Northwest", "GreatBasin", "GulfOfMexico", "EastCoast"]

    for year in years:
        for type in cyclGroups:
            meanPressure = 0
            count = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == type and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] == year):
                    lowestP = np.nanmin(stormObj[ed]['amp']) / 100
                    meanPressure += lowestP
                    count += 1
            if(count != 0):
                mP = meanPressure / count
                plt.plot(year, mP, '-o', color=zoneColors[type])

    # Regression Lines
    for type in cyclGroups:
        PList = []
        for year in years:
            pres = 0
            meanPres = 0
            mP = 0
            count = 0
            for ed in range(len(stormObj)):
                if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['classification'] == type and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] == year):
                    lowestP = np.nanmin(stormObj[ed]['amp']) / 100
                    meanPres += lowestP
                    count += 1
            if count != 0:
                mP = meanPres / count
            else:
                mP = np.nan

            PList.append(mP)
            print(type + " (" + str(year) + "): " + str(mP))

        finiteYmask = np.isfinite(np.array(PList))
        Yclean = np.array(PList)[finiteYmask]
        Xclean = np.array(years)[finiteYmask]

        slope, intercept, r_value, p_value, std_err = stats.linregress(Xclean, Yclean)
        line = slope*years+intercept
        plt.plot(years, line, color=zoneColors[type], linestyle='--')

        print("Statistics " + type + ": m: %.3f, b: %.3f, r: %.3f, p: %.3f, stdErr: %.3f, stdDev: %.3f" % (slope, intercept, r_value, p_value, std_err, np.std(Yclean)))
    # Add "Fake" Lines for Legend
    for type in cyclGroups:	
        plt.plot(0, 0, linestyle='-', alpha=1, color=zoneColors[type], label=typeStr[type])
    plt.xlim(1979, 2020)
    plt.ylim(950, 1020)
    plt.xlabel("Year")
    plt.ylabel("Pressure (hPa)")
    plt.legend(framealpha = 1.0)
    plt.title('Cyclone Mean Lowest Pressure')
    plt.savefig('figures/mean_lowest_pressure.pdf', bbox_inches='tight', pad_inches=0.05, dpi=300)
    print("plot_mean_pressure(): Done")
plot_mean_pressure(storms)<jupyter_output><empty_output><jupyter_text>Lifespan Multi-PanelPanel 1: Lifespan of all cyclonesPanel 2: Lifespan of bomb cyclones<jupyter_code>def plot_lifespan_panel(stormObj):
    print("plot_lifespan_panel()")

    fig = plt.figure(figsize=(12,9))
    gs = fig.add_gridspec(1, 2)
    ax1 = fig.add_subplot(gs[0, 0])
    ax2 = fig.add_subplot(gs[0, 1])

    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    dt = 3

    longest = 0
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            count = len(stormObj[ed]['amp'])
            if(count > longest):
                longest = count
                
    print(longest * dt)

    # Panel A
    lifespan_array1 = np.zeros((longest))
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            lifespan = len(stormObj[ed]['amp'])
            lifespan_array1[lifespan-1] += 1

    idx_of_last_non_zero = np.max(np.nonzero(lifespan_array1))
    splitArray1 = lifespan_array1[0:idx_of_last_non_zero+1]
    times1 = np.arange(0, (len(splitArray1)*3), 3)
    
    print(splitArray1[8:])
    print(splitArray1.shape)
    print(times1[8:])
    print(times1.shape)

    #plt.bar(times, splitArray)
    ax1.scatter(times1[8:], splitArray1[8:], s=2)
    ax1.set(title="Cyclone Lifespan (All Cyclones)", xlabel='Hours', ylabel='Count')

    #ax1.

    # Panel B
    lifespan_array2 = np.zeros((longest))
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            bergeron = calculate_bergeron(stormObj[ed])
            if(any(b >= 1 for b in bergeron)):
                lifespan = len(stormObj[ed]['amp'])
                lifespan_array2[lifespan-1] += 1

    idx_of_last_non_zero = np.max(np.nonzero(lifespan_array2))
    splitArray2 = lifespan_array2[0:idx_of_last_non_zero+1]
    times2 = np.arange(0, (len(splitArray2)*3), 3)

    ax2.bar(times2, splitArray2)
    #ax1.scatter(times2, splitArray2)
    ax2.set(title="Cyclone Lifespan (Bomb Cyclones Only)", xlabel='Hours', ylabel='Count')

    plt.suptitle('Lifespan of Cyclones')
    plt.savefig('figures/lifespan_panel', bbox_inches='tight', pad_inches=0.05, dpi=300)	
    print("plot_lifespan_panel(): Done")
plot_lifespan_panel(storms)<jupyter_output><empty_output><jupyter_text>Extras Statistics, Value Outputs, Etc<jupyter_code>def getPctLifespan(stormObj):
    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    dt = 3

    longest = 0
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            count = len(stormObj[ed]['amp'])
            if(count > longest):
                longest = count
                
    lastTimeStep = longest*dt            
    
    print(lastTimeStep, lastTimeStep%24, (lastTimeStep + (lastTimeStep%24) %24))
    
    timeBins = np.arange(1, 16+1, 1)
    timeBins *= 24
    
    counts = np.zeros(timeBins.shape)
    
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            stormLife = len(stormObj[ed]['amp']) * dt
            
            closest = (np.abs(timeBins - stormLife)).argmin()
            if(stormLife > timeBins[closest]):
                closest += 1
                
            counts[closest] += 1
            
    #print(counts)
    
    maxIdx = 12
    totalStorms = 0
    for i in range(maxIdx+1):
        totalStorms += counts[i]
        
    pctOfWhole = (counts[0:maxIdx] / totalStorms) * 100
    print(timeBins[0:maxIdx] / 24)
    print(np.array_str(pctOfWhole, precision = 5, suppress_small=True))
    
getPctLifespan(storms)
def check_bombs(stormObj, cName="Other", onlyShowBomb = False):
    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2021, 1)    
    
    i = 0
    bCt = 0
    
    for ed in range(len(stormObj)):
        isBomb = False
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            if(stormObj[ed]['classification'] == cName):
                bergeron = calculate_bergeron(stormObj[ed])
                if(any(b >= 1 for b in bergeron)):
                    isBomb = True
                    bCt += 1
                if(onlyShowBomb):
                    if(isBomb):
                        print("Storm " + str(i+1) + "\n" + str(stormObj[ed]) + "\n" + str(bergeron) + "\n Bomb? " + str(isBomb))
                else:
                    print("Storm " + str(i+1) + "\n" + str(stormObj[ed]) + "\n" + str(bergeron) + "\n Bomb? " + str(isBomb))
                i += 1
    print("Total Bombs: " + str(bCt))
check_bombs(storms, "Colorado", False)
def make_bomb_lines(stormObj):
    print("make_bomb_lines()")
    
    clipper_zone = [[-115,50],
                    [-105,50],
                    [-105,55],
                    [-110,55],
                    [-110,60],
                    [-125,60],
                    [-125,55],
                    [-115,55]]
    northwst_zone = [[-125, 60],
                    [-115,60],
                    [-115,65],
                    [-125,65]]
    colorado_zone = [[-105, 40],
                    [-100,40],
                    [-100,35],
                    [-105,35]]
    basin_zone =    [[-120, 45],
                    [-115,45],
                    [-115,35],
                    [-120,35]]
    gom_zone =      [[-100, 25],
                    [-90,25],
                    [-90,30],
                    [-100,30]]
    eastcst_zone =  [[-80, 30],
                    [-75,30],
                    [-75,35],
                    [-65,35],
                    [-65,45],
                    [-70,45],
                    [-70,40],
                    [-80,40]]    
    
    
    plotExtent4P = [-125, -70, 23, 60] #plotExtent
    
    
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    
    clipperPoly = mplPath.Path(clipper_zone)
    northwestPoly = mplPath.Path(northwst_zone)
    coloradoPoly = mplPath.Path(colorado_zone)
    basinPoly = mplPath.Path(basin_zone)
    gomPoly = mplPath.Path(gom_zone)
    eastCoastPoly = mplPath.Path(eastcst_zone)

    patch1 = matplotlib.patches.PathPatch(clipperPoly, facecolor=zoneColors["Clipper"], transform=ccrs.PlateCarree(), label="Alberta Clipper")
    patch2 = matplotlib.patches.PathPatch(northwestPoly, facecolor=zoneColors["Northwest"], transform=ccrs.PlateCarree(), label="Northwestern")
    patch3 = matplotlib.patches.PathPatch(coloradoPoly, facecolor=zoneColors["Colorado"], transform=ccrs.PlateCarree(), label="Colorado Low")
    patch4 = matplotlib.patches.PathPatch(basinPoly, facecolor=zoneColors["GreatBasin"], transform=ccrs.PlateCarree(), label="Great Basin Low")
    patch5 = matplotlib.patches.PathPatch(gomPoly, facecolor=zoneColors["GulfOfMexico"], transform=ccrs.PlateCarree(), label="Gulf of Mexico Low")
    patch6 = matplotlib.patches.PathPatch(eastCoastPoly, facecolor=zoneColors["EastCoast"], transform=ccrs.PlateCarree(), label="East Coast Low")
    ax.add_patch(patch1)
    ax.add_patch(patch2)
    ax.add_patch(patch3)
    ax.add_patch(patch4)
    ax.add_patch(patch5)
    ax.add_patch(patch6)    
    
    monthList = [1, 2, 3, 4, 10, 11, 12]
    yearList = np.arange(1979, 2020, 1)
    cycloneList = ['Other']#['Clipper', 'Northwest', 'Colorado', 'GreatBasin', 'GulfOfMexico', 'EastCoast']    
       
    ax.set_extent(plotExtent4P, crs=ccrs.PlateCarree())   
    ax.add_feature(states, edgecolor='k', linewidth=0.5)    
    ax.coastlines()
    
    i = 0
    
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic'):
            if(stormObj[ed]['month'][0] in monthList):
                if(storms[ed]['classification'] in cycloneList):
                    i+=1
                    bergeron = calculate_bergeron(stormObj[ed])
                    if(any(b >= 1 for b in bergeron)):                 
                        lon, lat = stormObj[ed]['lon'], stormObj[ed]['lat']
                        lon[lon <  0] = 360 + lon[lon < 0]
                        ax.plot(lon, lat, color=zoneColors[stormObj[ed]['classification']], linewidth=0.2, alpha=0.35, transform=ccrs.PlateCarree())
                        ax.plot(lon[0], lat[0], 'bo', transform=ccrs.PlateCarree())
                        #ax.text(lon[0], lat[0]+1, str(i), transform=ccrs.PlateCarree())
make_bomb_lines(storms)
def count_before42(stormObj):
    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    dt = 3

    cycloneList = ['Clipper', 'Northwest', 'Colorado', 'GreatBasin', 'GulfOfMexico', 'EastCoast']  
    
    # Find the longest cyclone.
    longest = 0
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' 
            and stormObj[ed]['month'][0] in months 
            and stormObj[ed]['year'][0] in years
            and stormObj[ed]['classification'] in cycloneList):
            count = len(stormObj[ed]['amp'])
            if(count > longest):
                longest = count

    totalCount = 0
    before42 = 0
                
    #totals_column = np.zeros((longest))
    for className in cycloneList:
        fall_array = np.zeros((longest))
        for ed in range(len(stormObj)):
            if (stormObj[ed]['type'] == 'cyclonic' 
                and stormObj[ed]['month'][0] in months 
                and stormObj[ed]['year'][0] in years 
                and stormObj[ed]['classification'] == className):
                dPressure = np.diff(stormObj[ed]['amp'])
                largestFall = np.nanmin(dPressure)
                index_of_largest_fall = np.argmin(dPressure)
                time_of_fall = index_of_largest_fall * dt
                fall_array[index_of_largest_fall] += 1
                
                if(time_of_fall <= 42):
                    before42 += 1
                totalCount += 1

    print(f"Totals of storms ({totalCount}), and counted ({before42}) => {(before42/totalCount)*100}%")
count_before42(storms)
def longtrack(stormObj):
    months = [1, 2, 3, 4, 10, 11, 12]
    years = np.arange(1979, 2020, 1)
    dt = 3

    longest = 0
    idx = 0
    for ed in range(len(stormObj)):
        if (stormObj[ed]['type'] == 'cyclonic' and stormObj[ed]['month'][0] in months and stormObj[ed]['year'][0] in years):
            count = len(stormObj[ed]['amp'])
            if(count > longest):
                longest = count
                idx = ed
                
    print(f"Longest Storm {longest * dt}hr, Storm Index {idx}")
    
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    #ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent(plotExtent, crs=ccrs.PlateCarree())   

    # Draw our plot, coastlines first, then the contours.
    ax.add_feature(states, edgecolor='k')    
    ax.coastlines()
    
    lon, lat = stormObj[idx]['lon'], stormObj[idx]['lat']
    lon[lon <  0] = 360 + lon[lon < 0]
    plt.plot(lon, lat, 'r-', linewidth=0.5, alpha=0.35, transform=ccrs.PlateCarree()) 
    
longtrack(storms)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.141087069307535, 'fraction_numerical': 0.02999969481490524, 'mean_word_length': 2.5731421405594026, 'pattern_counts': {'":': 31, '<': 97, '<?xml version=': 0, '>': 109, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '404238', 'license': 'MIT', 'n_tokens_mistral': 23509, 'n_tokens_neox': 21740, 'n_words': 4319, 'path': 'Plotting.ipynb', 'repo_name': 'niu-wcs/Storm-Climatology-1979-2019', 'size': '65534'}