import pandas as pd	# pandas is a powerful data manipulation and analysis library for structured data.
import numpy as np	# numpy is a comprehensive library for numerical computations and working with arrays.
from datetime import datetime, timedelta	# datetime provides classes for manipulating dates and times.
from pathlib import Path	# pathlib offers an object-oriented approach for handling filesystem paths

import geopandas as gpd	# geopandas extends pandas to allow spatial operations on geometric types.
from shapely.geometry import Point, LineString, box	# shapely.geometry provides geometric objects such as points, lines, and rectangles for spatial analysis.
from shapely.ops import nearest_points	# shapely.ops.nearest_points finds the nearest geometric points between objects.
import fiona	# fiona is a library for reading and writing vector data files in GIS formats.

import matplotlib.pyplot as plt	# matplotlib.pyplot is the main plotting library for creating static visualizations in Python.
import matplotlib.dates as mdates	# matplotlib.dates provides functions to format and plot dates on matplotlib graphs.
from matplotlib.patches import Patch	# matplotlib.patches.Patch is used to create simple shapes for legend entries in plots.
from matplotlib.lines import Line2D	# matplotlib.lines.Line2D is used for drawing and customizing lines in plots.
import folium	# folium creates interactive web-based maps for geospatial data visualization.
from folium import plugins	# folium.plugins offers additional plugins for enhanced interactive map features.

from sklearn.ensemble import IsolationForest	# sklearn.ensemble.IsolationForest is an unsupervised anomaly detection algorithm.
from sklearn.preprocessing import StandardScaler	# sklearn.preprocessing.StandardScaler standardizes numerical features by removing the mean and scaling to unit variance.

from joblib import Parallel, delayed	# joblib.Parallel and delayed enable easy-and-efficient parallel computation.
import multiprocessing	# multiprocessing allows for spawning processes for parallel execution in Python.

import warnings	# warnings provides a way to manage and filter Python warning messages.
warnings.filterwarnings('ignore')	# Disables the display of warning messages globally.

plt.style.use('seaborn-v0_8-whitegrid')	# Sets the plot style to 'seaborn-v0_8-whitegrid' for improved aesthetics in matplotlib.
plt.rcParams['figure.figsize'] = (12, 8)	# Sets default output figure size for matplotlib visualizations.
plt.rcParams['font.size'] = 11	# Configures the default font size for matplotlib plots.

print("All libraries imported successfully!")

All libraries imported successfully!

# Configuration: File Paths and Parameters
import os	# import the os module for operating system interactions
BASE_DIR = Path.cwd()

# Data paths
AIS_GDB_PATH = BASE_DIR / "GulfOfAmerica.gdb"	# path to the Gulf of America geodatabase file
print(AIS_GDB_PATH)
HARVEY_SHAPEFILES_DIR = BASE_DIR / "Harvey track points, lines, radii, windswath (shapefiles)"	# directory containing Hurricane Harvey shapefiles

# Hurricane Harvey study period (UTC timezone for AIS data)
HARVEY_START = pd.Timestamp("2017-08-15", tz='UTC')	# study period start date (UTC)
HARVEY_END = pd.Timestamp("2017-09-05", tz='UTC')	# study period end date (UTC)
HARVEY_LANDFALL = pd.Timestamp("2017-08-26 03:00:00", tz='UTC')	# Harvey first landfall datetime (UTC)
HARVEY_SECOND_LANDFALL = pd.Timestamp("2017-08-30 00:00:00", tz='UTC')	# Harvey second landfall datetime (UTC)
HARVEY_IMPACT_START = pd.Timestamp("2017-08-25 00:00:00", tz='UTC')	# beginning of main impact period (UTC)
HARVEY_IMPACT_END = HARVEY_SECOND_LANDFALL  # end of main impact period (use second landfall)
HARVEY_RECOVERY_START = HARVEY_IMPACT_END + pd.Timedelta(days=1)    # start of recovery period (1 day after second landfall)

# Study area bounding boxes
HOUSTON_BBOX = {'minx': -95.6, 'maxx': -94.8, 'miny': 29.5, 'maxy': 29.9}	# bounding box for Houston area study region
CORPUS_BBOX = {'minx': -97.6, 'maxx': -96.8, 'miny': 27.5, 'maxy': 28.0}	# bounding box for Corpus Christi area
TEXAS_COAST_BBOX = {'minx': -98.0, 'maxx': -93.5, 'miny': 26.0, 'maxy': 30.5}	# bounding box for entire Texas coast

print("Configuration loaded successfully!")	# notify configuration loaded without errors
print(f"Study period: {HARVEY_START.date()} to {HARVEY_END.date()}")	# print the dates of the study period

/scratch/user/u.jz324199/2026_cybertraining_copy/Maritime Traffic/GulfOfAmerica.gdb
Configuration loaded successfully!
Study period: 2017-08-15 to 2017-09-05

layers = fiona.listlayers(AIS_GDB_PATH)

# List all layers in the geodatabase	layers = fiona.listlayers(AIS_GDB_PATH)	# get list of all layer names in the AIS geodatabase
print(f"Found {len(layers)} layers in the AIS geodatabase:")
for i, layer in enumerate(layers, 1):	# loop through layers
    print(f"  {i}. {layer}")	# print each layer's index and name

Found 12 layers in the AIS geodatabase:
  1. Tracks_2017_03
  2. Tracks_2017_04
  3. Tracks_2017_05
  4. Tracks_2017_06
  5. Tracks_2017_07
  6. Tracks_2017_08
  7. Tracks_2017_09
  8. Tracks_2017_10
  9. Tracks_2017_11
  10. Tracks_2017_12
  11. Tracks_2017_01
  12. Tracks_2017_02

print("Loading AIS data for August and September 2017...")	# print message indicating which data is loading

ais_august = gpd.read_file(AIS_GDB_PATH, layer='Tracks_2017_08')	# read August 2017 AIS track data
ais_september = gpd.read_file(AIS_GDB_PATH, layer='Tracks_2017_09')	# read September 2017 AIS track data

print(f"August 2017: {len(ais_august):,} vessel track segments")	# print number of vessel track segments in August
print(f"September 2017: {len(ais_september):,} vessel track segments")	# print number of vessel track segments in September

Loading AIS data for August and September 2017...
August 2017: 277,473 vessel track segments
September 2017: 251,496 vessel track segments

ais_data = pd.concat([ais_august, ais_september], ignore_index=True)	# combine August and September data into one DataFrame

