1.6.3. Validation of ocean colour observations for harmful algal bloom application

Production Date: 15-11-2025
Produced by: Chiara Volta (ENEA, Italy)

🌍 Use case: Detection of harmful algal bloom

❓ Quality assessment question

How well do ocean colour observations detect Karenia brevis harmful bloom on the West Florida Shelf?

The Ocean Colour dataset version 6.0, distributed by the Copernicus Climate Change Service (C3S), includes two variables: the chlorophyll-a mass concentration (chlor_a) and the remote sensing reflectance (Rrs) at six visible wavelengths, from October 1997 to the present [1].
Ocean colour observations have long supported the detection of algal blooms, enabling broader and more frequent coverage than in situ sampling [2] [3]. Moreover, as harmful algal blooms (HABs) are expected to increase in frequency, intensity and duration due to both climatic and non-climatic drivers (e.g., increasing temperature, higher nutrient inputs from river runoff) [4], these observations are increasingly valuable for environmental monitoring and decision-making.
Here, the goal is to evaluate how accurately the Ocean Colour products can identify the HABs of Karenia brevis (K. brevis), a toxic dinoflagellate responsible for severe ecological and economic damage on the West Florida Shelf (WFS) [5] [6] [7].

📢 Quality assessment statement

These are the key outcomes of this assessment

• The Ocean Colour dataset version 6.0 enables the calculation of two satellite-derived indicators for Karenia brevis detection on the West Florida Shelf: chlorophyll anomaly (chl_anom) and particulate backscattering ratio at 560 nm (bbp_560_ratio).

• The satellite-derived indicators maintain stable performance from 2000 to 2024, with an average total accuracy above 63% relative to in situ observations, confirming the dataset as a robust tool for long-term monitoring.

• The two indicators show different sensitivities to bloom conditions: chl_anom provides the highest HAB detection accuracy, whereas bbp_560_ratio is more susceptible to interference from non-algal particles in coastal waters.

• While both satellite-derived indicators generally agree, particularly in identifying non-bloom conditions, disagreements occur when optical signals are weak or optical coastal complexity is high.

• Combining the two satellite-derived indicators reduces false positive blooms and improves the overall accuracy, but increases missed detections (false negative) by excluding bloom identified by only a single indicator, illustrating the limitations of fixed thresholds.

📋 Methodology

This notebook assesses of the ability of the Ocean Colour dataset version 6.0 to identify the HABs of K. brevis on the WFS over a 25-year period (2000-2024). Following the approaches described in [8] and [9], chlor_a and Rrs at 560 nm data are used to derive two satellite-derived indicators: chlorophyll-a concentration anomaly (chl_anom) and total particulate backscattering ratio at 560 nm (bbp_560_ratio). To validate these indicators, satellite pixels are compared with in situ K. brevis cell concentrations obtained from the Florida Fish and Wildlife Conservation Commission (FWC) Open Data portal, allowing each pixel to be classified as either HAB or non-HAB. Accuracy metrics are calculated for each indicator individually to determine their sensitivity to bloom conditions, as well as for a combined approach in which a HAB is flagged only when both indicators are simultaneously triggered. An inter-algorithm analysis is also performed to quantify the agreement/disagreement between the indicators in detecting the same biological signal. In addition, the stabilty of the detection performance of the two indicators is examined across the 25-year period to evaluate the consistency of the dataset over time. Finally, the spatial patterns of the indicators are analysed to evaluate how the optical complexity of the region influences the detection accuracy of each indicator.

The analysis and results are detailed in the sections below:

1. Set up the code and choose the satellite data to use

• Import libraries

• Setup of satellite data request

2. Data retrieval and filtering

• Retrieval and filtering of satellite data

• Download and extraction of in situ K. brevis cell concentrations

3. Satellite and in situ data pairing, and calculation of satellite-derived HAB indicators

• Satellite-in situ data match-up

• Computation of satellite-derived indicators of K. brevis HABs

4. Accuracy computation and results organisation

• Accuracy calculation

• Results organisation for visualisation

5. Display and discuss results

• Display results

• Discussion

📈 Analysis and results

1. Set up the code and choose the satellite data to use

Import libraries

Besides the standard libraries used for managing and visualising multidimensional arrays, the C3S EQC custom function c3s_eqc_automatic_quality_control is imported to retrieve satellite data. The requests module is used to download and manage observed K. brevis cell concentrations from the FWC website, while the confusion_matrix is used to quantify the accuracy of satellite-derived HAB/non-HAB flags against in situ observations, and to evaluate agreement between the two satellite-derived indicators.

Show code cell source

Hide code cell source

import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from IPython.display import display, HTML
from c3s_eqc_automatic_quality_control import download, utils
import requests
from sklearn.metrics import confusion_matrix
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)

Setup of satellite data request

The analysis in this notebook uses chlor_a and Rrs_560 nm data acquired over the WFS over 25 years (2000-2024). The data request is configured to also retrieve data from 1999 in order to provide the two-month baseline required to calculate the chl_anom, which is defined relative to a period ending two weeks before each in situ observed HAB event (see ‘Computation of satellite-derived indicators of K. brevis HABs’ section below).

Show code cell source

Hide code cell source

region_map = {
    "custom_region": {
        "lon_slice": slice(-85.5, -81.5),
        "lat_slice": slice(27.5, 30.5),
    }
}
start = "1999-01"
stop = "2024-12"
variables = ["chlor_a", "Rrs_560"]
collection_id = "satellite-ocean-colour"
request = {
    "projection": "regular_latitude_longitude_grid",
    "version": "6_0",
    "format": "zip",
}

2. Data retrieval and filtering

Retrieval and filtering of satellite data

The regionalised_filtred_data function is defined and applied to retrieve chlor_a and Rrs_560 data from the region of interest. Only chlor_a values within the valid range of 0.01–100 mg m^-3 are downloaded [10]. All satellite data are merged in a single array.

Show code cell source

Hide code cell source

