# Core libraries
import pandas as pd	# pandas is a library for data manipulation and analysis using dataframes
import numpy as np	# numpy provides support for large, multi-dimensional arrays and matrices, along with mathematical functions
import geopandas as gpd	# geopandas extends pandas to allow spatial operations on geometric data types
from shapely.geometry import Point	# shapely is a library for manipulation and analysis of planar geometric objects, with Point representing a single coordinate
from pathlib import Path	# pathlib offers an object-oriented approach for handling filesystem paths
import warnings	# warnings is a built-in module to manage warning messages
warnings.filterwarnings('ignore')	# disables showing warnings in the output

# Machine Learning
from sklearn.model_selection import train_test_split	# train_test_split in scikit-learn is used to split data arrays into random train and test subsets
from sklearn.metrics import roc_auc_score, classification_report, confusion_matrix, roc_curve	# these functions from scikit-learn provide tools for evaluating classification models
import xgboost as xgb	# xgboost is an optimized distributed gradient boosting library designed for high performance and efficiency
import shap	# shap is a library for interpreting machine learning models using SHAP values (SHapley Additive exPlanations)

# Visualization
import matplotlib.pyplot as plt	# matplotlib.pyplot is a collection of functions for creating static, interactive, and animated visualizations in Python
import seaborn as sns	# seaborn is a statistical data visualization library built on top of matplotlib

import plotly.express as px	# plotly.express is a high-level interface for creating interactive visualizations easily

# Set random seed for reproducibility
np.random.seed(42)

# STUDY AREA CONFIGURATION
# Load Essential Fish Habitat (EFH) shapefile
EFH_SHAPEFILE = Path('shapefile-reef-fish-efh-gomex-sero/ReefFish_EFH_GOM.shp')	# defines the path to the Essential Fish Habitat (EFH) shapefile
EFH_GDF = gpd.read_file(EFH_SHAPEFILE)	# reads the EFH shapefile into a GeoDataFrame using geopandas
EFH_GDF = EFH_GDF.to_crs('EPSG:4326')	# converts the GeoDataFrame coordinate reference system to WGS84 latitude/longitude (EPSG:4326)

# Study bounds from EFH polygon
efh_bounds = EFH_GDF.total_bounds	# extracts the bounding coordinates (minx, miny, maxx, maxy) of all geometries in the GeoDataFrame
BOUNDS = {	# creates a dictionary containing the bounding box values for longitude and latitude
    'lon_min': efh_bounds[0], 'lon_max': efh_bounds[2],
    'lat_min': efh_bounds[1], 'lat_max': efh_bounds[3]
}

# Time period and target species
TIME_PERIOD = {'start': '2000-01-01', 'end': '2023-12-31'}	# specifies the time period for the study
TARGET_SPECIES = 'Lutjanus campechanus'	# scientific name of the target species (Red Snapper)
COMMON_NAME = 'Red Snapper'	# common name of the target species

# Create output directories
DATA_DIR = Path('data')	# defines the path where data files will be saved
FIG_DIR = Path('figures')	# defines the path where figures will be saved
for d in [DATA_DIR, FIG_DIR]:	# iterates over the list of output directory paths
    d.mkdir(parents=True, exist_ok=True)	# creates the required directories if they do not exist

print("Setup complete!")
print(f"Study area: Gulf of Mexico Reef Fish EFH")
print(f"Target: {COMMON_NAME} ({TARGET_SPECIES})")
print(f"Period: {TIME_PERIOD['start']} to {TIME_PERIOD['end']}")

Setup complete!
Study area: Gulf of Mexico Reef Fish EFH
Target: Red Snapper (Lutjanus campechanus)
Period: 2000-01-01 to 2023-12-31

"STATIONID","CRUISEID","DECSLAT","DECSLON","DEPTH_SSTA","START_DATE"
"1","581","26.495","-96.505","147.2","11/10/2003 02:49:00"
"2","581","26.049","-96.471","119.4","11/10/2003 06:19:00"
"3","581","25.989","-96.985","51.8","11/10/2003 10:29:00"

SEAMAP_DIR = Path('SEAMAPDATAV3CSV')	# Path is a pathlib object that defines the directory containing SEAMAP data files
RED_SNAPPER_BIOCODE = 170151107	# Sets the unique biocode integer for Red Snapper species

print("LOADING SEAMAP DATA")

# Step 1: Load Station Records
print("\n1. Loading station records...")
starec = pd.read_csv(SEAMAP_DIR / 'starec.csv')	# pandas is a data analysis library; reads the station records CSV into a DataFrame
starec = starec[['STATIONID', 'CRUISEID', 'DECSLAT', 'DECSLON', 'DEPTH_SSTA', 'START_DATE']].copy()	# Selects and copies columns relevant to station metadata from the DataFrame
starec.columns = ['station_id', 'cruise_id', 'latitude', 'longitude', 'depth_m', 'start_date']	# Renames the selected columns to use standardized lowercase names