# After pd.concat, the result is often: Plain pandas DataFrame, because pandas doesn’t preserve spatial metadata reliably.
ais_data = gpd.GeoDataFrame(ais_data, geometry='geometry', crs=ais_august.crs)	# convert combined data to GeoDataFrame with appropriate geometry and crs

#Raw AIS timestamps are often Strings, so we need to create a new parsed datetime column
ais_data['datetime'] = pd.to_datetime(ais_data['TrackStartTime'])	# create new datetime column from TrackStartTime

#print out the feature name

print(f"\nFeature name:\n  {ais_data.columns} \n")

#print out the first 3 rows:
print(f"\nThe first three row:\n {ais_data.head(3)}")


print(f"\nTotal records: {len(ais_data):,}")	# print total number of combined records
print(f"Date range: {ais_data['datetime'].min()} to {ais_data['datetime'].max()}")	# print overall date range in combined data
print(f"CRS: {ais_data.crs}")	# print the coordinate reference system
print(f"\nVessel types: {ais_data.VesselGroup.unique()}")	# print unique vessel types in the data set

Feature name:
  Index(['MMSI', 'TrackStartTime', 'TrackEndTime', 'VesselType', 'Length',
       'Width', 'Draft', 'DurationMinutes', 'VesselGroup', 'Shape_Length',
       'geometry', 'datetime'],
      dtype='object')


The first three row:
         MMSI            TrackStartTime              TrackEndTime  VesselType  \
0  209004000 2017-08-01 00:00:08+00:00 2017-08-03 16:27:30+00:00      1004.0
1  209004000 2017-08-03 16:57:33+00:00 2017-08-04 01:12:29+00:00      1004.0
2  209004000 2017-08-04 01:42:29+00:00 2017-08-04 03:12:28+00:00      1004.0

   Length  Width  Draft  DurationMinutes VesselGroup  Shape_Length  \
0   229.0  32.26   14.5             3867       Cargo      0.000844
1   229.0  32.26   14.5              494       Cargo      0.000000
2   229.0  32.26   14.5               89       Cargo      0.000000

                                            geometry                  datetime