def regionalised_filtred_data(
    ds, variable, lon_slice, lat_slice, vmin, vmax
):
    da = ds[variable]
    da = utils.regionalise(da, lon_slice=lon_slice, lat_slice=lat_slice)    
    if vmin is not None:
        da = da.where(da > vmin)
    if vmax is not None:
        da = da.where(da < vmax)        
    return da.to_dataset(name=variable)
datasets = []
for variable in variables:
    for region, slices in region_map.items():
        current_request = download.update_request_date(
            request | {
                "variable": "remote_sensing_reflectance" 
                if variable.startswith("Rrs") 
                else "mass_concentration_of_chlorophyll_a"
            },
            start=start,
            stop=stop,
            stringify_dates=True,
        )        
        ds_var = download.download_and_transform(
            collection_id,
            current_request,
            transform_func=regionalised_filtred_data,
            transform_func_kwargs=slices | {
                "variable": variable,
                "vmin": 1.0e-2 if variable == "chlor_a" else None,
                "vmax": 1.0e2 if variable == "chlor_a" else None,
            },
            chunks={"year": 1, "month": 1}, 
            quiet=True,
        )
        datasets.append(ds_var.expand_dims(region=[region]))
        
ds_sat = xr.merge(datasets).compute()

Download and extraction of in situ K. brevis cell concentrations

In situ data are retrieved from the FWC Open Data portal, which provides laboratory-analysed K. brevis cell concentrations from offshore and inshore water samples collected since 1953. As sampling efforts increased substantially after the year 2000, only data from 2000 to 2024 are downloaded and used in the analysis. Additionally, downloaded data are restricted to samples collected over the WFS at depths shallower than 2 m. Despite the progressive increase in the number of samples, observations remain sparse in time and space and are typically collected during K. brevis bloom periods (from late summer through winter [9]) or in response to severe HAB events. Additional details on sampling and monitoring protocols are available on the FWC HAB monitoring webpage (see “Further readings” below).

Show code cell source

Hide code cell source

datasets = [
    {"name": "2000-2006", "layer": "5", "is_recent": False},
    {"name": "2007-2014", "layer": "6", "is_recent": False},
    {"name": "2015-2023", "layer": "12", "is_recent": False},
    {"name": "2024", "layer": "7", "is_recent": True},
]
all_dfs = []
for ds in datasets:
    if ds["is_recent"]:
        url = f"https://gis.myfwc.com/mapping/rest/services/Open_Data/Recent_Harmful_Algal_Bloom__HAB__Events_2024_present/MapServer/{ds['layer']}/query"
    else:
        base_path = "https://gis.myfwc.com/mapping/rest/services/Open_Data/Historic_Harmful_Algal_Bloom_Events_"
        url = f"{base_path}{ds['name'].replace('-', '___')}/MapServer/{ds['layer']}/query"
    all_records = []
    offset = 0
    step = 2000
    while True:
        params = {
            "where": "1=1",
            "outFields": "*",
            "f": "json",
            "resultOffset": offset,
            "resultRecordCount": step,
            "returnGeometry": "false",
        }
        try:
            r = requests.get(url, params=params, timeout=30)
            r.raise_for_status()
            data = r.json()
            feats = data.get("features", [])
            if not feats:
                break
            all_records.extend([f["attributes"] for f in feats])
            offset += step
        except Exception:
            break
    temp_df = pd.DataFrame(all_records)
    if not temp_df.empty:
        date_col = next(
            (col for col in temp_df.columns if "date" in col.lower()), None
        )
        if date_col:
            temp_df = temp_df.rename(columns={date_col: "SAMPLE_DATE"})
            temp_df["SAMPLE_DATE"] = pd.to_datetime(
                temp_df["SAMPLE_DATE"], unit="ms", errors="coerce"
            )
        if ds["is_recent"] and "SAMPLE_DATE" in temp_df.columns:
            temp_df = temp_df[temp_df["SAMPLE_DATE"] <= "2024-12-31 23:59:59"]
        all_dfs.append(temp_df.assign(DATASET_SOURCE=ds["name"]))
df = pd.concat(all_dfs, ignore_index=True) if all_dfs else pd.DataFrame()

# Data Subsetting and Filtering
lat_min, lat_max = 27.5, 30.5
lon_min, lon_max = -85.5, -81.5
if not df.empty:
    df["LATITUDE"] = pd.to_numeric(df["LATITUDE"], errors="coerce")
    df["LONGITUDE"] = pd.to_numeric(df["LONGITUDE"], errors="coerce")
    df["DEPTH"] = pd.to_numeric(df["DEPTH"], errors="coerce")
    df["COUNT_"] = pd.to_numeric(df["COUNT_"], errors="coerce")
    if "NAME" in df.columns and df["NAME"].dtype == "object":
        df["NAME"] = df["NAME"].str.strip()
    mask = (
        (df["LATITUDE"].between(lat_min, lat_max))
        & (df["LONGITUDE"].between(lon_min, lon_max))
        & (df["DEPTH"] <= 2)
        & (df["COUNT_"] > 0)
        & (df["NAME"] == "Karenia brevis")
    )
    df_filt = df.loc[mask].copy()
else:
    df_filt = pd.DataFrame()

3. Satellite and in situ data pairing, and calculation of satellite-derived HAB indicators

Satellite-in situ data match-up

The extract_hab_pairs_final function is defined to resolve the mismatch in scale between discrete in situ samples and continuous satellite imagery, ensuring that the extracted satellite data accurately represent the sampled K. brevis concentrations. First, a spatio-temporal binning is applied by averaging all in situ samples collected within the same 0.04°x0.04° grid cell on a single day to produce a single representative data point. This step ensures that the biological ground-truth corresponds to the spatial resolution of the Ocean Colour product ([1] [11]). Next, the match-up procedure pairs the nearest satellite pixel within a maximum distance of 4 km from each averaged in situ location, ensuring that each binned in situ observation is paired with its closest satellite data. The pairing is further constrained to a \(\pm7\)-day temporal window, following [11]. This window is ecologically relevant given the typical multi-day persistence and slow-growing nature of K. brevis blooms ([12][13]), thereby minimizing the impact of water-mass movement and/or bloom evolution. The resulting dataset contains only the coupled satellite–in situ pairs and is used for subsequent calculations of HAB indicators and accuracy.

