Where Hurricanes Hit Hardest: A County-Level Evaluation with Python

been retained by an insurance company to assist refine home insurance premiums across the southeastern United States. Their query is easy but high stakes: ? They usually don’t just mean , they need to account for storms that keep driving , delivering damaging rain and spawning tornadoes.

To tackle this, you’ll need two key ingredients:

• A reliable storm track database
• A county boundary shapefile

With these, the workflow becomes clear: discover and count every hurricane track that intersects a county boundary, then visualize the leads to each map and list form for max insight.

Python is an excellent fit for this job, due to its wealthy ecosystem of geospatial and scientific libraries:

• Tropycal for pulling open-source government hurricane data
• GeoPandas for loading and manipulating geospatial files
• Plotly Express for constructing interactive, explorable maps

Before diving into the code, let’s examine the outcomes. We’ll concentrate on the period 1975 to 2024, when global warming, to influence Atlantic hurricanes, became firmly established.

During the last 49 years, 338 hurricanes have struck 640 counties within the southeastern US. Coastal counties bear the brunt of wind and storm surge, while inland regions suffer from torrential rain and the occasional hurricane-spawned tornado. It’s a posh, far-reaching hazard, and with the appropriate tools, you possibly can map it county by county.

The next map, built using the Tropycal library, records the tracks of all of the hurricanes that made landfall within the US from 1975 through 2024.

While interesting, this map isn’t much use to an insurance adjuster. We want to quantify it by adding county-level resolution and counting the variety of tracks that cross into each county. Here’s how that appears:

Now we have now a greater idea of which counties act as “hurricane magnets.” Across the Southeast, hurricane “hit” counts range from zero to 12 per county — however the storms are removed from evenly distributed. Hotspots cluster along the Louisiana coast, in central Florida, and along the shorelines of the Carolinas. The East Coast really takes it on the chin, with Brunswick County, North Carolina, holding the unwelcome record for essentially the most hurricane strikes.

A look on the track map makes the pattern clear. Florida, Georgia, South Carolina, and North Carolina sit within the crosshairs of two storm highways — one from the Atlantic and one other from the Gulf of Mexico. The prevailing westerlies, which begin just north of the Gulf Coast, often bend northward-tracking storms toward the Atlantic seaboard. Fortunately for Georgia and the Carolinas, a lot of these systems lose strength over land, slipping below hurricane force before sweeping through.

For insurers, these visualizations aren’t just weather curiosities; they’re decision-making tools. And layering in historical loss data can provide a more complete picture of the true financial cost of living by the water’s edge.

The Choropleth Code

The next code, written in JupyterLab, creates a choropleth map of hurricane track counts per county. It uses geospatial data from the Plotly graphing library and pulls open-source weather data from the National Oceanic and Atmospheric Administration (NOAA) using the Tropycal library.

The code uses the next packages:

python 3.10.18
numpy 2.2.5
geopandas 1.0.1
plotly 6.0.1 (plotly_express 0.4.1)
tropical 1.4
shapely 2.0.6

Importing Libraries

Start by importing the next libraries.

import json
import numpy as np
import geopandas as gpd
import plotly.express as px
from tropycal import tracks
from shapely.geometry import LineString

Configuring Constants

Now, we arrange several constants. The primary is a set of the state “FIPS” codes. Short for , these “zip codes for states” are commonly utilized in geospatial files. On this case, they represent the southeastern states of Alabama, Florida, Georgia, Louisiana, Mississippi, North Carolina, South Carolina, and Texas. Later, we’ll use these codes to filter a single file of your complete USA.

# CONFIGURE CONSTANTS
# State: AL, FL, GA, LA, MS, NC, SC, TX:
SE_STATE_FIPS = {'01', '12', '13', '22', '28', '37', '45', '48'}  
YEAR_RANGE = (1975, 2024)
INTENSITY_THRESH = {'v_min': 64}  # Hurricanes (>= 64 kt)
COUNTY_GEOJSON_URL = (
 'https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json'
)

Next, we define a 12 months range (1975-2024) as a tuple. Then, we assign an intensity threshold constant for wind speed. Tropycal will filter storms based on wind speeds, and people with speeds of 64 knots or greater are classified as hurricanes.

Finally, we offer the URL address for the Plotly library’s counties geospatial shapefile. Later, we’ll use GeoPandas to load this as a GeoDataFrame, which is basically a pandas DataFrame with a Geometry column for geospatial mapping information.

