import subprocess
import sys

COLAB = "google.colab" in sys.modules


def _install(package):
    if COLAB:
        ans = input(f"Install { package }? [y/n]:")
        if ans.lower() in ["y", "yes"]:
            subprocess.check_call(
                [sys.executable, "-m", "pip", "install", "--quiet", package]
            )
            print(f"{ package } installed!")


def _colab_install_missing_deps(deps):
    import importlib

    for dep in deps:
        if importlib.util.find_spec(dep) is None:
            if dep == "iris":
                dep = "scitools-iris"
            _install(dep)


deps = ["pyobis", "joblib"]
_colab_install_missing_deps(deps)

Investigating PACE and OBIS as a data source for phytoplankton observations¶

Created: 2026-08-25

The Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) and the Ocean Biodiversity Information System (OBIS) both provide information about phytoplankton distribution on a global scale. PACE Chlorophyll is the near-surface concentration of the green pigment, chlorophyll a, at the ocean surface. Found in tiny photosynthetic organisms – phytoplankton – chlorophyll is a measure of ocean ecosystem health (link). On the other hand, OBIS is the world’s largest global open-access infrastructure for marine biodiversity data, providing occurrence records at species levels.

The former will have a much higher global coverage while the former will certainly be more accurate in its data. While researches can ue one to calibrate the other or mix-and-match these data to create a better final product, we wonder how easy it is to fetch them and display on a common grid.

We will use the pyobis and earthaccess libraries. Earthaccess requires a login to get download access the data. Please navigate to https://earthaccess.readthedocs.io/en/latest/user/ and follow the instructions there on how to create one.

import earthaccess

auth = earthaccess.login(persist=True)

We can not define our search bounding box, time span, and filter for the 1 degree dataset only, the 4 km results will also be fetched. Because it can be a lot of data, we will avoid the 4km high resolution for now with a simple filter to the results and get only the files that end with 1deg.nc string.

results = earthaccess.search_data(
    short_name="PACE_OCI_L3M_CHL",
    temporal=("2025-01-01", "2025-12-31"),
    bounding_box=(-90.0, 40.0, -75.0, 47.0),
)

fileset = earthaccess.open(results)
one_degree = [f for f in fileset if f.full_name.endswith("1deg.nc")]

The next loop downloads and save the data. We will re-download the data if we do not have it already. The local cache will save us some time and ensure we always have the data even if the server is down.

Local caching is always a good practice to preserver the data and avoid putting a lot of strain on servers with multiple identical requests.

from pathlib import Path

import pandas as pd
import xarray as xr

fname = "chlor_a.nc"

if not Path(fname).exists():

    def fake_time(ds):
        time = ds.attrs["product_name"].split(".")[1].split("_")
        time = pd.to_datetime(time).mean()
        ds = ds.assign_coords(time=time)
        return ds["chlor_a"].expand_dims(time=1)

    ds = xr.open_mfdataset(one_degree, preprocess=fake_time)
    chlor_a = ds["chlor_a"].where(
        (ds["chlor_a"] < ds["chlor_a"].display_max)
        & (ds["chlor_a"] > ds["chlor_a"].display_min)
    )

    chlor_a = chlor_a.mean(dim="time").compute()
else:
    chlor_a = xr.open_dataset(fname)


# DataArray for the chlor_a variable only.
chlor_a = chlor_a["chlor_a"]

Now we can query OBIS for appropriate phytoplankton taxonomic groupings. Here we list some generic phytoplankton groupings which would be equivalent to the chlorophyll satellite proxy data. Luckily OBIS allows users to use the WoRMS taxonomic backbone as a way to search for more exapansize taxonomic ranks, like family.