0  MULTILINESTRING ((-90.78767 30.00750, -90.7878... 2017-08-01 00:00:08+00:00
1                                               None 2017-08-03 16:57:33+00:00
2                                               None 2017-08-04 01:42:29+00:00

Total records: 528,969
Date range: 2017-08-01 00:00:00+00:00 to 2017-09-30 23:58:42+00:00
CRS: EPSG:4269

Vessel types: ['Cargo' 'Tanker' 'Others' 'Passenger' 'TugTow' None
 'Pleasure Craft/Sailing' 'Fishing' 'Not Available']

# Filter to blue economy vessels (Cargo, Tanker, Passenger, Fishing)
all_data = ais_data
ais_data = ais_data[ais_data.VesselGroup.isin(['Cargo', 'Tanker', 'Passenger', 'Fishing'])]

print(f"Blue economy vessels: {len(ais_data):,} tracks")
print(f"Percentage of total: {len(ais_data)/len(all_data)*100:.1f}%")

Blue economy vessels: 275,205 tracks
Percentage of total: 52.0%

# Load Hurricane Harvey shapefiles
# below is equivalent to: HARVEY_SHAPEFILES_DIR / "al092017_pts" / "al092017_pts.shp"
harvey_pts_path = os.path.join(HARVEY_SHAPEFILES_DIR, "al092017_pts", "al092017_pts.shp")	# Construct path to hurricane track points shapefile
harvey_lin_path = os.path.join(HARVEY_SHAPEFILES_DIR, "al092017_lin", "al092017_lin.shp")	# Construct path to hurricane track line shapefile
harvey_radii_path = os.path.join(HARVEY_SHAPEFILES_DIR, "al092017_radii", "al092017_radii.shp")	# Construct path to hurricane wind radii shapefile

print("Loading Hurricane Harvey shapefiles...")
harvey_points = gpd.read_file(harvey_pts_path, engine='fiona')	# Load hurricane track points into GeoDataFrame
harvey_line = gpd.read_file(harvey_lin_path, engine='fiona')	# Load hurricane track line into GeoDataFrame
harvey_radii = gpd.read_file(harvey_radii_path, engine='fiona')	# Load hurricane wind radii into GeoDataFrame

print(f"  Track points: {len(harvey_points)} records")	# Print the number of records in track points
print(f"  Track line: {len(harvey_line)} records")	# Print the number of records in track line
print(f"  Wind radii: {len(harvey_radii)} records")	# Print the number of records in wind radii

Loading Hurricane Harvey shapefiles...
  Track points: 74 records
  Track line: 17 records
  Wind radii: 61 records

# Process hurricane track data
# Example output: Timestamp('2017-08-26 03:00:00+0000', tz='UTC')
def parse_dtg(dtg):	# Define function to parse DTG (date-time group) format to datetime
    """Parse DTG format (e.g., 2017082603) to datetime."""
    dtg_str = str(dtg)	# Convert DTG to string
    year = int(dtg_str[:4])	# Extract year from DTG string
    month = int(dtg_str[4:6])	# Extract month from DTG string
    day = int(dtg_str[6:8])	# Extract day from DTG string
    hour = int(dtg_str[8:10]) if len(dtg_str) >= 10 else 0	# Extract hour from DTG string if available, otherwise set to 0
    return pd.Timestamp(year=year, month=month, day=day, hour=hour, tz='UTC')	# Return pandas Timestamp object with UTC timezone


# How apply works:
# for value in column:
#     parse_dtg(value)
harvey_points['datetime'] = harvey_points['DTG'].apply(parse_dtg)	# Apply parse_dtg to DTG column to create datetime column


# Rename columns for consistency
harvey_points = harvey_points.rename(columns={	# Rename columns in harvey_points DataFrame for consistency
    'LAT': 'lat', 'LON': 'lon', 'INTENSITY': 'wind_kt',	# Rename LAT to lat, LON to lon, INTENSITY to wind_kt
    'MSLP': 'pressure_mb', 'STORMTYPE': 'status', 'SS': 'saffir_simpson'	# Rename MSLP to pressure_mb, STORMTYPE to status, SS to saffir_simpson
})

# Before:
    # DTG = 2017082603
    # INTENSITY = 120
    # MSLP = 950
# After
    # datetime = Timestamp
    # wind_kt = 120
    # pressure_mb = 950

# Mark landfall
harvey_points['is_landfall'] = False	# Initialize 'is_landfall' column as False
#.loc is the Pandas label-based indexing.
harvey_points.loc[harvey_points['datetime'] == pd.Timestamp('2017-08-26 03:00:00', tz='UTC'), 'is_landfall'] = True	# Set 'is_landfall' to True at landfall time

# Set CRS
harvey_points = harvey_points.set_crs('EPSG:4326', allow_override=True)	# Set coordinate reference system to EPSG:4326 for points
harvey_line = harvey_line.set_crs('EPSG:4326', allow_override=True)	# Set coordinate reference system to EPSG:4326 for line

harvey_hurdat = harvey_points.copy()	# Make a copy of harvey_points as harvey_hurdat

print(f"\nHurricane track processed:")	# Print message that hurricane track has been processed
print(f"  Date range: {harvey_hurdat['datetime'].min()} to {harvey_hurdat['datetime'].max()}")	# Print the date range in hurricane data
print(f"  Peak intensity: {harvey_hurdat['wind_kt'].max():.0f} kt (Category {harvey_hurdat['saffir_simpson'].max()})")	# Print the peak wind intensity and category

Hurricane track processed:
  Date range: 2017-08-16 06:00:00+00:00 to 2017-09-02 12:00:00+00:00
  Peak intensity: 115 kt (Category 4)

# Create a bounding box for the Texas coast region
texas_bbox = box(TEXAS_COAST_BBOX['minx'], TEXAS_COAST_BBOX['miny'],	TEXAS_COAST_BBOX['maxx'], TEXAS_COAST_BBOX['maxy'])	# Create a bounding box geometry from specified coordinates

print("Filtering AIS data to Texas coast region...")	# Print message indicating start of filtering
print(f"Bounding box: {TEXAS_COAST_BBOX}")	# Print the bounding box coordinates being used

# Remove rows with null geometry
ais_data_valid = ais_data[~ais_data.geometry.isna()].copy()	# Filter out rows where geometry is null and copy the result
print(f"Removed {len(ais_data) - len(ais_data_valid):,} rows with null geometry")	# Print the number of rows removed due to null geometry

# Convert CRS to WGS84
ais_data_valid = ais_data_valid.to_crs('EPSG:4326')	# Convert the coordinate reference system to WGS84

# Spatial filter
intersects_mask = ais_data_valid.geometry.intersects(texas_bbox)	# Create a boolean mask for rows where geometry intersects the Texas bbox
ais_texas = ais_data_valid.loc[intersects_mask].copy()	# Select and copy only the rows that intersect the bounding box

print(f"\nFiltered from {len(ais_data_valid):,} to {len(ais_texas):,} vessel tracks")	# Print before/after counts for filtering
print(f"Date range: {ais_texas['datetime'].min()} to {ais_texas['datetime'].max()}")	# Print the date range of the resulting data

Filtering AIS data to Texas coast region...
Bounding box: {'minx': -98.0, 'maxx': -93.5, 'miny': 26.0, 'maxy': 30.5}
Removed 25,760 rows with null geometry

Filtered from 249,445 to 79,224 vessel tracks
Date range: 2017-08-01 00:00:00+00:00 to 2017-09-30 23:58:26+00:00

ais_texas['date'] = ais_texas['datetime'].dt.normalize()	# Create a new column 'date' by normalizing 'datetime' to midnight (preserves UTC)

daily_stats = ais_texas.groupby('date').agg({	# Group ais_texas by date and calculate aggregate metrics
    'MMSI': 'nunique',  # Number of unique vessels
    'geometry': 'count',  # Number of track segments
}).rename(columns={'MMSI': 'unique_vessels', 'geometry': 'track_segments'})	# Rename columns for clarity

daily_stats['avg_segments_per_vessel'] = daily_stats['track_segments'] / daily_stats['unique_vessels']	# Calculate average segments per vessel

print("Daily Vessel Traffic Metrics:")	# Print description for metrics
print(f"Unique vessels in study area: {len(ais_texas.MMSI.unique()):,}")	# Print the count of unique vessels

print(f"\nTraffic around Hurricane Harvey landfall ({HARVEY_LANDFALL.date()}):")	# Print landfall traffic analytics range
landfall_window = daily_stats[	# Filter data for window surrounding Hurricane Harvey landfall
    (daily_stats.index >= pd.Timestamp('2017-08-20', tz='UTC')) & (daily_stats.index <= HARVEY_END)
]
display(landfall_window)	# Show the filtered dataset

Daily Vessel Traffic Metrics:
Unique vessels in study area: 1,994

Traffic around Hurricane Harvey landfall (2017-08-26):

# Create bounding boxes for specific ports
houston_bbox = box(HOUSTON_BBOX['minx'], HOUSTON_BBOX['miny'],	#Houston bounding box from coordinates
                   HOUSTON_BBOX['maxx'], HOUSTON_BBOX['maxy'])
corpus_bbox = box(CORPUS_BBOX['minx'], CORPUS_BBOX['miny'],	#Corpus Christi bounding box from coordinates
                  CORPUS_BBOX['maxx'], CORPUS_BBOX['maxy'])

# Filter AIS data to each port region
ais_houston = ais_texas[ais_texas.geometry.intersects(houston_bbox)].copy()	# Filter tracks that intersect Houston bounding box
ais_corpus = ais_texas[ais_texas.geometry.intersects(corpus_bbox)].copy()	# Filter tracks that intersect Corpus Christi bounding box

print("Port Region Filtering Results:")
print(f"Houston Ship Channel: {len(ais_houston):,} track segments")
print(f"Corpus Christi: {len(ais_corpus):,} track segments")

# Calculate daily statistics for each port
def get_daily_port_stats(gdf, port_name):
    gdf = gdf.copy()	# Make a copy of the input GeoDataFrame
    gdf['date'] = gdf['datetime'].dt.normalize()	# Add a 'date' column normalized to midnight
    stats = gdf.groupby('date').agg({'MMSI': 'nunique', 'geometry': 'count'})	# Group by date and count vessels and segments
    stats = stats.rename(columns={'MMSI': f'{port_name}_vessels', 'geometry': f'{port_name}_segments'})	# Rename columns with port name
    stats.index = pd.DatetimeIndex(stats.index)	# Ensure index is a DatetimeIndex
    return stats	# Return the daily stats

houston_daily = get_daily_port_stats(ais_houston, 'houston')	# Get daily stats for Houston port region
corpus_daily = get_daily_port_stats(ais_corpus, 'corpus')	# Get daily stats for Corpus Christi port region
port_stats = houston_daily.join(corpus_daily, how='outer').fillna(0)	# Join Houston and Corpus stats, fill missing with 0

print("\nPort traffic change (unique vessels, daily avg):")
for port_name, col in [("Houston", "houston_vessels"), ("Corpus Christi", "corpus_vessels")]:	# Loop over each port and its corresponding column name
    pre = port_stats[port_stats.index < HARVEY_IMPACT_START][col].mean()	# Calculate the mean vessel count before the impact start for the given port
    during = port_stats[(port_stats.index >= HARVEY_IMPACT_START) & (port_stats.index <= HARVEY_IMPACT_END)][col].mean()	# Calculate the mean vessel count during the impact window for the given port
    post = port_stats[port_stats.index > HARVEY_IMPACT_END][col].mean()	# Calculate the mean vessel count after the impact end for the given port
    drop_pct = ((pre - during) / pre * 100) if pre > 0 else 0	# Compute the percent drop in vessel count from pre-storm to during-storm, avoiding division by zero
    print(f"  {port_name}: pre {pre:.1f}, during {during:.1f}, post {post:.1f}, drop {drop_pct:.1f}%")	# Output the pre, during, post storm means and percentage drop for the port

Port Region Filtering Results:
Houston Ship Channel: 2,574 track segments
Corpus Christi: 2,317 track segments

Port traffic change (unique vessels, daily avg):
  Houston: pre 32.6, during 1.3, post 29.2, drop 95.9%
  Corpus Christi: pre 31.5, during 11.3, post 25.3, drop 64.0%

GALVESTON_BBOX = {'minx': -95.0, 'maxx': -94.6, 'miny': 29.2, 'maxy': 29.5}
galveston_bbox = box(GALVESTON_BBOX['minx'], GALVESTON_BBOX['miny'],
                     GALVESTON_BBOX['maxx'], GALVESTON_BBOX['maxy'])

# Step 1: Define Galveston bounding box and filter AIS data
GALVESTON_BBOX = {'minx': -95.0, 'maxx': -94.6, 'miny': 29.2, 'maxy': 29.5}
galveston_bbox = box(GALVESTON_BBOX['minx'], GALVESTON_BBOX['miny'],
                     GALVESTON_BBOX['maxx'], GALVESTON_BBOX['maxy'])

ais_galveston = ais_texas[ais_texas.geometry.intersects(galveston_bbox)].copy()

# Step 2: Calculate daily statistics
galveston_daily = get_daily_port_stats(ais_galveston, 'galveston')

# Step 3: Calculate pre-storm and during-storm averages
pre_storm = galveston_daily[galveston_daily.index < HARVEY_IMPACT_START]['galveston_vessels'].mean()
during_storm = galveston_daily[
    (galveston_daily.index >= HARVEY_IMPACT_START) &
    (galveston_daily.index <= HARVEY_IMPACT_END)
]['galveston_vessels'].mean()

# Calculate traffic reduction
reduction_pct = ((pre_storm - during_storm) / pre_storm * 100)

print(f"Galveston Port Traffic Analysis:")
print(f"Pre-storm average: {pre_storm:.1f} vessels/day")
print(f"During storm average: {during_storm:.1f} vessels/day")
print(f"Traffic reduction: {reduction_pct:.1f}%")

Galveston Port Traffic Analysis:
Pre-storm average: 54.8 vessels/day
During storm average: 9.0 vessels/day
Traffic reduction: 83.6%

# Time series plot of unique vessels with Harvey timeline
fig, ax = plt.subplots(figsize=(14, 6))	# Create a matplotlib figure and axis with set size

# Filter to Harvey study period
plot_data = daily_stats[	# Filter daily_stats DataFrame to include only dates during the Harvey study period
    (daily_stats.index >= HARVEY_START) & (daily_stats.index <= HARVEY_END)
]

# Plot unique vessels
ax.plot(plot_data.index, plot_data['unique_vessels'],	# Plot the unique_vessels time series on the axis
        color='blue', linewidth=2, marker='o', markersize=4)	# Set plot line/marker styles for the plot
ax.fill_between(plot_data.index, plot_data['unique_vessels'], alpha=0.3, color='blue')	# Fill area under the unique_vessels line with color

# Add Harvey timeline markers
ax.axvline(HARVEY_LANDFALL, color='red', linestyle='--', linewidth=2,	# Draw a vertical line for Harvey's first landfall
           label='First Landfall (Aug 26)')
ax.axvline(HARVEY_SECOND_LANDFALL, color='orange', linestyle='--',	# Draw a vertical line for Harvey's second landfall
           linewidth=2, label='Second Landfall (Aug 30)')

# Add phase annotations
ax.axvspan(HARVEY_START, HARVEY_LANDFALL - pd.Timedelta(days=2),	# Shade the pre-storm period on the plot
           alpha=0.1, color='green', label='Pre-Storm')
ax.axvspan(HARVEY_LANDFALL - pd.Timedelta(days=2), HARVEY_IMPACT_END,	# Shade the storm impact period on the plot
           alpha=0.2, color='red', label='Storm Impact')
ax.axvspan(HARVEY_RECOVERY_START, HARVEY_END,	# Shade the recovery period on the plot
           alpha=0.1, color='blue', label='Recovery')

ax.set_ylabel('Unique Vessels', fontsize=12)	# Set the y-axis label
ax.set_title('Daily Vessel Traffic in Texas Coast Region During Hurricane Harvey',	# Set the plot title
             fontsize=14, fontweight='bold')
ax.legend(loc='upper right', fontsize=9)	# Show legend in the upper right corner
ax.grid(True, alpha=0.3)	# Add a grid to the plot with transparency
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %d'))	# Format x-axis ticks as month and day
ax.xaxis.set_major_locator(mdates.DayLocator(interval=2))	# Set major ticks to every two days
plt.xticks(rotation=45)	# Rotate x-axis tick labels for readability
plt.tight_layout()	# Adjust layout for better spacing
plt.show()	# Show the plot

# Calculate traffic reduction
pre_storm = plot_data[plot_data.index < HARVEY_LANDFALL]['unique_vessels'].mean()	# Compute average unique vessels before landfall
during_storm = plot_data[	# Compute average unique vessels during the storm
    (plot_data.index >= HARVEY_LANDFALL.normalize()) &
    (plot_data.index <= HARVEY_IMPACT_END)
]['unique_vessels'].mean()

print(f"\nTraffic Analysis:")	# Print header for analysis output
print(f"  Average vessels before storm: {pre_storm:.1f}")	# Print average before storm
print(f"  Average vessels during storm: {during_storm:.1f}")	# Print average during storm
print(f"  Traffic reduction: {((pre_storm - during_storm) / pre_storm * 100):.1f}%")	# Print calculated traffic reduction percentage

Traffic Analysis:
  Average vessels before storm: 244.8
  Average vessels during storm: 138.6
  Traffic reduction: 43.4%

fig, axes = plt.subplots(2, 1, figsize=(14, 10), sharex=True)	# Create a figure with 2 subplots (one for each port), share x-axis

port_plot = port_stats[	# Filter port_stats DataFrame to rows between HARVEY_START and HARVEY_END
    (port_stats.index >= HARVEY_START) & (port_stats.index <= HARVEY_END)
]

ax1 = axes[0]	# Select the first subplot axis for Houston
ax1.plot(port_plot.index, port_plot['houston_vessels'], 		# Plot Houston vessel count time series as line with markers
         color='blue', linewidth=2, marker='o', markersize=5, label='Houston')
ax1.fill_between(port_plot.index, port_plot['houston_vessels'], alpha=0.3, color='blue')	# Fill under Houston curve with light blue
ax1.axvline(HARVEY_LANDFALL, color='red', linestyle='--', linewidth=2, alpha=0.7)	# Draw Harvey landfall vertical line
ax1.axvspan(HARVEY_IMPACT_START, HARVEY_IMPACT_END, 		# Highlight the storm impact period for Houston
            alpha=0.2, color='red', label='Peak Impact Period')
ax1.set_ylabel('Unique Vessels', fontsize=12)	# Label y-axis for Houston plot
ax1.set_title('Port of Houston - Vessel Traffic During Hurricane Harvey', 		# Set Houston plot title
              fontsize=13, fontweight='bold')
ax1.legend(loc='upper right')	# Display legend for Houston plot
ax1.grid(True, alpha=0.3)	# Add a light grid to Houston plot

ax2 = axes[1]	# Select the second subplot axis for Corpus Christi
ax2.plot(port_plot.index, port_plot['corpus_vessels'], 		# Plot Corpus Christi vessel count time series
         color='green', linewidth=2, marker='s', markersize=5, label='Corpus Christi')
ax2.fill_between(port_plot.index, port_plot['corpus_vessels'], alpha=0.3, color='green')	# Fill under Corpus curve with light green
ax2.axvline(HARVEY_LANDFALL, color='red', linestyle='--', linewidth=2, alpha=0.7,	# Draw Harvey landfall line for Corpus
            label='Landfall (Aug 26)')
ax2.axvspan(HARVEY_IMPACT_START, HARVEY_IMPACT_END, alpha=0.2, color='red')	# Highlight storm impact period on Corpus plot
ax2.set_xlabel('Date', fontsize=12)	# Label x-axis for Corpus plot
ax2.set_ylabel('Unique Vessels', fontsize=12)	# Label y-axis for Corpus plot
ax2.set_title('Port of Corpus Christi - Vessel Traffic During Hurricane Harvey', 		# Set Corpus plot title
              fontsize=13, fontweight='bold')
ax2.legend(loc='upper right')	# Display legend for Corpus plot
ax2.grid(True, alpha=0.3)	# Add a light grid to Corpus plot
ax2.xaxis.set_major_formatter(mdates.DateFormatter('%b %d'))	# Format date ticks as month and day
ax2.xaxis.set_major_locator(mdates.DayLocator(interval=2))	# Set major x-ticks to every two days
plt.xticks(rotation=45)	# Rotate x-tick labels for readability
plt.tight_layout()	# Adjust subplot spacing to prevent overlap
plt.show()	# Show the figure

print("\nPort Impact Comparison:")	# Print header for port comparison
print("=" * 60)	# Print a separator line
for port, col in [('Houston', 'houston_vessels'), ('Corpus Christi', 'corpus_vessels')]:	# Iterate over each port and its data column
    pre = port_plot[port_plot.index < HARVEY_IMPACT_START][col].mean()	# Compute pre-storm average for the port
    during = port_plot[	# Compute storm-period average for the port
        (port_plot.index >= HARVEY_IMPACT_START) &
        (port_plot.index <= HARVEY_IMPACT_END)
    ][col].mean()
    post = port_plot[port_plot.index > HARVEY_IMPACT_END][col].mean()	# Compute post-storm average for the port

    print(f"\n{port}:")	# Print port name
    print(f"  Pre-storm average: {pre:.1f} vessels/day")	# Print pre-storm average vessels
    print(f"  During storm: {during:.1f} vessels/day")	# Print during-storm average vessels
    print(f"  Post-storm: {post:.1f} vessels/day")	# Print post-storm average vessels
    if pre > 0:	# Only calculate reduction if pre-storm is valid
        print(f"  Traffic reduction during storm: {((pre - during) / pre * 100):.1f}%")	# Print percent reduction

Port Impact Comparison:
============================================================

Houston:
  Pre-storm average: 21.3 vessels/day
  During storm: 1.3 vessels/day
  Post-storm: 23.5 vessels/day
  Traffic reduction during storm: 93.7%

Corpus Christi:
  Pre-storm average: 24.0 vessels/day
  During storm: 11.3 vessels/day
  Post-storm: 29.7 vessels/day
  Traffic reduction during storm: 52.8%

# Interactive storm track and vessel movements (Folium)
center_lat = (TEXAS_COAST_BBOX['miny'] + TEXAS_COAST_BBOX['maxy']) / 2	# Calculate center latitude of study area
center_lon = (TEXAS_COAST_BBOX['minx'] + TEXAS_COAST_BBOX['maxx']) / 2	# Calculate center longitude of study area

storm_map = folium.Map(location=[center_lat, center_lon], zoom_start=6, tiles='cartodbpositron')	# Create Folium map centered on study area

# Study area outline
folium.Rectangle(
    bounds=[
        [TEXAS_COAST_BBOX['miny'], TEXAS_COAST_BBOX['minx']],	# Set southwest corner of rectangle
        [TEXAS_COAST_BBOX['maxy'], TEXAS_COAST_BBOX['maxx']],	# Set northeast corner of rectangle
    ],
    color='black', weight=2, fill=False, tooltip='Study area'	# Style rectangle appearance and add tooltip
).add_to(storm_map)	# Add rectangle to the map

# Storm track line
harvey_line_wgs84 = harvey_line.to_crs('EPSG:4326')	# Convert Harvey track line to WGS84 coordinate system
folium.GeoJson(
    harvey_line_wgs84,	# Add Harvey track line as GeoJSON layer
    name='Harvey track',	# Set layer name
    style_function=lambda _: {'color': 'darkred', 'weight': 3, 'opacity': 0.9}	# Style GeoJSON line appearance
).add_to(storm_map)	# Add GeoJSON to the map

# Storm advisory points
harvey_pts_wgs84 = harvey_points.to_crs('EPSG:4326')	# Convert advisory points to WGS84
for _, row in harvey_pts_wgs84.iterrows():	# Iterate over advisory points rows
    status = row['status'] if 'status' in harvey_pts_wgs84.columns else ''	# Get status if available, else empty
    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],	# Set marker position
        radius=4,	# Set marker size
        color='red',	# Set marker color
        fill=True,	# Fill marker
        fill_opacity=0.8,	# Set fill opacity
        popup=f"{row['datetime']:%Y-%m-%d %HZ} | {status}"	# Set popup text
    ).add_to(storm_map)	# Add marker to the map