starec['latitude'] = pd.to_numeric(starec['latitude'], errors='coerce')	# Converts latitude values to numeric, coercing errors to NaN for invalid entries
starec['longitude'] = pd.to_numeric(starec['longitude'], errors='coerce')	# Converts longitude values to numeric, coercing errors to NaN for invalid entries
starec['depth_m'] = pd.to_numeric(starec['depth_m'], errors='coerce')	# Converts depth values to numeric (meters), coercing errors to NaN
starec['date'] = pd.to_datetime(starec['start_date'], format='mixed', errors='coerce')	# Parses the start_date column to pandas datetime, allowing mixed date formats and coercing errors to NaT
starec = starec.dropna(subset=['latitude', 'longitude', 'date'])	# Removes rows with missing latitude, longitude, or date fields
print(f"   Loaded {len(starec):,} station records")

LOADING SEAMAP DATA

1. Loading station records...
   Loaded 61,384 station records

# Step 2: Filter by time period
print("\n2. Filtering by time period...")
start_date = pd.to_datetime(TIME_PERIOD['start'])	# Converts the start date string in TIME_PERIOD to pandas datetime
end_date = pd.to_datetime(TIME_PERIOD['end'])	# Converts the end date string in TIME_PERIOD to pandas datetime
time_mask = (starec['date'] >= start_date) & (starec['date'] <= end_date)	# Creates a boolean mask to filter station records within the time period
starec = starec[time_mask].copy()	# Applies the time mask and makes a copy of the filtered DataFrame
print(f"   Stations in time period: {len(starec):,}")

2. Filtering by time period...
   Stations in time period: 31,132

# Step 3: Filter by EFH polygon
print("\n3. Filtering by EFH polygon...")
geometry = [Point(lon, lat) for lon, lat in zip(starec.longitude, starec.latitude)]	# Point is a shapely geometry object; constructs a list of Point objects from longitude and latitude
starec_gdf = gpd.GeoDataFrame(starec, geometry=geometry, crs='EPSG:4326')	# geopandas reads tabular data as a GeoDataFrame with coordinate reference system set to WGS84
stations = gpd.sjoin(starec_gdf, EFH_GDF, predicate='within', how='inner')	# Performs a spatial join to retain only stations within the Essential Fish Habitat polygon
stations = stations.drop(columns=['index_right', 'Area_SqKm', 'Perim_M'], errors='ignore')	# Drops nonessential columns resulting from the spatial join if present
print(f"   Stations within EFH: {len(stations):,}")

3. Filtering by EFH polygon...
   Stations within EFH: 27,881

# Step 4: Load Catch Records
print("\n4. Loading catch records...")
bgsrec = pd.read_csv(SEAMAP_DIR / 'bgsrec.csv', low_memory=False)	# Reads the catch records CSV file into a DataFrame, disabling low_memory mode for consistency
bgsrec = bgsrec[['STATIONID', 'BIO_BGS', 'CNT', 'CNTEXP']].copy()	# Selects relevant columns pertaining to catch counts and biological code
bgsrec.columns = ['station_id', 'biocode', 'count', 'count_expanded']	# Renames the columns for clarity and ease of use
bgsrec['abundance'] = bgsrec['count_expanded'].fillna(bgsrec['count'])	# Fills missing expanded counts with raw counts to compute total abundance
bgsrec['abundance'] = pd.to_numeric(bgsrec['abundance'], errors='coerce').fillna(0).astype(int)	# Ensures abundance is stored as numeric, filling invalid entries with zero and converting to integer
print(f"   Loaded {len(bgsrec):,} catch records")

4. Loading catch records...
   Loaded 891,280 catch records

# Step 5: Filter for Red Snapper
print(f"\n5. Filtering for Red Snapper...")
target_catch = bgsrec[bgsrec['biocode'] == RED_SNAPPER_BIOCODE].copy()	# Filters catch records for rows with the Red Snapper biocode
target_by_station = target_catch.groupby('station_id')['abundance'].sum().reset_index()	# Aggregates abundance of Red Snapper by station using groupby and sum
print(f"   Found {len(target_catch):,} Red Snapper records")

5. Filtering for Red Snapper...
   Found 12,583 Red Snapper records

# Step 6: Create presence/absence dataset
print("\n6. Create presence/absence dataset...")
fish_gdf = stations.merge(target_by_station, on='station_id', how='left')	# Merges the spatially filtered stations with Red Snapper catch data on station_id using a left join
fish_gdf['abundance'] = fish_gdf['abundance'].fillna(0).astype(int)	# Assigns a value of 0 to stations with missing abundance (no Red Snapper caught) and ensures integer type
fish_gdf['presence'] = (fish_gdf['abundance'] > 0).astype(int)	# Creates a binary column indicating presence (1) or absence (0) of Red Snapper

6. Create presence/absence dataset...

# Step 7: Load Environmental Data
print("\n7. Loading environmental data...")
envrec = pd.read_csv(SEAMAP_DIR / 'envrec.csv', low_memory=False)	# Reads the environmental measurements CSV file into a DataFrame
env_cols = ['STATIONID', 'TEMPSURF', 'TEMPMID', 'TEMPMAX', 	# Specifies the list of columns to extract for environmental variables
            'SALSURF', 'SALMID', 'SALMAX',
            'OXYSURF', 'OXYMID', 'OXYMAX',
            'CHLORSURF', 'CHLORMID', 'CHLORMAX']
