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

been retained by an insurance company to help refine home insurance premiums across the southeastern United States. Their question is simple but high stakes: which counties are hit most often by hurricanes? And they don’t just mean landfall, they want to account for storms that keep driving inland, 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: identify and count every hurricane track that intersects a county boundary, then visualize the results in both map and list form for maximum insight.

Python is an ideal fit for this job, thanks to its rich ecosystem of geospatial and scientific libraries:

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

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

Over the last 49 years, 338 hurricanes have struck 640 counties in 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 complex, far-reaching hazard, and with the right tools, you can map it county by county.

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

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

Now we have a better idea of which counties act as “hurricane magnets.” Across the Southeast, hurricane “hit” counts range from zero to twelve per county — but the storms are far 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 the most hurricane strikes.

A glance at the track map makes the pattern clear. Florida, Georgia, South Carolina, and North Carolina sit in the crosshairs of two storm highways — one from the Atlantic and another 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, many 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 following 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 following 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 following 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 set up several constants. The first is a set of the state “FIPS” codes. Short for Federal Information Processing Series, these “zip codes for states” are commonly used in geospatial files. In 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 the entire 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 year 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 those with speeds of 64 knots or greater are classified as hurricanes.

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

NOTE: Hurricanes quickly become tropical storms and depressions after making landfall. These are still destructive, however, so we’ll continue to track 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 that may 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 interactive visuals. We’ll be able to hover the cursor over counties in the choropleth map and launch a pop-up window of the county name and the number of hurricanes that have passed through. The get_hover_name_column(df) function selects the most readable column name for map hover labels. It checks for 'NAME' or 'name' in 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 into a LineString geometry by extracting valid longitude–latitude pairs. If the storm has fewer than two valid points, it returns ‘None’. This is essential 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 maximum hurricane count. It dynamically adjusts the step size to keep the legend readable, whether the range is small or large.

Prepare the County Map

The next cell loads the geospatial data and filters out the southeastern states, using our prepared set of FIPS codes. In the process, 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 begin, we load the county-level GeoJSON file using GeoPandas and prepare it for analysis. Each county is identified by a FIPS code, which we format as a 5-digit string to ensure consistency (the first two digits represent the state code). We then extract the state portion of each FIPS code and filter the dataset to include only counties in our eight southeastern states. Finally, we select a column for labeling counties in the hover text and confirm the number of counties that have been loaded.

Fetching and Processing Hurricane Data

Now it’s time to use Tropycal to fetch and process the hurricane data from the National Hurricane Center. This 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 is not 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 find 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 major 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 within a given year range.
• Extract Tracks: Iterates through each storm ID, converts its path to a LineString geometry, and stores both the track and storm name. Progress is printed every 50 storms.
• Create GeoDataFrame: Combines storm names and geometries into a GeoDataFrame with WGS84 coordinates.
• Spatial Filtering: Clips hurricane tracks to the bounding box of southeastern counties to improve performance.
• Assert CRS: Checks that the county and hurricane data use the same coordinate reference system (in case you want 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 can be tricky. For example, if a track doubles back and re-enters a county, you don’t want to count it twice.

To handle this, the code first identifies unique name pairs and then 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. This is for learning the process, only!

Errors loading the hurricane data are common, so keep an eye on the printout. If the data fails to load, keep rerunning the cell until it does.

A successful run will yield the following confirmation:

Building the Choropleth Map

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

# Build 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 key steps here include:

• GeoJSON Conversion: Converts the GeoDataFrame of counties to GeoJSON format for easy mapping with Plotly Express.
• Color Scaling: Determines the maximum hurricane count and calls the helper function to create tick values for the colorbar.
• Map Rendering:
  • Uses px.choropleth to visualize hurricane_count per county.
    • The locations='id' argument tells Plotly which column in the GeoDataFrame contains the unique identifiers for each county (county-level FIPS codes). These values match each row of data to the corresponding shape in the GeoJSON file.
    • The featureidkey='properties.id' argument specifies where to find the matching identifier inside the GeoJSON structure. GeoJSON features have a properties dictionary containing an 'id' field. This ensures that each county’s geometry is correctly 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 base map for a cleaner plot.
  • Refines hover tooltips for clarity.
  • Customizes the colorbar with tick marks and labels.
• Annotation: Adds a data source note referencing HURDAT2 and the analysis metric.
• Display: Shows the final interactive map with fig.show().

The deciding factor in using Plotly Express over static tools like Matplotlib is the addition of the dynamic hover data. Since there’s no practical way to label hundreds of counties, the hover data lets you 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 would be a nice touch, as well as a way to check the choropleth results. This map can be generated entirely with the Tropycal library, as shown below.

# Plot tracks colored 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 to the boundaries of the map. While you can use our previous xmin, xmax, ymin, and ymax variables, I’ve adjusted them slightly 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 want to cursor through a map to extract data. Because GeoDataFrames are a form of pandas DataFrame, it’s easy to slice and dice the data and present it as tables. The following code sorts the counties by hurricane count and then, for brevity, displays the top 20 counties based on their count.

Here’s the quick and easy way 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 select 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 can improve it using an HTML approach:

# Print out the top 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 new column for state abbreviation:
gdf_sorted['State Name'] = gdf_sorted['STATE'].map(fips_to_abbrev)

# Rename Existing Columns:
# Multiple columns at once
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 an attractive 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 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 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 top 20 rows are selected.
• 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 can be further enhanced by including interactive sorting or by embedding it directly into a dashboard.

Summary

In this project, we addressed a question on every actuary’s desk: Which counties get hit the hardest by hurricanes, year after year? Python’s rich ecosystem of third-party packages was key to making this easy and effective. Tropycal made accessing government hurricane data a breeze, Plotly provided the county boundaries, and GeoPandas merged the two datasets and counted the number of hurricanes per county. Finally, Plotly Express produced a dynamic, interactive map that made it easy to visualize and explore the county-level hurricane data.