NOTE: Hurricanes quickly develop into tropical storms and depressions after making landfall. These are still destructive, nonetheless, so we’ll proceed to trace them.

Defining Helper Functions

To streamline the hurricane mapping workflow, we’ll define three lightweight helper functions. These will help keep the code modular, readable, and adaptable, especially when working with real-world geospatial data which will vary in structure or scale.

# Define Helper Functions:
def get_hover_name_column(df: gpd.GeoDataFrame) -> str:
    # Prefer proper-case county name if available:
    if 'NAME' in df.columns:
        return 'NAME'
    if 'name' in df.columns:
        return 'name'
    # Fallback to id if no name column exists:
    return 'id'

def storm_to_linestring(storm_obj) -> LineString | None:
    df = storm_obj.to_dataframe()
    if len(df) < 2:
        return None
    coords = [(lon, lat) for lon, lat in zip(df['lon'], df['lat'])
              if not (np.isnan(lon) or np.isnan(lat))]
    return LineString(coords) if len(coords) > 1 else None

def make_tickvals(vmax: int) -> list[int]:
    if vmax <= 10:
        step = 2
    elif vmax <= 20:
        step = 4
    elif vmax <= 50:
        step = 10
    else:
        step = 20
    return list(range(0, int(vmax) + 1, step)) or [0]

Plotly Express creates visuals. We’ll have the ability to hover the cursor over counties within the choropleth map and launch a pop-up window of the county name and the variety of hurricanes which have passed through. The get_hover_name_column(df) function selects essentially the most readable column name for map hover labels. It checks for 'NAME' or 'name' within the GeoDataFrame and defaults to 'id' if neither is found. This ensures consistent labeling across datasets.

The storm_to_linestring(storm_obj) function converts a storm’s track data right into a LineString geometry by extracting valid longitude–latitude pairs. If the storm has fewer than two valid points, it returns ‘None’. This is important for spatial joins and visualizing storm paths.

Finally, the make_tickvals(vmax) function generates a clean set of tick marks for the choropleth colorbar based on the utmost hurricane count. It dynamically adjusts the step size to maintain the legend readable, whether the range is small or large.

Prepare the County Map

The following cell loads the geospatial data and filters out the southeastern states, using our prepared set of FIPS codes. In the method, it creates a GeoDataFrame and adds a column for the Plotly Express hover data.

# Load and filter county boundary data:
counties_gdf = gpd.read_file(COUNTY_GEOJSON_URL)

# Ensure FIPS id is string with leading zeros:
counties_gdf['id'] = counties_gdf['id'].astype(str).str.zfill(5)

# Derive state code from id's first two digits:
counties_gdf['STATE_FIPS'] = counties_gdf['id'].str[:2]
se_counties_gdf = (counties_gdf[counties_gdf['STATE_FIPS'].
                   isin(SE_STATE_FIPS)].copy())
hover_col = get_hover_name_column(se_counties_gdf)

print(f"Loading county data...")
print(f"Loaded {len(se_counties_gdf)} southeastern counties")

To start, we load the county-level GeoJSON file using GeoPandas and prepare it for evaluation. Each county is identified by a FIPS code, which we format as a 5-digit string to make sure consistency (the primary two digits represent the state code). We then extract the state portion of every FIPS code and filter the dataset to incorporate only counties in our eight southeastern states. Finally, we select a column for labeling counties within the hover text and make sure the variety of counties which have been loaded.

Fetching and Processing Hurricane Data

Now it’s time to make use of Tropycal to fetch and process the hurricane data from the . That is where we programmatically overlay the counties with the hurricane tracks and count the unique occurrences of tracks in each county.