Show code cell source

Hide code cell source

def extract_hab_pairs_final(dataset, df_input, max_dist_km=4.0, max_days=7):
    df_work = df_input.copy()
    df_work['SAMPLE_DATE'] = pd.to_datetime(df_work['SAMPLE_DATE']).dt.normalize()
    df_work['lat_bin'] = df_work['LATITUDE'].round(2)
    df_work['lon_bin'] = df_work['LONGITUDE'].round(2)    
    df_grouped = df_work.groupby(['SAMPLE_DATE', 'lat_bin', 'lon_bin']).agg({
        'LATITUDE': 'mean',
        'LONGITUDE': 'mean',
        'COUNT_': 'mean'
    }).reset_index()
    extracted = []    
    for _, row in df_grouped.iterrows():
        sample_time = np.datetime64(row["SAMPLE_DATE"])
        time_diff = np.abs(dataset.time.values - sample_time)
        valid_window = dataset.isel(time=(time_diff / np.timedelta64(1, 'D')) <= max_days)        
        if len(valid_window.time) == 0:
            continue            
        # Selection of the nearest pixel
        nearest_pixel = valid_window.sel(
            latitude=row["LATITUDE"], 
            longitude=row["LONGITUDE"], 
            method="nearest"
        ).isel(region=0)        
        # Calculate actual distance to center of nearest pixel
        lat_dist = abs(nearest_pixel.latitude.values - row["LATITUDE"]) * 111
        lon_dist = abs(nearest_pixel.longitude.values - row["LONGITUDE"]) * (111 * np.cos(np.radians(row["LATITUDE"])))
        total_dist = np.sqrt(lat_dist**2 + lon_dist**2)        
        if total_dist > max_dist_km:
            continue            
        clear_days = nearest_pixel.dropna(dim='time', how='any', subset=['chlor_a', 'Rrs_560'])        
        if len(clear_days.time) == 0:
            continue             
        best_match = clear_days.sortby(np.abs(clear_days.time - sample_time)).isel(time=0)
        extracted.append({
            "sample_date": row["SAMPLE_DATE"],
            "sat_date": pd.to_datetime(best_match.time.values),
            "lat": row["LATITUDE"],
            "lon": row["LONGITUDE"],
            "avg_in_situ_count": row["COUNT_"],
            "sat_chlor_a": float(best_match.chlor_a),
            "sat_Rrs_560": float(best_match.Rrs_560),
            "dist_km": round(total_dist, 2),
            "days_diff": round(abs((best_match.time.values - sample_time) / np.timedelta64(1, 'D')), 2)
        })
    return pd.DataFrame(extracted)
    
matched_df = extract_hab_pairs_final(ds_sat, df_filt)

Computation of satellite-derived indicators of K. brevis HABs

Two satellite-derived K. brevis HAB indicators are computed for each satellite-in situ pair: chl_anom and bbp_560_ratio.
The function get_dynamic_baseline is defined to calculate the mean chlor_a from the nearest satellite pixel to the in situ sampling location over a 2-month period ending 2 weeks before the sampling date, and the chl_anom is calculated by subtracting this baseline from the chlor_a observed on the date of the pair, following the method by [8], then adopted by [9] and [15]. Note that chlor_a averages are computed using log-transformed daily data, and then unlogged [14].
The bbp_560_ratio is computed by adapting the equation for bbp ratio at 555 nm from [9]. This approach relies on the same assumption made by [9], namely that the difference between Rrs at 555 nm and Rrs at 560 nm is negligible, and it is further justified because during K. brevis HABs, Rrs strongly decreases across all wavelengths below 600 nm [15]. Ultimately, this choice is operationally motivated by data availability, as neither Rrs at 555 nm nor at 550 nm used by [9] and [15], respectively, is available in the Ocean Colour dataset. The bbp_560_ratio is calculated as: \((-0.00182+2.058⋅Rrs\_560)/((0.3⋅chl^{0.62})⋅(0.002+0.02⋅(0.5-0.25⋅log_{10}chl)))\), where Rrs_560 and chl are the satellite data corresponding to each satellite-in situ pair.

Show code cell source

Hide code cell source

def get_dynamic_baseline(row, full_dataset):
    # Define the 60-day window ending 14 days before the event
    sample_time = np.datetime64(row["sample_date"])
    end_baseline = sample_time - np.timedelta64(14, 'D')
    start_baseline = sample_time - np.timedelta64(74, 'D')        
    # Select the nearest spatial pixel from the xarray dataset
    try:
        pixel_ts = full_dataset.sel(
            latitude=row["lat"], 
            longitude=row["lon"], 
            method="nearest"
        ).isel(region=0)            
        # Slice the time range for the baseline window
        baseline_window = pixel_ts.sel(time=slice(start_baseline, end_baseline))            
        # Extract chlor_a values and remove NaNs
        chl_values = baseline_window.chlor_a.values
        chl_values = chl_values[~np.isnan(chl_values)]
        if len(chl_values) == 0:
            return np.nan
        chl_values = chl_values[chl_values > 0]
        if len(chl_values) == 0:
            return np.nan            
        log_mean = np.mean(np.log10(chl_values))
        mean_val = 10**log_mean         
        return mean_val    
    except Exception:
        return np.nan

# Calculate the moving baseline
matched_df['chl_baseline'] = matched_df.apply(lambda row: get_dynamic_baseline(row, ds_sat), axis=1)

# Calculate Chlorophyll Anomaly
matched_df['chl_anomaly'] = matched_df['sat_chlor_a'] - matched_df['chl_baseline']