envrec = envrec[env_cols].copy()	# Selects only the relevant environmental columns
envrec.columns = ['station_id', 'temp_surf', 'temp_mid', 'temp_bot',	# Renames columns to concise, standardized names for downstream analysis
                  'sal_surf', 'sal_mid', 'sal_bot',
                  'oxy_surf', 'oxy_mid', 'oxy_bot',
                  'chl_surf', 'chl_mid', 'chl_bot']
for col in envrec.columns[1:]:	# Iterates over all environmental columns except station_id
    envrec[col] = pd.to_numeric(envrec[col], errors='coerce')	# Converts each environmental variable column to numeric, coercing errors to NaN

fish_gdf = fish_gdf.merge(envrec, on='station_id', how='left')	# Merges the environmental data into the fish GeoDataFrame based on station_id

# Summary
n_total = len(fish_gdf)	# Computes the total number of survey stations in the merged dataset
n_presence = fish_gdf['presence'].sum()	# Computes the total number of stations with Red Snapper present
print("DATA LOADED SUCCESSFULLY")
print(f"Total stations: {n_total:,}")
print(f"Presence: {n_presence:,} ({100*n_presence/n_total:.1f}%)")
print(f"Absence: {n_total - n_presence:,} ({100*(n_total-n_presence)/n_total:.1f}%)")
fish_gdf.head(5)	# Displays the first five rows of the final fish_gdf DataFrame for inspection

7. Loading environmental data...
DATA LOADED SUCCESSFULLY
Total stations: 27,922
Presence: 6,488 (23.2%)
Absence: 21,434 (76.8%)

# Visualize survey data
fig, axes = plt.subplots(1, 2, figsize=(14, 5))	# Creates a matplotlib figure with two subplot axes, with specified dimensions for plotting
presence = fish_gdf[fish_gdf['presence'] == 1]	# Filters the survey data for sites where Red Snapper were present
absence = fish_gdf[fish_gdf['presence'] == 0]	# Filters the survey data for sites where Red Snapper were absent

# Plot 1: Spatial distribution
ax1 = axes[0]	# Selects the first subplot axis for spatial plotting
EFH_GDF.boundary.plot(ax=ax1, color='navy', linewidth=1.5, label='EFH Boundary')	# Plots the Essential Fish Habitat boundary as a line on the map
ax1.scatter(absence.longitude, absence.latitude, c='lightgray', s=15, alpha=0.4, label='Absence')	# Plots absence points as light gray circles on the map
ax1.scatter(presence.longitude, presence.latitude, c='crimson', s=25, alpha=0.7, label='Presence')	# Plots presence points as crimson circles on the map
ax1.set_xlabel('Longitude'); ax1.set_ylabel('Latitude')	# Sets the x and y axis labels for spatial plot
ax1.set_title(f'{COMMON_NAME} Survey Locations')	# Sets the subplot title showing the fish common name
ax1.legend(loc='lower right')	# Adds a legend for the spatial plot in the lower right corner

# Plot 2: Depth distribution
ax2 = axes[1]	# Selects the second subplot axis for depth distribution
ax2.hist(absence['depth_m'], bins=30, alpha=0.5, label='Absence', color='gray', density=True, range=(0, 200))	# Plots a normalized histogram of depths for absence locations
ax2.hist(presence['depth_m'], bins=30, alpha=0.7, label='Presence', color='crimson', density=True, range=(0, 200))	# Plots a normalized histogram of depths for presence locations
ax2.set_xlabel('Depth (m)'); ax2.set_ylabel('Density')	# Sets x and y labels for the histogram plot
ax2.set_title('Depth Distribution by Presence/Absence')	# Assigns a title to the depth distribution plot
ax2.legend()	# Adds a legend to the depth distribution subplot

plt.tight_layout()	# Adjusts the layout to prevent overlap of subplots and labels
plt.savefig(FIG_DIR / 'fish_survey_exploration.png', dpi=150)	# Saves the generated figure as a PNG file in the figure output directory
plt.show()	# Displays the plots

print(f"Total surveys: {len(fish_gdf):,} | Presence: {len(presence):,} ({100*len(presence)/len(fish_gdf):.1f}%)")

Total surveys: 27,922 | Presence: 6,488 (23.2%)

# Environmental conditions by presence/absence
fig, axes = plt.subplots(1, 3, figsize=(14, 4))	# matplotlib is a plotting library for Python, this line creates a figure with three subplots horizontally and sets the figure size.
presence_data = fish_gdf[fish_gdf['presence'] == 1]	# Selects rows from the GeoDataFrame where Red Snapper are present, creating a new DataFrame for these locations.
absence_data = fish_gdf[fish_gdf['presence'] == 0]	# Selects rows from the GeoDataFrame where Red Snapper are absent, creating a corresponding DataFrame.

