<jupyter_start><jupyter_text>Storm Census Part of our modified GIT repository, adjusted from the original package to work with Python 3, and be more optimized compared to the original versionRequired Modules:* NumPY* Matplotlib* Cartopy* SciPyRequired File:* storm_track_slp.npz<jupyter_code>import numpy as np
from matplotlib import pyplot as plt
import cartopy
import cartopy.crs as ccrs
from shapely.geometry import MultiLineString
from shapely.geometry import LineString
from shapely.geometry import Point
from shapely.geometry import Polygon
from shapely.ops import polygonize
import geopandas as gpd<jupyter_output><empty_output><jupyter_text>Data LoadRun the below code block once to load the data<jupyter_code>data = np.load('storm_track_slp.npz', encoding='latin1', allow_pickle=True)
storms = data['storms']<jupyter_output><empty_output><jupyter_text>Program SettingsThese are the settings for the census code, you can adjust the extent of the census domain, the count of the grid cells, and adjust a flag for debug prints<jupyter_code># testMode: Set this to True to run the code blocks as debug-only (Will only print results, and not save to a file)
testMode = False

# plotExtent: The domain of the census box (Format is: Lower Lon, Upper Lon, Lower Lat, Upper Lat)
plotExtent = [-160, -40, 20, 70]

# gridResolution: The factor of resolution for the lon x lat box, the count is currently set to:
#  (upperlon - lowerlon) * gridResolution & (upperlat - lowerlat) * gridResolution
gridResolution = 1

# timeRange: The start and end years of the analyzed data
timeRange = [1979, 2020]<jupyter_output><empty_output><jupyter_text>Helper FunctionsRun this code block prior to running the census below, it contains helper functions used by the census code.<jupyter_code># A quick and simple point-in-box code
def storm_in_grid(grid, stormLat, stormLon):
    gridLon_l = grid[0]
    gridLon_u = grid[1]
    gridLat_l = grid[2]
    gridLat_u = grid[3]

    return ((stormLat >= gridLat_l and stormLat <= gridLat_u) and (stormLon >= gridLon_l and stormLon <= gridLon_u))

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 = sum(val for val in diffs if val < 0)

        bFactor = (-1 * 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>SetupRun this code block after applying the settings to store global variables which are needed (If you alter the settings at any point, rerun this block)<jupyter_code>tRange = np.arange(timeRange[0], timeRange[1]+1, 1)

polyIn = gpd.read_file("aea_grid_climo.shp")
polyFrame = polyIn.to_crs('epsg:4326') 

polyFrame["ID"] = polyFrame.index

# Populate the dataframe with the years
for y in tRange:
    polyFrame["count_" + str(y)] = 0
    polyFrame["count_bomb_" + str(y)] = 0
    
print(polyFrame)<jupyter_output><empty_output><jupyter_text>Grid CheckRun this block of code to test the grid to make sure it looks good.<jupyter_code>projObj = get_projection_object()    
# Create our figure and axis object
fig = plt.figure(figsize=(12,9))
ax = plt.axes(projection=projObj)
ax.set_extent(plotExtent, crs=ccrs.PlateCarree())

states = cartopy.feature.NaturalEarthFeature(category='cultural',
                                            name='admin_1_states_provinces_lakes',
                                            scale='110m',
                                            facecolor='none')
ax.add_feature(states, edgecolor='k')    
ax.coastlines()

for(index, row) in polyFrame.iterrows():
    ax.add_geometries([row.geometry], crs=ccrs.PlateCarree(), edgecolor='black', alpha=0.8)<jupyter_output><empty_output><jupyter_text>Census CodeRun the below code to calculate the storm census (This will take some time depending on your computer's CPU / RAM).Currently, the code bundles all storms together, however since we're loading the storm_track_slp.npz file you can easily access objects such as *storms[i]['classification']* to split by other fields in the storm object if desired.<jupyter_code># Regardless of test mode or not, we need to start by resetting the counts.
for y in tRange:
    polyFrame["count_" + str(y)] = 0
    polyFrame["count_bomb_" + str(y)] = 0

if testMode:
    print("Debug mode enabled")

    storm_lats = storms[0]['lat']
    storm_lons = storms[0]['lon']
    
    storm_lons[storm_lons > 180] = storm_lons[storm_lons > 180] - 360

    polyLine = LineString(list(zip(storm_lons, storm_lats)))
    coord_tests = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polyLine, crs='EPSG:4326')) 
    
    I = gpd.sjoin(polyFrame, coord_tests)
    
    I["count_" + str(int(storms[0]['year'][0]))] += 1
    
    polyFrame.update(I)
    
    projObj = get_projection_object()    
    # Create our figure and axis object
    fig = plt.figure(figsize=(12,9))
    ax = plt.axes(projection=projObj)
    ax.set_extent(plotExtent, crs=ccrs.PlateCarree())
   
    states = cartopy.feature.NaturalEarthFeature(category='cultural',
                                                name='admin_1_states_provinces_lakes',
                                                scale='110m',
                                                facecolor='none')
    ax.add_feature(states, edgecolor='k')    
    ax.coastlines()
    
    lx, ly = polyLine.xy
    ax.plot(lx, ly, transform=ccrs.PlateCarree(), color='green')
    
    #for(index, row) in polyFrame.iterrows()
    #ax.add_geometries(polyFrame.geometry, crs=ccrs.PlateCarree(), edgecolor='black', alpha=0.8)
    
    for index, row in I.iterrows():
        p = row[0]
        
        ax.add_geometries([p], crs=ccrs.PlateCarree(), facecolor='b', edgecolor='red', alpha=0.8)