term	WoRMS	OBIS Map
Diatoms (Bacillariophyceae)	https://www.marinespecies.org/aphia.php?p=taxdetails&id=148899	https://mapper.obis.org/?taxonid=148899
Dinoflagellates (Dinophyceae)	https://www.marinespecies.org/aphia.php?p=taxdetails&id=19542	https://mapper.obis.org/?taxonid=19542
Cyanobacteria (Blue-green algae)	https://www.marinespecies.org/aphia.php?p=taxdetails&id=146537	https://mapper.obis.org/?taxonid=146537
Coccolithophyceae	https://www.marinespecies.org/aphia.php?p=taxdetails&id=592906	https://mapper.obis.org/?taxonid=592906

Below we do an initial test to ensure our search criterias are working with pyobis.

from pyobis import occurrences

# IDs from the table above.
aphia_ids = ["148899", "19542", "146537", "592906"]

query = occurrences.search(taxonid=aphia_ids, size=100)
data = query.execute()

data

2026-08-26 11:10:04 - pyobis.cache.cache - INFO - Cache initialized at C:\Users\Mathew.Biddle\AppData\Local\pyobis
2026-08-26 11:10:04 - pyobis.obisutils - INFO - 100 to be fetched. Estimated time =0.064834473133087160 seconds
2026-08-26 11:10:04 - pyobis.obisutils - INFO - Fetching: [████████████████████████████████████████████████████████████████████████████████████████████████████] 100/100
2026-08-26 11:10:04 - pyobis.cache.cache - INFO - Cache initialized at C:\Users\Mathew.Biddle\AppData\Local\pyobis
2026-08-26 11:10:04 - pyobis.obisutils - INFO - Fetched 100 records.

pd.to_datetime(data.date_year, format="%Y").describe()

count                            91
mean     2014-05-21 06:35:36.263736
min             1972-01-01 00:00:00
25%             2013-01-01 00:00:00
50%             2016-01-01 00:00:00
75%             2019-01-01 00:00:00
max             2024-01-01 00:00:00
Name: date_year, dtype: object

The next step is to query for gridded occurrence data for the aphia_ids we selected above. We search for gridded occurrences because this reduces the amount of data sent across the web. We could collect specific point observations, but that would take a lot of bandwidth and processing power to ingest and convert to appropriate georeferenced data. The .grid() response also provides total occurrence counts per cell, which is what we want to collect anyways.

features = occurrences.grid(3, taxonid=aphia_ids).execute()

2026-08-26 11:12:20 - pyobis.cache.cache - INFO - Cache initialized at C:\Users\Mathew.Biddle\AppData\Local\pyobis
2026-08-26 11:12:20 - pyobis.cache.cache - INFO - Cache initialized at C:\Users\Mathew.Biddle\AppData\Local\pyobis

In order to match with PACE data later, we need to organize the data into geo-aware constructs. We will use geopandas.GeoDataFrames as the unified construct to align our data into.

import geopandas as gpd

gdf = gpd.GeoDataFrame.from_features(features)
gdf

Note that the occurence rages from 1 to very large numbers (~10e6). We need to log-normalize the colors for better visualization. After that, we can plot the two side.

import cartopy.crs as ccrs
import matplotlib.colors as colors
import matplotlib.pyplot as plt

norm = colors.LogNorm()

projection = ccrs.PlateCarree()
fig, (ax0, ax1) = plt.subplots(
    nrows=2,
    subplot_kw={"projection": projection},
    sharex=True,
    sharey=True,
    figsize=(11, 11),
)

ax = gdf.plot(
    column="n",
    norm=colors.LogNorm(),
    legend=True,
    legend_kwds={"label": "OBIS Occurrences"},
    ax=ax0,
)


cs = chlor_a.plot(norm=colors.LogNorm(), ax=ax1)

We can see some features, like a higher density near the coast, currents, and the oligotrophic areas at the gyres centers. As expected, OBIS data density falls off quickly outside of the North Atlantic. However, there is a clear “Gulf Stream signal” in the dataset that is absent in the PACE data. We can envison some sort of interpolation that will use the OBIS data to increase precision and PACE to augument the spatial resolution.