# Bottom Temperature
ax = axes[0]	# Assigns the first subplot axis (for bottom temperature histograms) to the variable ax.
ax.hist(absence_data['temp_bot'].dropna(), bins=25, alpha=0.5, label='Absence', color='gray', density=True)	# Plots a normalized histogram of bottom temperature for absence locations, ignoring missing values.
ax.hist(presence_data['temp_bot'].dropna(), bins=25, alpha=0.7, label='Presence', color='crimson', density=True)	# Plots a normalized histogram of bottom temperature for presence locations, ignoring missing values.
ax.set_xlabel('Bottom Temperature (C)'); ax.set_ylabel('Density')	# Sets the x-axis label as bottom temperature in Celsius and the y-axis label as density for the first subplot.
ax.set_title('Bottom Temperature'); ax.legend()	# Sets the title for the first subplot and adds a legend explaining labels.

# Bottom Oxygen  
ax = axes[1]	# Assigns the second subplot axis (for bottom oxygen histograms) to the variable ax.
ax.hist(absence_data['oxy_bot'].dropna(), bins=25, alpha=0.5, label='Absence', color='gray', density=True)	# Plots a normalized histogram of bottom oxygen for absence locations, skipping missing values.
ax.hist(presence_data['oxy_bot'].dropna(), bins=25, alpha=0.7, label='Presence', color='crimson', density=True)	# Plots a normalized histogram of bottom oxygen for presence locations, skipping missing values.
ax.axvline(x=2, color='red', linestyle='--', label='Hypoxia')	# Draws a vertical dashed red line at oxygen = 2 ml/L to indicate the hypoxia threshold.
ax.set_xlabel('Bottom Oxygen (ml/L)'); ax.set_ylabel('Density')	# Sets the x-axis label as bottom oxygen in milliliters per liter and the y-axis label as density for the second subplot.
ax.set_title('Bottom Oxygen'); ax.legend()	# Sets the title for the second subplot and adds a legend.

# Depth
ax = axes[2]	# Assigns the third subplot axis (for depth histograms) to the variable ax.
ax.hist(absence_data['depth_m'].dropna(), bins=25, alpha=0.5, label='Absence', color='gray', density=True, range=(0,200))	# Plots a normalized histogram of depth for absence locations, restricting the range to 0-200 meters.
ax.hist(presence_data['depth_m'].dropna(), bins=25, alpha=0.7, label='Presence', color='crimson', density=True, range=(0,200))	# Plots a normalized histogram of depth for presence locations, restricting the range to 0-200 meters.
ax.set_xlabel('Depth (m)'); ax.set_ylabel('Density')	# Sets the x-axis label as depth in meters and the y-axis label as density for the third subplot.
ax.set_title('Depth'); ax.legend()	# Sets the title for the third subplot and adds a legend.

plt.suptitle(f'{COMMON_NAME}: Environmental Conditions', fontweight='bold')	# Sets an overall title for the figure showing the common name of the fish, in bold font weight.
plt.tight_layout()	# Automatically adjusts subplot parameters to give specified padding and prevent overlaps.
plt.savefig(FIG_DIR / 'environmental_distributions.png', dpi=150)	# Saves the current figure as a PNG image with 150 dpi resolution to the designated figure directory.
plt.show()	# Displays the generated matplotlib figure to the screen.

print(f"Key finding: Red Snapper occur at warmer bottom temps ({presence_data['temp_bot'].mean():.1f}C vs {absence_data['temp_bot'].mean():.1f}C) and shallower depths ({presence_data['depth_m'].mean():.0f}m vs {absence_data['depth_m'].mean():.0f}m)")

Key finding: Red Snapper occur at warmer bottom temps (24.8C vs 23.1C) and shallower depths (39m vs 48m)

# Add temporal features and prepare training data
fish_gdf['month'] = pd.to_datetime(fish_gdf['date']).dt.month	# Converts the 'date' column to datetime and extracts the month as a new column in the GeoDataFrame.
fish_gdf['year'] = pd.to_datetime(fish_gdf['date']).dt.year	# Converts the 'date' column to datetime and extracts the year as a new column in the GeoDataFrame.

fish_gdf['month_sin'] = np.sin(2 * np.pi * fish_gdf['month'] / 12)	# Calculates the sine of the month (cyclical encoding) to represent seasonal patterns for modeling.
fish_gdf['month_cos'] = np.cos(2 * np.pi * fish_gdf['month'] / 12)	# Calculates the cosine of the month (cyclical encoding) to represent seasonal patterns for modeling.

# Define features for the model
FEATURE_COLS = ['temp_bot', 'temp_surf', 'sal_bot', 'oxy_bot', 'depth_m', 'month_sin', 'month_cos']	# Lists the selected feature column names that will be used as input variables for the model.
TARGET_COL = 'presence'	# Specifies the target variable column name for model prediction.

# Remove rows with missing environmental data
df_clean = fish_gdf.dropna(subset=FEATURE_COLS + [TARGET_COL]).copy()	# Creates a new DataFrame by removing rows that have missing values in any of the selected feature or target columns.

print(f"Training data prepared:")
print(f"Features: {FEATURE_COLS}")
print(f"Records: {len(df_clean):,} (dropped {len(fish_gdf)-len(df_clean):,} with missing data)")
print(f"Presence rate: {df_clean[TARGET_COL].mean()*100:.1f}%")

Training data prepared:
Features: ['temp_bot', 'temp_surf', 'sal_bot', 'oxy_bot', 'depth_m', 'month_sin', 'month_cos']
Records: 24,092 (dropped 3,830 with missing data)
Presence rate: 23.6%