# Prepare vessel track subsets for specific days (UTC)
landfall_day = HARVEY_LANDFALL.normalize()+ pd.Timedelta(days=2)
pre_day = landfall_day - pd.Timedelta(days=2)
post_day = landfall_day + pd.Timedelta(days=5)

pre_mask = ais_texas['datetime'].dt.normalize() == pre_day
landfall_mask = ais_texas['datetime'].dt.normalize() == landfall_day
post_mask = ais_texas['datetime'].dt.normalize() == post_day

def line_to_latlon(geom):	# Define function to extract [lat, lon] pairs from LineString or MultiLineString
    """Return list of [lat, lon] pairs from LineString or MultiLineString."""
    if geom is None or geom.is_empty:	# If geometry is empty
        return []	# Return empty list
    if geom.geom_type == 'LineString':	# If geometry is LineString
        return [[lat, lon] for lon, lat in geom.coords]	# Swap lon,lat to lat,lon for each coordinate pair
    if geom.geom_type == 'MultiLineString':	# If geometry is MultiLineString
        coords = []	# Initialize coords list
        for part in geom.geoms:	# For each part in MultiLineString
            coords.extend([[lat, lon] for lon, lat in part.coords])	# Add each lat,lon pair to coords
        return coords	# Return full list of lat,lon pairs
    return []	# For other geometry types, return empty list