# 3. Calculate bbp_560_ratio
chl = matched_df['sat_chlor_a']
rrs = matched_df['sat_Rrs_560']
numerator = -0.00182 + (2.058 * rrs)
denom_part1 = 0.3 * (chl**0.62)
denom_part2 = 0.002 + 0.02 * (0.5 - 0.25 * np.log10(chl))
matched_df['bbp_560_ratio'] = numerator / (denom_part1 * denom_part2)

final_analysis_df = matched_df.dropna(subset=['chl_baseline', 'bbp_560_ratio'])

4. Accuracy computation and results organisation

Accuracy calculation

Before the accuracy calculation, each satellite-in situ pair is classified as HAB or non-HAB using the thresholds adapted from [9] and [15]. An in situ sample is considered a K. brevis HAB if cell concentration is equal to or greater than 100000 cells/l. For both satellite-derived indicators, conservative thresholds are applied to optimise the balance between detection sensitivity and classification accuracy in the optically complex waters of the WFS. For the bbp_560_ratio, the stricter threshold from [15] is used to minimise coastal interference, and a match-up is considered a K. brevis HAB when the ratio is lower than 1. For the chl_anom, the threshold applied to detect a K. brevis HAB is increased from >1 mg m⁻³, used by [9], to >2 mg m⁻³ in order to focus on mature, high-concentration blooms while reducing the detection of false positive events. A combined indicator is also created, identifying a HAB only when both satellite-derived indicators detect it simultaneously.
Accuracy statistics are then calculated using the accuracy_metrics function. The function compares each satellite-derived indicator with in situ observations and extracts the number of true positives (TP) and true negatives (TN), which represent cases where satellite predictions and in situ observations agree; false positives (FP) when blooms are predicted but not observed in situ; and false negatives (FN) when blooms are observed in situ but not predicted. From these values, the function calculates the accuracy of the indicators (in %) in correctly identifying in situ blooms and non-blooms, following [9]. In particular, HAB accuracy is calculated as the proportion of correctly predicted blooms relative to the total number of HAB predicted by the satellite indicator (true and false), whereas non-HAB accuracy is the ratio between the number of correctly detected non-HAB cases and the total number of non-HAB cases identified by the satellite indicator. The overall accuracy, defined as the proportion of all correct predictions (positive and negative) relative to the total number of in situ observations, is also calculated.
Furthermore, a cross-comparison is performed to evaluate the consistency between the two satellite algorithms. For this analysis, chl_anom is treated as the reference against which the bbp_560_ratio is tested. Here, the agreement (in %) is defined as the proportion of the total observations where both methods simultaneously identify a HAB (positive agreement) or a non-HAB (negative agreement). Conversely, the disagreement represents cases where the two indicators yield conflicting classifications, whereas the overall agreement is calculated as the proportion of all co-occurring detections (positive and negative) relative to the total number of observations.

Show code cell source

Hide code cell source

final_analysis_df = final_analysis_df.copy()

THRESHOLDS = {
    "cell_count": 100000,   # in situ bloom threshold
    "chl_anomaly": 2.0,     # satellite chlorophyll anomaly threshold
    "bbp_ratio": 1.0        # satellite bbp_560_ratio threshold
}

# Classify bloom/non-bloom for each pair
final_analysis_df["obs_bloom"] = (final_analysis_df["avg_in_situ_count"] >= THRESHOLDS["cell_count"]).astype(int)
final_analysis_df["sat_bloom_chl"] = (final_analysis_df["chl_anomaly"] > THRESHOLDS["chl_anomaly"]).astype(int)
final_analysis_df["sat_bloom_bbp"] = (final_analysis_df["bbp_560_ratio"] < THRESHOLDS["bbp_ratio"]).astype(int)

# Create combined indicator (both satellite indicators detect bloom)
final_analysis_df['sat_combined'] = ((final_analysis_df['sat_bloom_chl'] == 1) & 
                                     (final_analysis_df['sat_bloom_bbp'] == 1)).astype(int)

# Accuracy metrics function
def accuracy_metrics(actual, predicted):
    cm = confusion_matrix(actual, predicted, labels=[0, 1])
    tn, fp, fn, tp = cm.ravel()
    ua_b = tp / (tp + fp) if (tp + fp) > 0 else np.nan  # HAB user accuracy
    ua_nb = tn / (tn + fn) if (tn + fn) > 0 else np.nan  # non-HAB user accuracy
    oa = (tp + tn) / len(actual)                          # overall accuracy
    return tp, tn, fp, fn, ua_b, ua_nb, oa

# Compute accuracy metrics for each indicator
tp_c, tn_c, fp_c, fn_c, ua_bc, ua_nbc, oac = accuracy_metrics(final_analysis_df["obs_bloom"], final_analysis_df["sat_bloom_chl"])
tp_b, tn_b, fp_b, fn_b, ua_bb, ua_nbb, oab = accuracy_metrics(final_analysis_df["obs_bloom"], final_analysis_df["sat_bloom_bbp"])
tp_comb, tn_comb, fp_comb, fn_comb, ua_b_comb, ua_nb_comb, oa_comb = accuracy_metrics(final_analysis_df["obs_bloom"], final_analysis_df["sat_combined"])

# Calculate the Inter-Algorithm Confusion Matrix
cm_inter = confusion_matrix(final_analysis_df["sat_bloom_chl"], 
                            final_analysis_df["sat_bloom_bbp"], labels=[0, 1])
tn_i, fp_i, fn_i, tp_i = cm_inter.ravel()

Result organisation for visualisation

Accuracy metrics are organized into three dataframes to facilitate their visualisation (see ‘Display results’ section below): one for the counts of correct and incorrect detections (TP, TN, FP, FN) across the chl_anom, bbp_560_ratio, and combined indicators; one for the derived accuracies in identifying HAB and non-HAB events; and one for the agreement and disagreement rates between the two satellite algorithms. Additionally, each satellite-in situ pair is labelled by its classification category to enable the spatial mapping of detection performance for each satellite indicator.