train_years = [year for year in range(2000, 2015)]	# Creates a list of years from 2000 to 2014 for the training set in a temporal split.
test_years = [year for year in range(2016, 2023)]	# Creates a list of years from 2016 to 2022 for the testing set in a temporal split.

train_mask = df_clean['year'].isin(train_years)	# Generates a Boolean mask identifying rows in the DataFrame where the 'year' value is in the training years.
test_mask = df_clean['year'].isin(test_years)	# Generates a Boolean mask for rows where the 'year' value belongs to the test years.

X_train = df_clean.loc[train_mask, FEATURE_COLS]	# Selects the feature columns for training samples from the DataFrame using the training mask.
y_train = df_clean.loc[train_mask, TARGET_COL]	# Selects the target column for the training set using the training mask.
X_test = df_clean.loc[test_mask, FEATURE_COLS]	# Selects the feature columns for test samples from the DataFrame using the test mask.
y_test = df_clean.loc[test_mask, TARGET_COL]	# Selects the target column for the test set using the test mask.

print(f"Train/Test Split (Temporal):")
print(f"Training: {len(X_train)} samples ({train_years[0]}-{train_years[-1]})")
print(f"Testing:  {len(X_test)} samples ({test_years[0]}-{test_years[-1]})")
print(f"\nTrain presence rate: {y_train.mean()*100:.1f}%")
print(f"Test presence rate:  {y_test.mean()*100:.1f}%")

Train/Test Split (Temporal):
Training: 18360 samples (2000-2014)
Testing:  4124 samples (2016-2022)

Train presence rate: 20.9%
Test presence rate:  33.0%

# Train XGBoost classifier
xgb_model = xgb.XGBClassifier(	# XGBoost is an optimized gradient boosting library used for supervised learning tasks.
    n_estimators=4000, max_depth=20, learning_rate=0.1,	# Sets the number of boosting rounds, the maximum depth of trees, and the learning rate.
    subsample=0.8, colsample_bytree=0.8, random_state=42,	# Specifies subsampling ratio of the training instances, subsampling ratio of columns, and ensures reproducibility by setting a random seed.
    eval_metric='auc',
    early_stopping_rounds=30	# Uses the AUC metric for evaluation and enables early stopping if the score does not improve for 20 rounds.
)

xgb_model.fit(X_train, y_train, eval_set=[(X_test, y_test)], verbose=False)	# Trains the XGBoost classifier on training data and monitors performance on test data for early stopping.
print(f"XGBoost trained with {xgb_model.best_iteration} trees")

XGBoost trained with 53 trees

# Defines a function to evaluate a binary classifier, visualize results, and optionally save the plot.
def evaluate_model(model, X_test, y_test, model_name='Model', save_fig=True, fig_filename='model_evaluation.png'):
    """
    Evaluate a binary classifier and visualize results.
    
    Parameters:
    -----------
    model : classifier with predict_proba method
        Trained model to evaluate
    X_test : DataFrame or array
        Test features
    y_test : Series or array
        Test labels (0/1)
    model_name : str
        Name of the model for plot labels
    save_fig : bool
        Whether to save the figure
    fig_filename : str
        Filename for saved figure
        
    Returns:
    --------
    dict : Dictionary containing AUC score and predictions
    """
    y_pred_prob = model.predict_proba(X_test)[:, 1]	# Generates probability predictions for the positive class using the model's predict_proba method.
    y_pred = (y_pred_prob > 0.5).astype(int)	# Binarizes the predicted probabilities using a threshold of 0.5 to obtain 0/1 predictions.
    auc_score = roc_auc_score(y_test, y_pred_prob)	# Computes the Area Under the ROC Curve (AUC) using true labels and predicted probabilities.

    print(f"\n{model_name} Performance: AUC-ROC = {auc_score:.3f}")

    fig, axes = plt.subplots(1, 2, figsize=(12, 4))	# Creates a matplotlib figure with two subplots for ROC curve and confusion matrix.

    ax = axes[0]	# Selects the first subplot axis for ROC curve plotting.
    fpr, tpr, _ = roc_curve(y_test, y_pred_prob)	# Computes the false positive rate and true positive rate for ROC curve plotting.
    ax.plot(fpr, tpr, 'b-', linewidth=2, label=f'{model_name} (AUC={auc_score:.3f})')	# Plots the ROC curve with AUC annotation.
    ax.plot([0, 1], [0, 1], 'k--', linewidth=1, label='Random')	# Plots a diagonal line as the random performance baseline.
    ax.set_xlabel('False Positive Rate')	# Sets the x-axis label for the ROC curve.
    ax.set_ylabel('True Positive Rate')	# Sets the y-axis label for the ROC curve.
    ax.set_title(f'{model_name} ROC Curve')	# Sets the plot title for the ROC curve.
    ax.legend()	# Displays the legend on the ROC plot.
    ax.grid(True, alpha=0.3)	# Enables a semi-transparent grid on the ROC plot.

    ax = axes[1]	# Selects the second subplot axis for the confusion matrix.
    cm = confusion_matrix(y_test, y_pred)	# Computes the confusion matrix from the true and predicted labels.
    sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=ax,	# Uses seaborn to plot the confusion matrix as a heatmap.
                xticklabels=['Pred Absent', 'Pred Present'],	# Sets the x-axis tick labels for predicted classes.
                yticklabels=['Actual Absent', 'Actual Present'])	# Sets the y-axis tick labels for actual classes.
    ax.set_title(f'{model_name} Confusion Matrix')	# Sets the plot title for the confusion matrix.

    plt.tight_layout()	# Automatically adjusts subplot parameters for better layout.

    if save_fig:	# Checks whether to save the figure to disk.
        plt.savefig(FIG_DIR / fig_filename, dpi=150)	# Saves the figure as a file in the specified directory with 150 dpi.

    plt.show()	# Displays the figure window with the evaluation plots.

    print(f"\n{model_name} Classification Report:")
    print(classification_report(y_test, y_pred, target_names=['Absent', 'Present']))

    return {	# Returns a dictionary of evaluation results including AUC, binary predictions, and predicted probabilities.
        'auc': auc_score,	# Stores the computed AUC score.
        'predictions': y_pred,	# Stores the binary predictions.
        'probabilities': y_pred_prob	# Stores the predicted probabilities for the positive class.
    }