# Get and process hurricane data using Tropycal library:
try:
    atlantic = tracks.TrackDataset(basin='north_atlantic', 
                                   source='hurdat', 
                                   include_btk=True)

    storms_ids = atlantic.filter_storms(thresh=INTENSITY_THRESH,
                                                  year_range=YEAR_RANGE)

    print(f"Found {len(storms_ids)} hurricanes from "
    f"{YEAR_RANGE[0]}–{YEAR_RANGE[1]}")

    storm_names = []
    storm_tracks = []

    for i, sid in enumerate(storms_ids, start=1):
        if i % 50 == 0 or i == 1 or i == len(storms_ids):
            print(f"Processing storm {i}/{len(storms_ids)}")
        try:
            storm = atlantic.get_storm(sid)
            geom = storm_to_linestring(storm)
            if geom shouldn't be None:
                storm_tracks.append(geom)
                storm_names.append(storm.name)
        except Exception as e:
            print(f"  Skipped {sid}: {e}")

    print(f"Successfully processed {len(storm_tracks)} storm tracks")

    hurricane_tracks_gdf = gpd.GeoDataFrame({'name': storm_names}, 
                                            geometry=storm_tracks,
                                            crs="EPSG:4326")

    # Pre-filter tracks to the bounding box of the SE counties for speed:
    xmin, ymin, xmax, ymax = se_counties_gdf.total_bounds
    hurricane_tracks_gdf = hurricane_tracks_gdf.cx[xmin:xmax, ymin:ymax]

    # Check that county data and hurricane tracks are same CRS:
    assert se_counties_gdf.crs == hurricane_tracks_gdf.crs, 
    f"CRS mismatch: {se_counties_gdf.crs} vs {hurricane_tracks_gdf.crs}"

    # Spatial join to search out counties intersecting hurricane tracks:
    print("Performing spatial join...")
    joined = gpd.sjoin(se_counties_gdf[['id', hover_col, 'geometry']],
                       hurricane_tracks_gdf[['name', 'geometry']],
                       how="inner",
                       predicate="intersects")

    # Count unique hurricanes per county:
    unique_pairs = joined[['id', 'name']].drop_duplicates()
    hurricane_counts = (unique_pairs.groupby('id', as_index=False).size().
                                    rename(columns={'size': 'hurricane_count'}))

    # Merge counts back
    se_counties_gdf = se_counties_gdf.merge(hurricane_counts, 
                                            on='id', 
                                            how='left')
    se_counties_gdf['hurricane_count'] = (se_counties_gdf['hurricane_count'].
                                          fillna(0).astype(int))

    print(f"Hurricane counts: Max: {se_counties_gdf['hurricane_count'].max()} | "
          f"Nonzero counties: {(se_counties_gdf['hurricane_count'] > 0).sum()}")

except Exception as e:
    print(f"Error loading hurricane data: {e}")
    print("Creating sample data for demonstration...")
    np.random.seed(42)
    se_counties_gdf['hurricane_count'] = np.random.poisson(2, 
                                                           len(se_counties_gdf))

Here’s a breakdown of the foremost steps:

• Load Dataset: Initializes the TrackDataset using HURDAT data, including best track (btk) points.
• Filter Storms: Selects hurricanes that meet a specified intensity threshold and fall inside a given 12 months range.
• Extract Tracks: Iterates through each storm ID, converts its path to a LineString geometry, and stores each the track and storm name. Progress is printed every 50 storms.
• Create GeoDataFrame: Combines storm names and geometries right into a GeoDataFrame with WGS84 coordinates.
• Spatial Filtering: Clips hurricane tracks to the bounding box of southeastern counties to enhance performance.
• Assert CRS: Checks that the county and hurricane data use the identical coordinate reference system (in case you desire to use different geospatial and/or hurricane track files).
• Spatial Join: Identifies which counties intersect with hurricane tracks using a spatial join.

Performing the spatial join could be tricky. For instance, if a track doubles back and re-enters a county, you don’t need to count it twice.

To handle this, the code first identifies unique name pairs after which drops duplicate rows from the GeoDataFrame before performing the count.

• Count Hurricanes per County:
  • Drops duplicate storm–county pairs.
  • Groups by county ID to count unique hurricanes.
  • Merges results back into the county GeoDataFrame.
  • Fills missing values with zero and converts to integer.
• Fallback Handling: If hurricane data fails to load, synthetic hurricane counts are generated using a Poisson distribution for demonstration purposes. That is for learning the method, only!

Errors loading the hurricane data are common, so keep watch over the printout. If the information fails to load, keep rerunning the cell until it does.

A successful run will yield the next confirmation:

Constructing the Choropleth Map

The following cell generates a customized choropleth map of hurricane counts per county within the Southeastern US using Plotly Express.

# Construct the choropleth map:
print("Creating choropleth map...")

se_geojson = json.loads(se_counties_gdf.to_json())
max_count = int(se_counties_gdf['hurricane_count'].max())
tickvals  = make_tickvals(max_count)

fig = px.choropleth(se_counties_gdf, 
                    geojson=json.loads(se_counties_gdf.to_json()),
                    locations='id',
                    featureidkey='properties.id',
                    color='hurricane_count',
                    color_continuous_scale='Reds',
                    range_color=[0, max_count],
                    title=(f"Southeastern US: Hurricane Counts Per County "
                           f"({YEAR_RANGE[0]}–{YEAR_RANGE[1]})"),
                    hover_name=hover_col,
                    hover_data={'hurricane_count': True, 'id': False})