Show code cell source

Hide code cell source

# Table 1 Data
total_n = len(final_analysis_df)
rows_cm = ["True Positive (TP)", "True Negative (TN)", "False Positive (FP)", "False Negative (FN)", 
           "Correctly Classified", "Incorrectly Classified", "Total"]
data_cm = [
    [tp_c, tp_b, tp_comb], [tn_c, tn_b, tn_comb], 
    [fp_c, fp_b, fp_comb], [fn_c, fn_b, fn_comb],
    [(tp_c + tn_c), (tp_b + tn_b), (tp_comb + tn_comb)],
    [(fp_c + fn_c), (fp_b + fn_b), (fp_comb + fn_comb)],
    [total_n]*3
]

# Table 2 Data
acc_data = [
    {"Indicator": "chl_anom", "HAB Accuracy (%)": ua_bc, "non-HAB Accuracy (%)": ua_nbc, "Total Accuracy (%)": oac},
    {"Indicator": "bbp_560_ratio", "HAB Accuracy (%)": ua_bb, "non-HAB Accuracy (%)": ua_nbb, "Total Accuracy (%)": oab},
    {"Indicator": "Combined (AND)", "HAB Accuracy (%)": ua_b_comb, "non-HAB Accuracy (%)": ua_nb_comb, "Total Accuracy (%)": oa_comb}
]
accuracy_df = pd.DataFrame(acc_data).set_index("Indicator")

# Table 3 Data
inter_data = {
    "Agreement type": [
    "Agreement (HAB)",
    "Agreement (no-HAB)",
    "Disagreement (chl_anom: HAB, bbp_560_ratio: no-HAB)",
    "Disagreement (chl_anom: no-HAB, bbp_560_ratio: HAB)",
    "Overall Agreement"],
    "Count": [tp_i, tn_i, fn_i, fp_i, (tp_i + tn_i)],
    "Percentage (%)": [(tp_i/total_n)*100, (tn_i/total_n)*100, (fn_i/total_n)*100, (fp_i/total_n)*100, (tp_i + tn_i) / total_n * 100]
}
df_inter = pd.DataFrame(inter_data)

# Create confusion labels for mapping
for col, pred_col in [("chl", "sat_bloom_chl"), ("bbp", "sat_bloom_bbp"), ("combined", "sat_combined")]:
    actual = final_analysis_df["obs_bloom"]
    pred = final_analysis_df[pred_col]   
    conditions = [
        (actual == 1) & (pred == 1), # TP
        (actual == 0) & (pred == 0), # TN
        (actual == 0) & (pred == 1), # FP
        (actual == 1) & (pred == 0)  # FN
    ]
    choices = ["TP", "TN", "FP", "FN"]
    final_analysis_df[f"sat_bloom_{col}_confusion"] = np.select(conditions, choices, default="Unknown")

5. Display and discuss results

Display results

Table 1 summarises the counts of true positives (TP), true negatives (TN), false positive (FP) and false negative (FN), as well as the total correctly and incorrectly classified observations detected by the chl_anom, bb_560_ratio and their combination (i.e., both indicators have to be simultaneously satisfied; see ‘Accuracy calculation’ section above), as well as the total number of paired satellite-in situ observations. Table 2 presents the accuracy metrics derived from these counts, reporting the performance of each satellite indicator and their combination as the percentage of correctly identified HABs and non-HABs, along with total accuracy. Table 3 details the inter-algorithm agreement, showing the counts and percentages of concordant and discordant detections between the two satellite indicators. Annual time series plots for each indicator, displaying the annual HAB, non-HAB, and total accuracy alongside the total number of satellite-in situ match-ups (grey bars), are also provided to evaluate their temporal consistency. These figures also display the number of in situ HAB events exceeding the 100000 cell/l threshold in parentheses. Finally, two maps are displayed to visually assess the spatial performance of each indicator. Here, coloured markers are used to indicate TP (green), FP (orange), TN (cyan) and FN (red).

Show code cell source

Hide code cell source

# Display Table 1
columns_cm = pd.MultiIndex.from_tuples([
    ("Satellite-derived indicator (counts):", "chl_anom"),
    ("Satellite-derived indicator (counts):", "bbp_560_ratio"),
    ("Satellite-derived indicator (counts):", "Combined (AND)")])
df_cm = pd.DataFrame(data_cm, index=rows_cm, columns=columns_cm)
display(HTML("<h3>Table 1: Satellite-in situ match-up counts</h3>"))
display(df_cm.style.set_table_styles([
    {'selector': 'th', 'props': [('text-align', 'center'), ('background-color', '2px solid #ccc')]},
    {'selector': 'tr:nth-child(5)', 'props': [('text-align', 'center'), ('border-top', '2px solid #ccc')]},
    {'selector': 'tr:nth-child(7)', 'props': [('text-align', 'center'), ('border-top', '2px solid #ccc')]}
]))

# Display Table 2
disp_acc = accuracy_df.copy()
for col in disp_acc.columns:
    disp_acc[col] = disp_acc[col].apply(lambda x: f"{round(x*100, 1)}%" if not pd.isna(x) else "NaN")
display(HTML("<h3 style='margin-top: 50px; margin-bottom: 10px;'>Table 2: Accuracy Percentage Comparison</h3>"))
display(disp_acc.style.set_table_styles([
    {'selector': 'th', 'props': [('background-color', '#4CAF50'), ('color', 'white'), ('text-align', 'center')]},
    {'selector': 'td', 'props': [('text-align', 'center'), ('padding', '8px')]}
]))

# Display Table 3
display(HTML("<h3 style='margin-top: 50px; margin-bottom: 10px;'>Table 3: Inter-Algorithm Agreement (chl_anom vs. bbp_560_ratio)</h3>"))
display(df_inter.style.format({"Percentage (%)": "{:.1f}%"}).hide(axis="index").set_table_styles([
    {'selector': 'th', 'props': [('text-align', 'center'), ('background-color', '#01579B'), ('color', 'white')]},
    {'selector': 'tr:nth-child(5)', 'props': [('text-align', 'center'), ('border-top', '2px solid #333'), ('font-weight', 'bold')]},
    {'selector': 'td', 'props': [('text-align', 'center'), ('padding', '10px')]}
]))