pre_tracks = ais_texas.loc[pre_mask]
landfall_tracks = ais_texas.loc[landfall_mask]
post_tracks = ais_texas.loc[post_mask]

for gdf, color, label in [
    (pre_tracks, '#1E88E5', f"Pre-storm ({pre_day.date()})"),
    (landfall_tracks, '#F44336', f"Landfall day ({landfall_day.date()})"),
    (post_tracks, '#4CAF50', f"Post-storm ({post_day.date()})"),
]:
    layer = folium.FeatureGroup(name=label)
    for geom in gdf.geometry:	# Loop through each geometry in the phase
        coords = line_to_latlon(geom)	# Get list of lat,lon pairs for the geometry
        if len(coords) >= 2:	# Only plot lines with at least two coordinate pairs
            folium.PolyLine(coords, color=color, weight=2, opacity=0.7).add_to(layer)	# Add track to layer as polyline
    if len(layer._children) > 0:
        layer.add_to(storm_map)

folium.LayerControl().add_to(storm_map)

# Print traffic counts for the three focus days
for day_label, mask in [
    (f"Pre-storm {pre_day.date()}", pre_mask),
    (f"Landfall {landfall_day.date()}", landfall_mask),
    (f"Post-storm {post_day.date()}", post_mask),
]:
    subset = ais_texas.loc[mask]
    print(f"{day_label}: {subset['MMSI'].nunique()} unique vessels, {len(subset)} track segments")