print("Model evaluation function created!")

Model evaluation function created!

# Use the evaluation function for XGBoost
xgb_results = evaluate_model(
    model=xgb_model,
    X_test=X_test,
    y_test=y_test,
    model_name='XGBoost',
    save_fig=True,
    fig_filename='xgboost_evaluation.png'
)

XGBoost Performance: AUC-ROC = 0.705

XGBoost Classification Report:
              precision    recall  f1-score   support

      Absent       0.73      0.85      0.79      2764
     Present       0.54      0.35      0.43      1360

    accuracy                           0.69      4124
   macro avg       0.64      0.60      0.61      4124
weighted avg       0.67      0.69      0.67      4124

# ANSWER

# Step 1: Import RandomForestClassifier
from sklearn.ensemble import RandomForestClassifier

# Step 2: Create and train the Random Forest model
rf_model = RandomForestClassifier(
    n_estimators=200,        # Number of trees in the forest
    max_depth=15,            # Maximum depth of each tree
    min_samples_split=10,    # Minimum samples required to split a node
    random_state=42,         # For reproducibility
    n_jobs=-1                # Use all CPU cores for faster training
)

# Train the model
print("Training Random Forest model...")
rf_model.fit(X_train, y_train)
print(f"Random Forest trained with {rf_model.n_estimators} trees")

# Step 3: Evaluate the model using the evaluate_model() function
rf_results = evaluate_model(
    model=rf_model,
    X_test=X_test,
    y_test=y_test,
    model_name='Random Forest',
    save_fig=True,
    fig_filename='random_forest_evaluation.png'
)

Training Random Forest model...
Random Forest trained with 200 trees

Random Forest Performance: AUC-ROC = 0.716

Random Forest Classification Report:
              precision    recall  f1-score   support

      Absent       0.73      0.86      0.79      2764
     Present       0.55      0.34      0.42      1360

    accuracy                           0.69      4124
   macro avg       0.64      0.60      0.61      4124
weighted avg       0.67      0.69      0.67      4124

# Calculate and visualize SHAP values
explainer = shap.TreeExplainer(xgb_model)	# shap is a library for model explainability; this line creates a SHAP TreeExplainer object for the trained XGBoost model.
shap_values = explainer.shap_values(X_test)	# Computes the SHAP values for the test set features, quantifying each feature's contribution to the predictions.

# SHAP Summary Plot - shows feature importance and direction of effects
plt.figure(figsize=(10, 6))	# matplotlib.pyplot is used for plotting; this initializes a new figure with a specified size.
shap.summary_plot(shap_values, X_test, show=False)	# Generates a summary plot visualizing feature importance and the effect direction, using SHAP values and the test set features.
plt.title(f'{COMMON_NAME} Habitat Model: Feature Importance (SHAP)')	# Sets the title of the plot using the species common name.
plt.tight_layout()	# Adjusts plot spacing to prevent overlap.
plt.savefig(FIG_DIR / 'shap_summary.png', dpi=150)	# Saves the current figure to the designated FIG_DIR as a PNG file with 150 dpi resolution.
plt.show()	# Displays the figure.

print("Interpretation: Features on top are most important. Red = high values, Blue = low values.")
print("Positive SHAP = increases habitat suitability, Negative SHAP = decreases suitability.")

Interpretation: Features on top are most important. Red = high values, Blue = low values.
Positive SHAP = increases habitat suitability, Negative SHAP = decreases suitability.

df_predict = df_clean.copy()	# Creates a copy of the cleaned dataframe to retain the original data and allow modifications for prediction.
df_predict['habitat_prob'] = xgb_model.predict_proba(df_predict[FEATURE_COLS])[:, 1]	# Uses the trained XGBoost model to predict the probability of suitable habitat for each row, storing the result in a new column.
df_predict['habitat_class'] = (df_predict['habitat_prob'] > 0.5).astype(int)	# Converts habitat probability to a binary class where values above 0.5 are suitable (1) and others are unsuitable (0).