# Display time series
display(HTML("<br><br>"))
fig, axes = plt.subplots(2, 1, figsize=(16, 12), sharex=True)
indicators = {
    'sat_bloom_chl': {'label': 'chl_anom', 'color': 'green'},
    'sat_bloom_bbp': {'label': 'bbp_560_ratio', 'color': 'green'}
}
final_analysis_df['year'] = pd.to_datetime(final_analysis_df['sample_date']).dt.year
for i, (col_name, info) in enumerate(indicators.items()):
    ax1 = axes[i]
    yearly_stats = []        
    for year, group in final_analysis_df.groupby('year'):
        valid_group = group.dropna(subset=[col_name, "obs_bloom"])        
        if not valid_group.empty:
            cm = confusion_matrix(valid_group["obs_bloom"], valid_group[col_name], labels=[0, 1])
            tn, fp, fn, tp = cm.ravel()                       
            ua_b = (tp / (tp + fp) * 100) if (tp + fp) > 0 else np.nan
            ua_nb = (tn / (tn + fn) * 100) if (tn + fn) > 0 else np.nan
            oa = ((tp + tn) / len(valid_group) * 100)            
            yearly_stats.append({
                'Year': year,
                'HAB_Accuracy': ua_b,
                'non_HAB_Accuracy': ua_nb,
                'Total_Accuracy': oa,
                'Sample_Count': len(valid_group),
                'HAB_Count': tp 
            })               
    y_df = pd.DataFrame(yearly_stats)    
    line1, = ax1.plot(y_df['Year'], y_df['HAB_Accuracy'], marker='o', color=info['color'], linewidth=2.5, label='HAB Accuracy (%)')
    line2, = ax1.plot(y_df['Year'], y_df['non_HAB_Accuracy'], marker='^', color='orange', linewidth=2, linestyle='--', label='non-HAB Accuracy (%)')
    line3, = ax1.plot(y_df['Year'], y_df['Total_Accuracy'], marker='s', color='blue', 
                      linewidth=1.5, alpha=0.7, label='Total Accuracy (%)')    
    ax1.set_ylabel('Accuracy (%)', fontsize=15, fontweight='bold')
    ax1.set_ylim(-5, 125)
    ax1.grid(True, alpha=0.2)
    ax1.set_title(f'Detection performance of {info["label"]}', fontsize=16, loc='left', pad=15, fontweight='bold')    
    ax2 = ax1.twinx()
    ax1.tick_params(axis='both', labelsize=13)
    ax2.tick_params(axis='y', labelsize=13)    
    bars = ax2.bar(y_df['Year'], y_df['Sample_Count'], alpha=0.15, color='gray', label='Total Match-ups (of which in situ HABs)')
    ax2.set_ylabel('Number of Match-ups', fontsize=15, fontweight='bold')
    for idx, bar in enumerate(bars):
        height = bar.get_height()
        hab_val = y_df['HAB_Count'].iloc[idx]
        ax2.text(bar.get_x() + bar.get_width()/2., height + 0.5,
                 f'{int(height)}\n({int(hab_val)})', 
                 ha='center', va='bottom', fontsize=10, color='black', fontweight='bold')        
    lines, labels = ax1.get_legend_handles_labels()
    bars_h, bars_l = ax2.get_legend_handles_labels()
axes[1].set_xlabel('Year', fontsize=18, fontweight='bold', labelpad=15)
fig.legend(lines + bars_h, labels + bars_l, 
           loc='lower center', 
           ncol=4, 
           fontsize=14, 
           frameon=False,
           bbox_to_anchor=(0.5, 0.05)) 
plt.tight_layout(rect=[0, 0.12, 1, 0.95])
plt.show()

# Display Maps
display(HTML("<br><br>"))
conf_colors = {"TP": "green", "FP": "orange", "TN": "cyan", "FN": "red"}
conf_markers = {"TP": "o", "FP": "s", "TN": "d", "FN": "v"}
fig, axes = plt.subplots(1, 2, figsize=(14, 8), subplot_kw={'projection': ccrs.PlateCarree()})
for i, (alg_col, title) in enumerate([("sat_bloom_chl_confusion", "chl_anom"), 
                                      ("sat_bloom_bbp_confusion", "bbp_560_ratio")]):
    ax = axes[i]
    ax.add_feature(cfeature.LAND, facecolor="lightgray")
    ax.coastlines(resolution="10m")
    ax.set_extent([-86, -81, 26, 31], crs=ccrs.PlateCarree())    
    for label in ["TP", "FP", "TN", "FN"]:
        subset = final_analysis_df[final_analysis_df[alg_col] == label]
        ax.scatter(subset["lon"], subset["lat"], 
                   color=conf_colors[label], marker=conf_markers[label],
                   label=label, transform=ccrs.PlateCarree(),
                   edgecolor="black", s=60, zorder=5)
    ax.set_title(f"{title} classification", fontsize=14)
fig.legend(["True Positive (TP)", "False Positive (FP)", "True Negative (TN)", "False Negative (FN)"], loc='lower center', ncol=4, bbox_to_anchor=(0.5, 0.1))
plt.show()

Table 1: Satellite-in situ match-up counts

	Satellite-derived indicator (counts):
	chl_anom	bbp_560_ratio	Combined (AND)
True Positive (TP)	468	352	188
True Negative (TN)	2200	2118	2527
False Positive (FP)	632	714	305
False Negative (FN)	566	682	846
Correctly Classified	2668	2470	2715
Incorrectly Classified	1198	1396	1151
Total	3866	3866	3866

Table 2: Accuracy Percentage Comparison

	HAB Accuracy (%)	non-HAB Accuracy (%)	Total Accuracy (%)