storm_map

Pre-storm 2017-08-26: 90 unique vessels, 361 track segments
Landfall 2017-08-28: 115 unique vessels, 422 track segments
Post-storm 2017-09-02: 225 unique vessels, 824 track segments

def create_anomaly_features(daily_df):	# Define a function to create anomaly detection features from daily data
    """
    Create features for anomaly detection from daily vessel statistics.
    
    Features include:
    - Raw counts (vessels, segments)
    - Rolling statistics (3-day and 7-day moving averages)
    - Day-over-day changes
    - Deviation from rolling mean
    """
    features = pd.DataFrame(index=daily_df.index)	# Initialize an empty DataFrame with the same index as input data

    features['vessel_count'] = daily_df['unique_vessels']	# Add vessel count feature from unique vessel column
    features['segment_count'] = daily_df['track_segments']	# Add segment count feature from track segments column
    features['segments_per_vessel'] = daily_df['avg_segments_per_vessel']	# Add average segments per vessel feature

    # Moving averages smooth noise so “normal” behavior is clearer.
    features['vessel_ma3'] = daily_df['unique_vessels'].rolling(window=3, min_periods=1).mean()	# Compute 3-day moving average of vessel count
    features['vessel_ma7'] = daily_df['unique_vessels'].rolling(window=7, min_periods=1).mean()	# Compute 7-day moving average of vessel count

    features['vessel_change'] = daily_df['unique_vessels'].diff()	# Calculate day-over-day change in vessel count
    features['vessel_pct_change'] = daily_df['unique_vessels'].pct_change()	# Calculate day-over-day percent change in vessel count

    # Deviation from rolling mean, 0 = normal, negative = lower than normal, positive = higher than normal
    features['vessel_deviation'] = (	# Calculate deviation of vessel count from 7-day moving average
        daily_df['unique_vessels'] - features['vessel_ma7']	# Subtract 7-day moving average from vessel count
    ) / features['vessel_ma7'].replace(0, 1)	# Divide by 7-day moving average, replacing zero with one to avoid division by zero

    features['day_of_week'] = daily_df.index.dayofweek	# Add day of week as a feature

    return features.fillna(0)	# Return features DataFrame with NaN values filled with zero

anomaly_features = create_anomaly_features(daily_stats)	# Call the function to create the anomaly feature matrix

print("Anomaly Detection Feature Matrix:")	# Print informational string for feature matrix
print(f"Shape: {anomaly_features.shape}")	# Print the shape of the feature matrix
print(f"Features: {list(anomaly_features.columns)}")	# Print the list of feature column names
display(anomaly_features.head(10))	# Display the first 10 rows of the feature matrix

Anomaly Detection Feature Matrix:
Shape: (61, 9)
Features: ['vessel_count', 'segment_count', 'segments_per_vessel', 'vessel_ma3', 'vessel_ma7', 'vessel_change', 'vessel_pct_change', 'vessel_deviation', 'day_of_week']

# Prepare data for Isolation Forest
feature_cols = ['vessel_count', 'segment_count', 'vessel_change',
                'vessel_pct_change', 'vessel_deviation', 'vessel_ma3']	# Columns to use as features for anomaly detection

X = anomaly_features[feature_cols].copy()	# Select feature columns from anomaly_features and make a copy