fig, ax = plt.subplots(1, 1, figsize=(12, 8))

# Add hexbin heatmap
hexbin = ax.hexbin(	# Creates a hexagonal binning plot on the axes to spatially aggregate predicted habitat probabilities.
    df_predict['longitude'],	# Supplies longitude data for the hexbin x-coordinates.
    df_predict['latitude'],	# Supplies latitude data for the hexbin y-coordinates.
    C=df_predict['habitat_prob'],	# Sets the value to be aggregated within each hex grid cell as predicted habitat probability.
    gridsize=35,	# Specifies the number of hexagons in the x-direction (resolution of spatial aggregation).
    cmap='RdYlGn',	# Applies the RdYlGn colormap for the hexbin color mapping, visualizing suitability from low to high.
    reduce_C_function=np.mean,	# Uses NumPy's mean function to compute the average predicted suitability in each hexagon.
    mincnt=1,	# Renders hexagons only where at least one point falls within a grid cell.
    alpha=0.8,	# Sets the transparency of the hexagons for visual effect.
    vmin=0,	# Sets the minimum colorbar value to 0 for normalization.
    vmax=1	# Sets the maximum colorbar value to 1 for normalization.
)

EFH_GDF.boundary.plot(ax=ax, color='navy', linewidth=1.5, alpha=0.9, zorder=2)	# Plots the Essential Fish Habitat (EFH) boundary on the axes, overlaying it above the hexbin layer for spatial context.

plt.colorbar(hexbin, ax=ax, label='Mean Habitat Suitability (0-1)')	# Adds a colorbar to the plot to indicate the scale of mean predicted habitat suitability.

ax.set_xlabel('Longitude')	# Sets the x-axis label to 'Longitude'.
ax.set_ylabel('Latitude')	# Sets the y-axis label to 'Latitude'.
ax.set_title(f'Red Snapper Habitat Suitability Heatmap') 	# Sets the plot title

plt.tight_layout()	# Adjusts subplot parameters to prevent overlapping plot elements.
plt.savefig(FIG_DIR / 'habitat_suitability_heatmap.png', dpi=150)	# Saves the figure to the designated output directory at 150 dots per inch resolution.
plt.show()	# Displays the plot in the notebook output.

# ENVIRONMENT CONDITION CHANGE SCENARIO SIMULATION

# Use a single recent year (2022) as the baseline for cleaner interpretation
SCENARIO_YEAR = 2022	# Defines the reference year for baseline conditions.
df_scenario = df_clean[df_clean['year'] == SCENARIO_YEAR].copy()	# Filters the cleaned data to only include observations from the scenario year.
print(f"Using {SCENARIO_YEAR} as baseline year: {len(df_scenario):,} survey locations")
print(f"Average bottom temp: {df_scenario['temp_bot'].mean():.1f}°C | Average bottom oxygen: {df_scenario['oxy_bot'].mean():.1f} ml/L")

# Baseline prediction (current conditions for the scenario year)
baseline_probs = xgb_model.predict_proba(df_scenario[FEATURE_COLS])[:, 1]	# Uses the XGBoost model to predict habitat suitability probabilities for the baseline year.

# Combined stress scenario: +2°C warming and -1.5 ml/L oxygen decline
df_stress = df_scenario.copy()	# Creates a copy of the scenario DataFrame
df_stress['temp_bot'] = df_stress['temp_bot'] + 2	# Increases the bottom temperature column by 2°C across all observations.
df_stress['temp_surf'] = df_stress['temp_surf'] + 2	# Increases the surface temperature column by 2°C for all locations.
df_stress['oxy_bot'] = np.clip(df_stress['oxy_bot'] - 1.5, 0.5, 10)	# Reduces bottom oxygen by 1.5 ml/L, ensuring values remain within the range 0.5-10.
stress_probs = xgb_model.predict_proba(df_stress[FEATURE_COLS])[:, 1]	# Predicts suitability probability using the modified stress scenario environmental features.

# Visualize comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 6))	# Matplotlib is used for plotting; creates a figure with two subplots side by side.

for ax, probs, title in [(axes[0], baseline_probs, f'Baseline ({SCENARIO_YEAR})'),
                          (axes[1], stress_probs, 'Environment Stress (+2°C, -1.5 O₂)')]:	# Iterates over subplots and their corresponding data to plot baseline and stress scenario.
    EFH_GDF.boundary.plot(ax=ax, color='navy', linewidth=1, alpha=0.5)	# Plots the Essential Fish Habitat boundary on the map using GeoPandas.
    scatter = ax.scatter(df_scenario['longitude'], df_scenario['latitude'],
                        c=probs, cmap='RdYlGn', s=25, alpha=0.8, vmin=0, vmax=1)	# Creates a scatter plot of survey locations colored by predicted suitability.
    plt.colorbar(scatter, ax=ax, label='Suitability')	# Adds a colorbar to each subplot to display suitability values.
    ax.set_title(title); ax.set_xlabel('Longitude'); ax.set_ylabel('Latitude')	# Sets the title and axis labels for each subplot.