Indicator
chl_anom	42.5%	79.5%	69.0%
bbp_560_ratio	33.0%	75.6%	63.9%
Combined (AND)	38.1%	74.9%	70.2%

Table 3: Inter-Algorithm Agreement (chl_anom vs. bbp_560_ratio)

Agreement type	Count	Percentage (%)
Agreement (HAB)	493	12.8%
Agreement (no-HAB)	2193	56.7%
Disagreement (chl_anom: HAB, bbp_560_ratio: no-HAB)	607	15.7%
Disagreement (chl_anom: no-HAB, bbp_560_ratio: HAB)	573	14.8%
Overall Agreement	2686	69.5%

Discussion

Performance of individual satellite indicators

According to Table 1 and Table 2, chl_anom emerges as the most sensitive individual satellite indicator, achieving the highest HAB detection accuracy (42.5%) and a total accuracy of 69.0%. While it correctly identifies 468 true positive cases, the results reveal a significant number of false positives (632) and false negatives (566). These errors suggest that although chl_anom is a useful indicator for bloom conditions, the application of a fixed threshold applied over two decades may not fully account for the natural baseline variability in the region, often flagging non-HAB biomass or missing lower-density blooms. Moreover, following periods in which HABs persist in a location for several months, the dynamic anomaly calculation method (see ‘Computation of satellite-derived indicators of K. brevis HABs’ section above) would inflate the baseline 2-month mean once the bloom exceeds the 2-week offset window, progressively masking the bloom signal.
The bbp_560_ratio shows relatively lower overall performance, with a total accuracy of 63.9%. It also has the lowest HAB detection accuracy among the tested methods (33.0%) and produces the highest number of both false positives and false negatives (714 and 682, respectively). This indicates that the optical signature of bbp_560_ratio, intended to highlight the low-backscattering nature of K. brevis [9], is frequently confounded, leading both to more false alarms and missed detections when the K. brevis signal is not optically dominant.

Temporal stability of detection performance

As shown in the annual performance time series, the total accuracy of both satellite-derived indicators remains rather stable across the 25-year period analysed, despite significant inter-annual variability in the number of available match-ups and in situ HABs events. This indicates a stable temporal behaviour in the Ocean Colour dataset, which provides consistent results despite the use of successive generations of satellite sensors [1], and suggests that the dataset is well suited for long-term monitoring in the region.
The analysis also indicates that chl_anom acts as a more conservative indicator than bbp_560_ratio. Indeed, the time series show that during certain years (e.g., 2007, 2019), the chl_anom remains below the detection threshold (2 mg m^-3), while the bbp_560_ratio continues to trigger. This demonstrates that chl_anom effectively prioritises high-confident detection, as also indicated by the fewer false alarms produced compared to the bbp_560_ratio (Table 1).

Agreement between satellite indicators

To determine whether the two different physical signals (i.e., chl concentration and particle backscattering) converge on the same biological signal over the long term, the agreement between the two satellite indicators is examined. Results summarised in Table 3 indicate an overall agreement of 69.5%, which suggests that the algorithms largely detect the same bloom-related optical signals across different sensor technologies and varying environmental conditions. Most of this agreement is driven by the dominance of clear-water conditions (56.7%), confirming the ability of the dataset in identifying non-HAB conditions (Tables 1 and 2). The active bloom agreement (12.8%) highlights the most robust detections, where the K. brevis signal is strong enough to trigger both indicators simultaneously. These cases therefore represent the highest-confidence satellite detections in the records. Results also indicate that the two indicators diverge in the 30.5% of the cases. This decoupling likely reflects multiple factors, including biological processes (e.g., vertical migration [11]) and/or the presence of optically active constituents (e.g., suspended matter, mixed phytoplankton assemblages, coloured dissolved organic matter (CDOM) [9] [15]). Additionally, a threshold effect likely contributes to this disagreement. Indeed, when a bloom falls below the detection limit of one indicator, the other may still detect it. These decoupled cases highlight the complementary nature of the two indicators and suggest that relying on a single indicator would result in missing up to 16% of the potential HABs detected by the alternative method.

Spatial patterns of detection performance

The variations in algorithm performance are closely linked to the geographical characteristics of the study area. As shown in the maps, true negatives are predominantly located offshore, whereas classification errors cluster along the coastline. This spatial distribution suggests that the optical complexity in nearshore waters, where the presence of material that can modify the pigment-to-backscattering ratio is higher (e.g., sediments, CDOM), exerts a strong influence on both indicators compared to clearer offshore regions. Additionally, satellite-in situ match-ups are more challenging in coastal environments owing to land adjacency effects on reflectance and the increased difficulty of atmospheric correction. While both indicators struggle with false negatives along the coast, the bbp_560_ratio exhibits a higher susceptibility to false positives in these areas. This further indicates that while bbp_560_ratio is sensitive to backscattering particles, it lacks the specificity of chl_anom in distinguishing K. brevis from other suspended solids in shallow, turbid waters. Ultimately, this spatial divergence also confirms that the disagreement between the two indicators is geographically driven, reinforcing the need for an integrated approach to address this optically complex region.

Implication of the combined detection approach

The optical interference and signal divergence observed in nearshore waters directly influence the performance of the “Combined (AND)” detection strategy (Tables 1 and 2). By requiring both indicators to trigger simultaneously, this strategy acts as a high-precision filter, achieving the highest total accuracy (70.2%) and reducing false positives to their lowest level (305). However, because it only accepts the 12.7% of cases where both algorithms agree on a HAB (Table 3), it inherently ignores the 30.5% of the data where the signals diverge. This results in the highest number of false negatives (846) and a HAB detection accuracy of only 38.1%. While effective for high-confidence confirmation, the “AND” logic overlooks the transitional or optically complex phases of a bloom where only one indicator is dominant. Interestingly, if a less conservative chl_anom threshold is applied following [9] (i.e., >1 mg m⁻³), the error rate of this indicator increases by about 4%, while marginally dropping the combined accuracy to 69.4% (data not shown). These results suggest that, although the choice of specific thresholds can influence the error dynamics, the overall stability of the combined approach is not compromised; instead, the threshold choice serves to selectively refine the error distribution.