scaler = StandardScaler()	# Create a StandardScaler object for feature standardization

#make each column have roughly a mean of 0 and a standard deviation of 1
X_scaled = scaler.fit_transform(X)	# Fit the scaler to the data and transform the features

iso_forest = IsolationForest(
    n_estimators=50000,	# Set the number of trees in the forest
    contamination=0.15,  # Expect ~15% of days to be anomalous
    random_state=42,	# Set a random seed for reproducibility
    max_samples='auto',	# Use all available samples for each tree
    n_jobs=8
)	# Create an IsolationForest object with specified parameters

iso_forest.fit(X_scaled)	# Fit the Isolation Forest model to the standardized feature data

anomaly_labels = iso_forest.predict(X_scaled)	# Predict anomaly labels for the data (-1 for anomaly, 1 for normal)
anomaly_scores = iso_forest.decision_function(X_scaled)	# Get anomaly scores (lower means more anomalous)

anomaly_features['is_anomaly'] = anomaly_labels == -1	# Add a boolean column to indicate if a day is anomalous
anomaly_features['anomaly_score'] = anomaly_scores	# Add the anomaly score to the dataframe

print("Isolation Forest Anomaly Detection Results:")	# Print header for the anomaly detection results display
print("=" * 60)	# Print a separator line
print(f"Total days analyzed: {len(anomaly_features)}")	# Print how many rows (days) were analyzed
print(f"Anomalies detected: {anomaly_features['is_anomaly'].sum()}")	# Print the number of detected anomalous days
print(f"Anomaly rate: {anomaly_features['is_anomaly'].mean()*100:.1f}%")	# Print the percentage of anomalies

print("\nDetected Anomaly Days:")	# Print label for the list of detected anomaly days
anomaly_days = anomaly_features[anomaly_features['is_anomaly']].sort_values('anomaly_score')	# Filter and sort detected anomaly days by anomaly score
display(anomaly_days[['vessel_count', 'vessel_change', 'anomaly_score']])	# Display vessel count, vessel change, and anomaly score for anomaly days

# Save model
# import joblib
# joblib.dump(iso_forest, 'iso_forest_model.pkl')
# print("\nModel saved to iso_forest_model.pkl")

Isolation Forest Anomaly Detection Results:
============================================================
Total days analyzed: 61
Anomalies detected: 9
Anomaly rate: 14.8%

Detected Anomaly Days:

# Visualize anomaly detection results
fig, axes = plt.subplots(3, 1, figsize=(14, 12), sharex=True)

plot_anomaly = anomaly_features[
    (anomaly_features.index >= HARVEY_START) &
    (anomaly_features.index <= HARVEY_END)
]

# Plot 1: Vessel count with anomalies highlighted
ax1 = axes[0]
ax1.plot(plot_anomaly.index, plot_anomaly['vessel_count'],
         color='#1E88E5', linewidth=2, label='Vessel Count')
anomaly_mask = plot_anomaly['is_anomaly']
ax1.scatter(plot_anomaly.index[anomaly_mask],
            plot_anomaly.loc[anomaly_mask, 'vessel_count'],
            color='red', s=80, zorder=5, label='Detected Anomaly', marker='X')
ax1.axvline(HARVEY_LANDFALL, color='red', linestyle='--', linewidth=2, alpha=0.7)
ax1.axvspan(HARVEY_IMPACT_START, HARVEY_IMPACT_END,
            alpha=0.15, color='red', label='Hurricane Impact Period')
ax1.set_ylabel('Vessel Count', fontsize=11)
ax1.set_title('GeoAI Anomaly Detection: Vessel Traffic During Hurricane Harvey',
              fontsize=13, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)

# Plot 2: Anomaly scores
ax2 = axes[1]
colors = ['red' if is_anom else '#1E88E5' for is_anom in plot_anomaly['is_anomaly']]
ax2.bar(plot_anomaly.index, plot_anomaly['anomaly_score'], color=colors, alpha=0.7)
ax2.axhline(0, color='black', linestyle='-', linewidth=0.5)
if anomaly_mask.any():
    threshold = plot_anomaly[plot_anomaly['is_anomaly']]['anomaly_score'].max()
    ax2.axhline(threshold, color='red', linestyle='--', linewidth=2, label='Anomaly Threshold')
ax2.axvline(HARVEY_LANDFALL, color='red', linestyle='--', linewidth=2, alpha=0.7)
ax2.set_ylabel('Anomaly Score\n(lower = more anomalous)', fontsize=11)
ax2.legend(loc='upper right')
ax2.grid(True, alpha=0.3)

# Plot 3: Day-over-day vessel change
ax3 = axes[2]
colors = ['red' if is_anom else '#4CAF50' for is_anom in plot_anomaly['is_anomaly']]
ax3.bar(plot_anomaly.index, plot_anomaly['vessel_change'], color=colors, alpha=0.7)
ax3.axhline(0, color='black', linestyle='-', linewidth=1)
ax3.axvline(HARVEY_LANDFALL, color='red', linestyle='--', linewidth=2, alpha=0.7, label='Landfall')
ax3.set_xlabel('Date', fontsize=11)
ax3.set_ylabel('Day-over-Day\nVessel Change', fontsize=11)
ax3.legend(loc='upper right')
ax3.grid(True, alpha=0.3)
ax3.xaxis.set_major_formatter(mdates.DateFormatter('%b %d'))
ax3.xaxis.set_major_locator(mdates.DayLocator(interval=2))
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# Evaluate model performance
print("\nModel Evaluation:")
print("=" * 60)
harvey_impact = (plot_anomaly.index >= HARVEY_IMPACT_START) & \
                (plot_anomaly.index <= HARVEY_IMPACT_END)
detected_during_impact = plot_anomaly[harvey_impact]['is_anomaly'].sum()
total_impact_days = harvey_impact.sum()
total_anomalies = plot_anomaly['is_anomaly'].sum()

print(f"Hurricane impact period: Aug 25-30 ({total_impact_days} days)")
print(f"Anomalies detected during impact: {detected_during_impact}")
print(f"Detection rate during impact: {detected_during_impact/max(total_impact_days, 1)*100:.1f}%")
print(f"\nTotal anomalies in study period: {total_anomalies}")
print(f"Precision: {detected_during_impact/max(total_anomalies, 1)*100:.1f}%")

Model Evaluation:
============================================================
Hurricane impact period: Aug 25-30 (6 days)
Anomalies detected during impact: 4
Detection rate during impact: 66.7%

Total anomalies in study period: 5
Precision: 80.0%