else:
    print("Beginning census counting")

    newN = 0    
    newNBomb = 0

    for ed in range(len(storms)):#range(1537, len(storms)):
        skip = False
        bomb = False
        bergeron = calculate_bergeron(storms[ed])
        if(any(b >= 1 for b in bergeron)):
            bomb = True
            
        coordList = []
    
        cYear = storms[ed]['year'][0]
        outYear = storms[ed]['year'][-1]
        for t in range(len(storms[ed]['year'])):
            storm_lat = storms[ed]['lat'][t]
            storm_lon = storms[ed]['lon'][t]
            
            if storm_lon > 180:
                storm_lon = storm_lon - 360
            
            pair = [storm_lon, storm_lat]
            
            storm_year = storms[ed]['year'][t]
            if(storm_year != cYear):
                # The year on this track has changed, break off the line segment here
                cLen = len(coordList)
                if(cLen > 0 and cLen <= 1):
                    singlePoint = Point(coordList[0])
                    coord_tests = gpd.GeoDataFrame(geometry=gpd.GeoSeries(singlePoint, crs='epsg:4326'))  
                    I = gpd.sjoin(polyFrame, coord_tests)
                else:                    
                    polyLine = LineString(coordList)
                    coord_tests = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polyLine, crs='epsg:4326'))
                    I = gpd.sjoin(polyFrame, coord_tests)
                    
                # storm_year is the current year, so we subtract one to go back to where these data occurred
                I["count_" + str(int(storm_year-1))] += 1
                newN += len(I.index)
                if(bomb):
                    I["count_bomb_" + str(int(storm_year-1))] += 1
                    newNBomb += len(I.index)

                polyFrame.update(I)
                # Clear the list
                coordList = []
                # Update the current year
                cYear = storm_year
            else:
                coordList.append(pair)
        # End of list, if the year has not changed, we're now ready to process.
        cLen = len(coordList)
        if(cLen == 0):
            # This shouldn't happen, but it's here to prevent an error...
            skip = True
        elif(cLen > 0 and cLen <= 1):
            # In the event the list is one point (Termination point on 00Z Jan-1)
            singlePoint = Point(coordList[0])
            coord_tests = gpd.GeoDataFrame(geometry=gpd.GeoSeries(singlePoint, crs='epsg:4326')) 
            I = gpd.sjoin(polyFrame, coord_tests)
        else:            
            polyLine = LineString(coordList)
            coord_tests = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polyLine, crs='epsg:4326')) 
            I = gpd.sjoin(polyFrame, coord_tests)
                
        if not skip:
            I["count_" + str(int(outYear))] += 1
            newN += len(I.index)
            if(bomb):
                I["count_bomb_" + str(int(outYear))] += 1
                newNBomb += len(I.index)

            polyFrame.update(I)            
        if(ed % 100 == 0):
            print("Completed " + str(ed) + " / " + str(len(storms)) + " (" + str((ed/len(storms)) * 100) + "%) -> Added " + str(newN) + " N Points " +
                  "| Bergeron -> Added " + str(newNBomb) + " N Points")
            newN = 0
            newNBomb = 0          

    print("Done processing")<jupyter_output><empty_output><jupyter_text>SavingRun the below block of code once the above processing is complete to save the output to a shape file.<jupyter_code>print(polyFrame)

print("Data Saved")
polyFrame.to_file('storm_census.geojson', driver='GeoJSON')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10322047894302229, 'fraction_numerical': 0.015506009725662904, 'mean_word_length': 2.3767038413878563, 'pattern_counts': {'":': 2, '<': 38, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '368227', 'license': 'MIT', 'n_tokens_mistral': 3498, 'n_tokens_neox': 3277, 'n_words': 860, 'path': 'Storm Census.ipynb', 'repo_name': 'niu-wcs/Storm-Climatology-1979-2019', 'size': '10899'}