# Adjust the map layout and clean the Plotly hover data:
fig.update_geos(fitbounds="locations", visible=False)
fig.update_traces(
    hovertemplate="%{hovertext}
Hurricanes: %{z}"
)

fig.update_layout(
    width=1400,
    height=1000,
    title=dict(
        text=(f"Southeastern US: Hurricane Counts Per County "
              f"({YEAR_RANGE[0]}–{YEAR_RANGE[1]})"),
        x=0.5,
        xanchor='center',
        y=0.85,
        yanchor='top',
        font=dict(size=24),
        pad=dict(t=0, b=10)
    ),
    coloraxis_colorbar=dict(
        x=0.96,
        y=0.5,
        len=0.4,
        thickness=16,
        title='Hurricane Count',
        outlinewidth=1,
        tickvals=tickvals,
        tickfont=dict(size=16)
    )
)

fig.add_annotation(
    text="Data: HURDAT2 via Tropycal | Metric: counties intersecting hurricane "
         f"tracks ({YEAR_RANGE[0]}–{YEAR_RANGE[1]})",
    x=0.521,
    y=0.89,
    showarrow=False,
    font=dict(size=16),
    xanchor='center'
)

fig.show()

The important thing steps here include:

• GeoJSON Conversion: Converts the GeoDataFrame of counties to GeoJSON format for simple mapping with Plotly Express.
• Color Scaling: Determines the utmost hurricane count and calls the helper function to create tick values for the colorbar.
• Map Rendering:
  • Uses px.choropleth to visualise hurricane_count per county.
    • The locations='id' argument tells Plotly which column within the GeoDataFrame accommodates the unique identifiers for every county (county-level FIPS codes). These values match each row of information to the corresponding shape within the GeoJSON file.
    • The featureidkey='properties.id' argument specifies where to search out the matching identifier contained in the GeoJSON structure. GeoJSON features have a properties dictionary containing an 'id' field. This ensures that every county’s geometry is accurately paired with its hurricane count.
    • Applies a red color scale, sets the range, and defines hover behavior.
• Layout & Styling:
  • Centers and styles the title.
  • Adjusts map bounds and hides geographic outlines.
    • The fig.update_geos(fitbounds="locations", visible=False) line turns off the bottom map for a cleaner plot.
  • Refines hover tooltips for clarity.
  • Customizes the colorbar with tick marks and labels.
• Annotation: Adds an information source note referencing HURDAT2 and the evaluation metric.
• Display: Shows the ultimate interactive map with fig.show().

The deciding think about using Plotly Express over static tools like Matplotlib is the addition of the dynamic hover data. Since there’s no practical strategy to label lots of of counties, the hover data permits you to query the map while keeping all that extra information out of sight until needed.

The Track Map Code

Although unnecessary, viewing the actual hurricane tracks can be a pleasant touch, in addition to a strategy to check the choropleth results. This map could be generated entirely with the Tropycal library, as shown below.

# Plot tracks coloured by category:
title = 'SE USA Hurricanes (1975-2024)'
ax = atlantic.plot_storms(storms=storms_ids,
                          title=title,
                          domain={'w':-97.68,'e':-70.3,'s':22,'n':ymax},
                          prop={'plot_names':False, 
                                'dots':False,
                                'linecolor':'category',
                                'linewidth':1.0}, 
                          map_prop={'plot_gridlines':False})

# plt.savefig('counties_tracks.png', dpi=600, bbox_inches='tight')

Note that the domain parameter refers back to the boundaries of the map. While you should use our previous xmin, xmax, ymin, and ymax variables, I’ve adjusted them barely for a more visually appealing map. Here’s the result:

For more on using the Tropycal library, see my previous article: Easy Hurricane Tracking with Tropycal | by Lee Vaughan | TDS Archive | Medium.

The Hurricane List Code

No insurance adjuster will need to cursor through a map to extract data. Because GeoDataFrames are a type of pandas DataFrame, it’s easy to slice and dice the information and present it as tables. The next code sorts the counties by hurricane count after which, for brevity, displays the highest 20 counties based on their count.

Here’s the short and straightforward strategy to generate this table; I’ve added some extra code for the state abbreviations:

# Map FIPS to state abbreviation:
fips_to_abbrev = {'01': 'AL', '12': 'FL', '13': 'GA', '22': 'LA', 
                  '28': 'MS', '37': 'NC', '45': 'SC', '48': 'TX'}

# Add state abbreviation column:
se_counties_gdf['state_abbrev'] = se_counties_gdf['STATE'].map(fips_to_abbrev)

# Sort and choose top 20 counties by hurricane count
top20 = (se_counties_gdf.sort_values(by='hurricane_count', 
                                     ascending=False) 
         [['state_abbrev', 'NAME', 'hurricane_count']].head(20))

# Display result
print(top20.to_string(index=False))

And here’s the result:

While this works, it’s not very professional-looking. We will improve it using an HTML approach:

# Print out the highest 20 counties based on hurricane impacts:
# Map FIPS to state abbreviation:
fips_to_abbrev = {'01': 'AL', '12': 'FL', '13': 'GA', '22': 'LA', 
                  '28': 'MS', '37': 'NC', '45': 'SC', '48': 'TX'}

gdf_sorted = se_counties_gdf.copy()

# Add latest column for state abbreviation:
gdf_sorted['State Name'] = gdf_sorted['STATE'].map(fips_to_abbrev)

# Rename Existing Columns:
# Multiple columns without delay
gdf_sorted = gdf_sorted.rename(columns={'NAME': 'County Name', 
                                        'hurricane_count': 'Hurricane Count'})

# Sort by hurricane_count:
gdf_sorted = gdf_sorted.sort_values(by='Hurricane Count', ascending=False)

# Create a horny HTML display:
df_display = gdf_sorted[['State Name', 'County Name', 'Hurricane Count']].head(20)
df_display['Hurricane Count'] = df_display['Hurricane Count'].astype(int)

# Create styled HTML table without index:
styled_table = (
    df_display
        .style
        .set_caption("Top 20 Counties by Hurricane Impacts")
        .set_table_styles([
            # Hide the index
            {'selector': 'th.row_heading',
             'props': [('display', 'none')]},
            {'selector': 'th.blank',
             'props': [('display', 'none')]},

            # Caption styling:
            {'selector': 'caption',
             'props': [('caption-side', 'top'),
                       ('font-size', '16px'),
                       ('font-weight', 'bold'),
                       ('text-align', 'center'),
                       ('color', '#333')]},

            # Header styling:
            {'selector': 'th.col_heading',
             'props': [('background-color', '#004466'),
                       ('color', 'white'),
                       ('text-align', 'center'),
                       ('padding', '6px')]},

            # Cell styling:
            {'selector': 'td',
             'props': [('text-align', 'center'),
                       ('padding', '6px')]}
        ])
        # Add zebra striping:
        .apply(lambda col: [
            'background-color: #f2f2f2' if i % 2 == 0 else ''
            for i in range(len(col))
        ], axis=0)
)

# Save styled HTML table to disk:
styled_table.to_html("top20_hurricane_table.html")
styled_table

This cell transforms our raw geospatial data right into a clean, publication-ready summary of hurricane exposure by county. It prepares and presents a ranked table of the 20 counties most affected by hurricanes:

• State Abbreviation Mapping: It starts by mapping each county’s FIPS state code to its two-letter abbreviation (e.g., '48' → 'TX') and adds this as a brand new column.
• Column Renaming: The county name ('NAME') and hurricane count ('hurricane_count') columns are renamed to 'County Name' and 'Hurricane Count' for clarity.
• Sorting and Selection: The GeoDataFrame is sorted in descending order by hurricane count, and the highest 20 rows are chosen.
• Styled Table Creation: Using pandas’ Styler, the code builds a visually formatted HTML table:
  • Adds a centered caption
  • Hides the index column
  • Applies custom header and cell styling
  • Adds zebra striping for readability
• Export to HTML: The styled table is saved as top20_hurricane_table.html, making it easy to embed in reports or share externally.

Here’s the result:

This table could be further enhanced by including interactive sorting or by embedding it directly right into a dashboard.

Summary

On this project, we addressed a matter on every actuary’s desk: Python’s wealthy ecosystem of third-party packages was key to creating this easy and effective. Tropycal made accessing government hurricane data a breeze, Plotly provided the county boundaries, and GeoPandas merged the 2 datasets and counted the variety of hurricanes per county. Finally, Plotly Express produced a dynamic, interactive map that made it easy to visualise and explore the county-level hurricane data.