Date	Event
Aug 17, 2017	Formation as tropical depression
Aug 24, 2017	Intensified to hurricane status
Aug 25-26, 2017	Texas landfall near Rockport (Category 4)
Aug 28-30, 2017	Looped back over Gulf, second landfall
Sep 2, 2017	Dissipation

	unique_vessels	track_segments	avg_segments_per_vessel
date
2017-08-20 00:00:00+00:00	288	1172	4.069444
2017-08-21 00:00:00+00:00	275	1191	4.330909
2017-08-22 00:00:00+00:00	258	1123	4.352713
2017-08-23 00:00:00+00:00	277	1206	4.353791
2017-08-24 00:00:00+00:00	299	1429	4.779264
2017-08-25 00:00:00+00:00	115	339	2.947826
2017-08-26 00:00:00+00:00	90	361	4.011111
2017-08-27 00:00:00+00:00	71	303	4.267606
2017-08-28 00:00:00+00:00	115	422	3.669565
2017-08-29 00:00:00+00:00	90	298	3.311111
2017-08-30 00:00:00+00:00	327	1015	3.103976
2017-08-31 00:00:00+00:00	213	582	2.732394
2017-09-01 00:00:00+00:00	566	2053	3.627208
2017-09-02 00:00:00+00:00	225	824	3.662222
2017-09-03 00:00:00+00:00	212	871	4.108491
2017-09-04 00:00:00+00:00	245	998	4.073469
2017-09-05 00:00:00+00:00	257	1002	3.898833

	vessel_count	segment_count	segments_per_vessel	vessel_ma3	vessel_ma7	vessel_change	vessel_pct_change	vessel_deviation	day_of_week
date
2017-08-01 00:00:00+00:00	615	1822	2.962602	615.000000	615.000000	0.0	0.000000	0.000000	1
2017-08-02 00:00:00+00:00	323	1243	3.848297	469.000000	469.000000	-292.0	-0.474797	-0.311301	2
2017-08-03 00:00:00+00:00	398	1828	4.592965	445.333333	445.333333	75.0	0.232198	-0.106287	3
2017-08-04 00:00:00+00:00	503	2526	5.021869	408.000000	459.750000	105.0	0.263819	0.094073	4
2017-08-05 00:00:00+00:00	353	1943	5.504249	418.000000	438.400000	-150.0	-0.298211	-0.194799	5
2017-08-06 00:00:00+00:00	302	1304	4.317881	386.000000	415.666667	-51.0	-0.144476	-0.273456	6
2017-08-07 00:00:00+00:00	295	1389	4.708475	316.666667	398.428571	-7.0	-0.023179	-0.259591	0
2017-08-08 00:00:00+00:00	276	1188	4.304348	291.000000	350.000000	-19.0	-0.064407	-0.211429	1
2017-08-09 00:00:00+00:00	315	1483	4.707937	295.333333	348.857143	39.0	0.141304	-0.097052	2
2017-08-10 00:00:00+00:00	477	2726	5.714885	356.000000	360.142857	162.0	0.514286	0.324474	3

	vessel_count	vessel_change	anomaly_score
date
2017-09-01 00:00:00+00:00	566	353.0	-0.152955
2017-08-30 00:00:00+00:00	327	237.0	-0.146924
2017-09-07 00:00:00+00:00	647	394.0	-0.132424
2017-08-01 00:00:00+00:00	615	0.0	-0.060690
2017-08-28 00:00:00+00:00	115	44.0	-0.056371
2017-08-25 00:00:00+00:00	115	-184.0	-0.051441
2017-08-10 00:00:00+00:00	477	162.0	-0.037558
2017-08-27 00:00:00+00:00	71	-19.0	-0.036667
2017-08-02 00:00:00+00:00	323	-292.0	-0.014254

GeoAI for Maritime Traffic Disruption During Hurricane Harvey (2017)¶

Learning Objectives¶

Overview¶

Why This Matters (Blue Economy Connection)¶

Hurricane Harvey Timeline¶

Section 1: Environment Setup¶

Section 2: Load AIS Vessel Track Data¶

We first examine the metadata to see what layers it has:¶

As the Hurricane event happened during August and September, we will extract the data from this two months¶

Combine August and September Data¶

Filter AIS data to only blue economy vessels (Cargo, Tanker, Passenger, Fishing) and compute their proportion of the dataset.¶

Section 3: Load Hurricane Harvey Track Data¶

Load Hurricane Harvey track shapefiles (points, path, and wind radii) into GeoDataFrames¶

Convert hurricane DTG timestamps into UTC datetimes and standardize storm attribute column names.¶

Mark the hurricane landfall time and assign a consistent CRS. This adds an explicit landfall event for temporal analysis and ensures the hurricane data uses the same coordinate system as AIS data for spatial comparisons.¶

Section 4: Spatial Analysis - Detecting Traffic Disruptions¶

4.1 Filter AIS Data to Texas Coast Study Area¶

Filter AIS data to the Texas coast study region using a spatial bounding box. This reduces noise outside the study area and ensures later hurricane impact analysis focuses only on vessels potentially affected by Harvey.¶

4.2 Calculate Daily Vessel Traffic Metrics¶

Aggregate AIS data into daily vessel traffic metrics (unique vessels and track density). This converts raw AIS tracks into time-series features so we can observe traffic patterns and detect disruptions around Hurricane Harvey landfall.¶

4.3 Port Region Analysis¶

Focus the Texas-coast AIS data on two key port areas (Houston and Corpus Christi), compute daily vessel traffic for each, and compare pre-/during-/post-Harvey averages. This helps detect localized port disruptions instead of only broad regional changes.¶

Exercise: Analyze Galveston Port Traffic¶

Section 5: Visualizations¶

5.1 Vessel Traffic Time Series¶

Visualize daily unique vessel traffic during the Harvey study period and overlay key hurricane events (landfalls and phases). This makes traffic disruptions easy to see and lets us quantify the average traffic drop during the storm compared to pre-storm conditions.¶

5.2 Port Comparison: Houston vs Corpus Christi¶

5.3 Storm Track and Vessel Movements¶

Section 6: GeoAI Anomaly Detection¶

6.1 Feature Engineering¶

Build a daily feature matrix for anomaly detection using traffic counts, rolling averages, and day-to-day changes. These features turn raw vessel metrics into signals that help a GeoAI model spot unusual drops or spikes (like disruptions caused by Hurricane Harvey).¶

6.2 Train Isolation Forest Anomaly Detector¶

Isolation Forest is an unsupervised anomaly detection algorithm that identifies unusual data points by measuring how easily they can be isolated from the rest of the data.¶

6.3 Visualize Anomaly Detection Results¶