plt.suptitle(f'{COMMON_NAME}: Environmental Stress Impact on Habitat Suitability ({SCENARIO_YEAR} Baseline)', fontweight='bold')	# Sets a super title summarizing the plot's context with targeted species name.
plt.tight_layout()	# Adjusts subplot parameters to improve layout and prevent overlap.
plt.savefig(FIG_DIR / 'Environmental_Stress.png', dpi=150)	# Saves the visualization as a PNG file in the specified directory at 150 dpi.
plt.show()	# Displays the figure.

# Summary statistics
baseline_suit = (baseline_probs > 0.5).sum()	# Counts the number of baseline locations with suitability probability higher than 0.5.
stress_suit = (stress_probs > 0.5).sum()	# Counts the number of locations under stress scenario with suitability higher than 0.5.
lost = baseline_suit - stress_suit	# Calculates the loss in suitable habitat locations due to the environment scenario.
print(f"\Environmental Stress Impact Summary ({SCENARIO_YEAR} baseline):")
print(f"  Baseline suitable: {baseline_suit:,} locations")
print(f"  After stress: {stress_suit:,} locations")
print(f"  Habitat lost: {lost:,} locations ({100*lost/baseline_suit:.1f}% reduction)")

Using 2022 as baseline year: 445 survey locations
Average bottom temp: 25.1°C | Average bottom oxygen: 5.5 ml/L

\Environmental Stress Impact Summary (2022 baseline):
  Baseline suitable: 87 locations
  After stress: 56 locations
  Habitat lost: 31 locations (35.6% reduction)

	station_id	cruise_id	latitude	longitude	depth_m	start_date	date	geometry	abundance	presence	...	temp_bot	sal_surf	sal_mid	sal_bot	oxy_surf	oxy_mid	oxy_bot	chl_surf	chl_mid	chl_bot
0	1	581	26.495	-96.505	147.2	11/10/2003 02:49:00	2003-11-10 02:49:00	POINT (-96.50500 26.49500)	0	0	...	20.98	34.38	36.13	36.43	5.7	5.3	4.4	1.90	15.99	3.57
1	2	581	26.049	-96.471	119.4	11/10/2003 06:19:00	2003-11-10 06:19:00	POINT (-96.47100 26.04900)	0	0	...	24.08	31.96	36.02	36.43	5.7	5.4	4.9	2.42	4.06	4.41
2	4	581	26.086	-97.094	18.3	11/10/2003 13:37:00	2003-11-10 13:37:00	POINT (-97.09400 26.08600)	26	1	...	27.00	30.21	31.28	32.58	5.6	5.3	4.1	3.36	4.34	9.04
3	5	581	26.121	-97.140	13.4	11/10/2003 15:19:00	2003-11-10 15:19:00	POINT (-97.14000 26.12100)	1	1	...	26.84	30.41	31.10	31.66	5.6	4.6	4.2	3.22	9.27	9.35
4	6	581	26.222	-97.155	13.5	11/10/2003 16:59:00	2003-11-10 16:59:00	POINT (-97.15500 26.22200)	0	0	...	26.69	29.94	30.22	30.94	5.6	5.3	4.6	3.42	5.45	11.58

GeoAI for Fish Habitat Suitability Modeling¶

Learning Objectives¶

Case Study: Red Snapper in the Gulf of America¶

Section 1: Setup¶

STUDY AREA CONFIGURATION¶

Section 2: Load & Explore Fish Survey Data¶

SEAMAP survey data overview¶

`starec.csv` – Station (tow) metadata¶

`bgsrec.csv` – Biological catch records¶

`envrec.csv` – Environmental measurements at stations¶

Load the Survey Data¶

Visualize Survey Data¶

Section 3: Visualize Environmental Data¶

Section 4: Prepare Training Data¶

Section 5: Build Habitat Model¶

Train XGBoost¶

Create Model Evaluation Function¶

Evaluate XGBoost Model¶

Section 5.5: Exercise - Build a Random Forest Model¶

Background¶

Your Task¶

Section 6: Model Interpretation with SHAP¶

Section 7: Habitat Suitability Predictions¶

Section 8: Environment Condition Change Analysis¶

GeoAI for Fish Habitat Suitability Modeling¶

Learning Objectives¶

Case Study: Red Snapper in the Gulf of America¶

Section 1: Setup¶

STUDY AREA CONFIGURATION¶

Section 2: Load & Explore Fish Survey Data¶

SEAMAP survey data overview¶

starec.csv – Station (tow) metadata¶

bgsrec.csv – Biological catch records¶

envrec.csv – Environmental measurements at stations¶

Load the Survey Data¶

Visualize Survey Data¶

Section 3: Visualize Environmental Data¶

Section 4: Prepare Training Data¶

Section 5: Build Habitat Model¶

Train XGBoost¶

Create Model Evaluation Function¶

Evaluate XGBoost Model¶

Section 5.5: Exercise - Build a Random Forest Model¶

Background¶

Your Task¶

Section 6: Model Interpretation with SHAP¶

Section 7: Habitat Suitability Predictions¶

Section 8: Environment Condition Change Analysis¶

`starec.csv` – Station (tow) metadata¶

`bgsrec.csv` – Biological catch records¶

`envrec.csv` – Environmental measurements at stations¶