Overall, this 25-year assessment shows that satellite indicators derived from the Ocean Colour dataset provide a valuable long-term monitoring tool for detecting K. brevis blooms on the West Florida Shelf. Ocean colour observations can detect K. brevis blooms with moderate accuracy, although performance varies between indicators and is reduced in optically complex nearshore waters. However, the reliance on fixed thresholds limits the ability of these indicators to capture the full spectrum of bloom conditions, and persistent, long-lasting blooms introduce mathematical constraints potentially underestimating bloom intensity. Future improvements may involve the use of multispectral and/or hyperspectral satellite data to develop species- or region-specific algorithms, and/or the application of machine-learning approaches that integrate multiple optical and environmental variables [16] [17] [18].

ℹ️ If you want to know more

Key resources

• CDS catalogue entry used in this notebook is the Ocean colour daily data from 1997 to present derived from satellite observations

• Data download is from CDS

• Product Documentation is available on the CDS website

Code libraries used:

• C3S EQC custom functions, c3s_eqc_automatic_quality_control, prepared by B-Open

Further readings:

• Ocean Colour Monitoring Portals: CMEMS, ESA-GlobColour Project, ESA-OC-CCI, ESA-Sentinel 3, IMOS, NASA, NOAA CoastWatch.

• Use cases: EUMETSAT, C3S.

• International Ocean Colour Coordinating Group (IOCCG).

• Harmful Algae Event Database (HAEDAT).

• Florida Fish and Wildlife Conservation Commission.

• Groom, S., et al. (2019). Satellite Ocean Colour: Current Status and Future Perspective. Frontiers in Marine Science, 6.

• IOCCG (2021). Observation of Harmful Algal Blooms with Ocean Colour Radiometry. Bernard, S., et al. (eds.), IOCCG Report Series, No. 20, International Ocean Colour Coordinating Group, Darthmouth, Canada.

References

[1] Jackson, T., et al. (2023). C3S Ocean Colour Version 6.0: Product User Guide and Specification. Issue 1.1. E.U. Copernicus Climate Change Service. Document ref. WP2-FDDP-2022-04_C3S2-Lot3_PUGS-of-v6.0-OceanColour-product.

[2] Stumpf, R.P., et al. (2007). Remote Sensing of Harmful Algal Blooms. In: Miller, R.L., et al. (eds.) Remote Sensing of Coastal Aquatic Environments. Remote Sensing and Digital Image Processing, vol. 7. Springer, Dordrecht, pp. 277-296.

[3] Caballero, C.B., et al. (2025). The need for advancing algal bloom forecasting using remote sensing and modeling: Progress and future directions. Ecological Indicators, 172.

[4] Bindoff, N.L., et al. (2019). Changing Ocean, Marine Ecosystems, and Dependent Communities. In: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate. H.-O. Pörtner, et al. (eds.). Cambridge University Press, Cambridge, UK and New York, NY, USA, pp. 447-587.

[5] Pierce, R.H., and Henry, M.S. (2008). Harmful algal toxins of the Florida red tide (Karenia Brevis): natural chemical stressors in South Florida coastal ecosystems. Ecotoxicology 17, pp. 623-631.

[6] Béchard, A. (2020). Economics losses to fishery and seafood related businesses during harmful algal blooms. Fisheries Research, 230.

[7] Alvarez, S., et al. (2024). Non-linear impacts of harmful algae blooms on the coastal tourism economy. Journal of Environmental Management, 351.

[8] Stumpf, R.P., et al. (2003). Monitoring Karenia brevis blooms in the Gulf of Mexico using satellite ocean colour imagery and other data. Harmful Algae, 2(2), pp. 147-160.

[9] Tomlinson, M.C., et al. (2009). An evaluation of remote sensing techniques for enhanced detection of the toxic dinoflagellate, Karenia brevis. Remote Sensing of Environment, 113(3), 598-609.

[10] Sathyendranath, S., et al. (2019). An Ocean-Colour time series for use in climate studies: the experience of the Ocean-Colour Climate Change Initiative (OC-CCI). Sensors, 19, 4285.

[11] Hu, C., et al. (2022). Karenia Brevis bloom patterns on the west Florida shelf between 2003 and 2019: Integration of Field and satellite observations. Harmful Algae, 17, 102289.

[12] Vargo, G.A. (2009). A brief summary of the physiology and ecology of Karenia brevis Davis (g. Hansen and Moestrup comb. nov.) red tides on the West Florida Shelf and of hypotheses posed for their initiation, growth, maintenance, and termination. Harmful Algae, 8(4), 573-584.

[13] Stumpf, R.P., et al. (2022). Quantifying Karenia brevis bloom severity and respiratory irritation impact along the shoreline of Southwest Florida. PLoS ONE, 17(1), e0260755.

[14] Campbell, J.W. (1995). The lognormal distribution as a model for bio-optical variability in the sea. Journal of Geophysical Research, 100, 13237-13254.

[15] Cannizzaro, J.P., et al. (2008). A novel technique for detection of the toxic dinoflagellate, Karenia brevis, in the Gulf of Mexico from remotely sensed ocean color data. Continental Shelf Research, 28(1), 137-158.

[16] Arias, F., et al. (2025). Mapping Harmful Algae Blooms: The Potential of Hyperspectral Imaging Technologies. Remote Sensing, 17(4), 608.

[17] Izadi, M., et al. (2021). A Remote Sensing and Machine Learning-Based Approach to Forecast the Onset of Harmful Algal Bloom. Remote Sensing, 13(19), 3863.

[18] Zahir, M. et al. (2024) A review on monitoring, forecasting, and early warning of harmful algal bloom. Aquaculture, 593, 741351.