We won’t be doing any of that here though. Our goal is to demostrate how to fetch matching data and display them side-by-side. In order to test a better “qualitative” way to compare the data we will average out the PACE data to the OBIS occurrences geometry, and put the data on an hexagonal hierarchical geospatial indexing system (H3) for plotting and comparing each hexagon.

def average_chlor_a(geom):
    xmin, ymin, xmax, ymax = geom

    c = chlor_a.where(
        (xmin <= chlor_a["lon"])
        & (chlor_a["lon"] <= xmax)
        & (ymin <= chlor_a["lat"])
        & (chlor_a["lat"] <= ymax),
        drop=True,
    )
    return float(c.mean().to_numpy())

We will use joblib’s Parallel to speed up the averages.

import multiprocessing

from joblib import Parallel, delayed

num_cores = multiprocessing.cpu_count()

bounds = gdf["geometry"].bounds.to_numpy()

chlor_avg = Parallel(n_jobs=num_cores)(
    delayed(average_chlor_a)(geom) for geom in bounds
)

gdf["chlor_avg"] = chlor_avg

Now make the H3 index for both datasets.

import h3pandas  # noqa I001

gdf["x"] = gdf["geometry"].centroid.x
gdf["y"] = gdf["geometry"].centroid.y

# The resolution goes from 0-15, using 2 b/c that is the closest to the OBIS squares.
resolution = 2

df = gdf[["x", "y", "n", "chlor_avg"]]

dfh3 = df.h3.geo_to_h3(
    resolution=resolution,
    lat_col="y",
    lng_col="x",
    set_index=True,
)


dfh3 = dfh3.drop(columns=["x", "y"]).groupby(dfh3.index.name).sum()
gdfh3 = dfh3.h3.h3_to_geo_boundary()

We need to correct the longitudes for better plotting at the dateline.

import numpy as np


def fix_geom(geom):
    from shapely import Polygon

    lons, lats = map(list, geom.boundary.coords.xy)
    if max(abs(np.diff(lons))) > 180:
        xx = [x + 360 if x < 0 else x for x in lons]
        return Polygon(zip(xx, lats))
    else:
        return geom


gdfh3["geometry"] = gdfh3["geometry"].apply(fix_geom)
gdfh3["log_n"] = np.log(gdfh3["n"])

Let’s take a quick look at the data to see if we are getting the structure we are looking for.

gdfh3

Finally, let’s plot the OBIS data at the top and PACE at the lower panel.

fig, (ax0, ax1) = plt.subplots(nrows=2, figsize=(13, 11), sharex=True, sharey=True)

ax = gdfh3.plot(
    column="n",
    ax=ax0,
    legend=True,
    norm=colors.LogNorm(),
)
gdfh3.plot(
    column="chlor_avg",
    ax=ax1,
    legend=True,
    norm=colors.LogNorm(),
)
ax0.axis([-180, 180, -90, 90])
ax0.set_title("OBIS: Number of Phytoplankton Occurrences")
ax1.set_title("PACE: Average Chlorophyll-a Concentration")
fig.tight_layout()

Done! Hopefully this example will be useful for those using these different data sources.

We would expect to see some rough alignment of higher concentrations of phytoplankton species occurrences (OBIS) with higher concentrations of Chlorophyll-a (PACE). And we do see some of that off Europe and the Gulf Stream. However, we do need to consider that the observing tool for PACE is satellite observations, so that is only presenting surface observations. Whereas OBIS occurrences are the total number of records for one geographic cell, regardless of the depth at which they observed the animal (think of a CTD profile where bottles were collected at the same station but different depths. The OBIS records would be the summation of all the occurrences for that animal at that station, regardless of depth). Additionally, PACE is a snapshot from 2025-01-01 to 2025-12-31 versus OBIS as a cummulative total from 1972-2025.