How fast the ocean surface moves: 197 million hourly positions of 19,396 drifting buoys

[jvdd]

The NOAA Global Drifter Program has tracked 19,396 drifting buoys since 1987. Its hourly dataset holds 197 million positions and velocities, from October 1987 to October 2022. The map shows the mean speed of the buoys per 0.2° cell.
-> Brush the speed histogram and the map keeps only the buoy-hours in that speed range. Zoom the map and it re-bins to the new view.

• Data: NOAA Global Drifter Program, hourly dataset v2.01.1, public on AWS: https://registry.opendata.aws/noaa-oar-hourly-gdp/ (file s3://noaa-oar-hourly-gdp-pds/latest/gdp-v2.01.1.parquet, doi:10.25921/x46c-3620)
• Install: pip install "flexviz>=0.1.0b5" pyarrow (and shapely for insight.py)
• Size: one 8.2 GB Parquet file, with one row per buoy and the observations as nested lists. etl.py flattens it to 197 million rows (2.1 GB of Parquet) in about 2 minutes, with a peak of 14 GB of RAM.

Steps

1. Download the file to data/: curl -o data/gdp-v2.01.1.parquet https://noaa-oar-hourly-gdp-pds.s3.amazonaws.com/latest/gdp-v2.01.1.parquet
2. python etl.py writes data/drifters.parquet, with one row per buoy-hour.
3. python dashboard.py opens the live dashboard. It keeps the 197 million rows in memory (about 3.2 GB), so a brush answers in about a second.

Speed

The hourly velocities include tides and inertial oscillations (circular motions that storms start). These can be fast where the mean current is slow. So the speed here is the mean velocity of each buoy over one UTC day. With the raw hourly speed, the map of fast water is noisy.

Buoys that lost their drogue (a sea anchor at 15 m depth) move a little with the wind. They are 60 % of the buoy-hours, but only 53 % of the hours faster than 1 m/s. The map with drogued buoys only shows the same currents, so all buoys are included.

What the data shows

• The median daily speed is 0.20 m/s (0.7 km/h). Only 0.6 % of the buoy-hours are faster than 1 m/s. Half of those are in 1.6 % of the ocean area (1° cells with data, weighted by latitude).
• The fast currents run along the east coasts of the continents: the Gulf Stream, the Kuroshio, the Agulhas Current and the East Australian Current. Within 2° of an east coast, 5.7 % of the buoy-hours are faster than 1 m/s. Within 2° of a west coast, it is 0.2 %: 27 times less (60° S to 60° N, see insight.py).
• Near the equator (5° S to 5° N) the mean speed is 0.40 m/s. That is about twice the 0.21 m/s of the subtropics (20° to 30° N and S).

The video

The video is a headless-browser capture of this dashboard at 1800 × 790 bins (0.2°). The brush moves from 0 to 1 m/s in steps of 0.013 m/s. During the sweep, the color range starts at the brush value: at "speed > X", the colors run from X to max(1, X + 0.4) m/s. So the strong currents stay bright and single fast hours stay dark. The live dashboard keeps the fixed range of 0 to 1 m/s. The title, the filter text and the current names are overlays.

Code

etl.py: flatten the buoy file to one row per buoy-hour

"""Flatten the NOAA Global Drifter Program hourly file to one row per buoy-hour.

Source: s3://noaa-oar-hourly-gdp-pds/latest/gdp-v2.01.1.parquet (one row per
buoy, observations as nested lists). Output: data/drifters.parquet with lon,
lat, time (days since 1970), buoy id, the hourly speed, the daily-mean speed and velocity, and the drogue flag.

Hourly velocities include tides and inertial oscillations (circular motions
that storms start), which can be fast where the mean current is slow. The mean
of ve and vn over one UTC day per buoy removes most of that, so speed_d shows
the current that carries the buoy.

Usage: python etl.py
"""

from pathlib import Path

import polars as pl
import pyarrow.parquet as pq

DATA = Path(__file__).parent / "data"
COLS = ["id", "time", "lon", "lat", "ve", "vn", "drogue_status"]

table = pq.read_table(DATA / "gdp-v2.01.1.parquet", columns=[f"obs.{c}" for c in COLS])
df = (
    pl.from_arrow(table)
    .explode(COLS[1:])
    .rename({"drogue_status": "drogued"})
    .with_columns(day=(pl.col("time") // 86400).cast(pl.Int32))
    .with_columns(
        ve_d=pl.col("ve").mean().over("id", "day"),
        vn_d=pl.col("vn").mean().over("id", "day"),
    )
    .select(
        pl.col("lon", "lat").cast(pl.Float32),
        speed_h=(pl.col("ve") ** 2 + pl.col("vn") ** 2).sqrt().cast(pl.Float32),
        speed_d=(pl.col("ve_d") ** 2 + pl.col("vn_d") ** 2).sqrt().cast(pl.Float32),
        ve_d=pl.col("ve_d").cast(pl.Float32),
        vn_d=pl.col("vn_d").cast(pl.Float32),
        year=pl.from_epoch(pl.col("time").cast(pl.Int64)).dt.year().cast(pl.Int16),
        # Days since 1970; Float32 keeps about 3 minutes of resolution.
        t=(pl.col("time") / 86400).cast(pl.Float32),
        id=pl.col("id"),
        drogued=pl.col("drogued"),
    )
)
df.write_parquet(DATA / "drifters.parquet")
print(df.height, "rows")
print(df.describe())

dashboard.py: the live dashboard

"""Live FlexViz dashboard: how fast the ocean surface moves.

Brush the speed histogram and the map keeps only the buoy-hours in that speed
range. Zoom the map and it re-bins to the new view.

Run: python dashboard.py   (data/drifters.parquet comes from etl.py)
"""

import polars as pl
from flexviz import Dashboard, GridItem, LayoutSpec

# 197 million buoy-hours. collect() keeps them in RAM (about 3.2 GB), so a
# brush answers in about a second. Pass the LazyFrame instead to stream them
# from disk.
drifters = (
    pl.scan_parquet("data/drifters.parquet")
    .select(
        "lon",
        "lat",
        # The histogram ends at 2 m/s ("2+"): 0.01 % of buoy-hours are faster.
        speed=pl.col("speed_d").clip(upper_bound=2.0),
        speed_ms=pl.col("speed_d"),
    )
    .collect()
)

dash = Dashboard(drifters, cache=True)

speed_map = dash.add_figure(title="Mean speed of drifting buoys, 1987–2022 (m/s)")
speed_map.add_histogram2d(
    x="lon",
    y="lat",
    x_bins=360,
    y_bins=160,
    # x_bins=1800,  # 0.2 degree per bin, as in the video
    # y_bins=790,
    z="speed_ms",
    histfunc="mean",
    color_scale="YlGnBu",
    color_range=(0.0, 1.0),
)
speed_map.update_layout(
    plot_bgcolor="black",
    xaxis={"visible": False},
    yaxis={"visible": False, "scaleanchor": "x"},
)

speeds = dash.add_figure(title="Buoy-hours by speed")
speeds.add_histogram(x="speed", bins=80).xlabel("daily-mean speed (m/s)")

map_uid, speeds_uid = (f.uid for f in dash.to_spec().figures)
dash.show(
    port=8050,
    source_name="drifters",
    # A live-brush cube for a fine map is too big for the browser, so brushes
    # use the normal request on mouse-up.
    live_brush="off",
    layout=LayoutSpec(
        grid_items=[
            GridItem(fig_uid=map_uid, x=0, y=0, w=12, h=8),
            GridItem(fig_uid=speeds_uid, x=0, y=8, w=12, h=3),
        ]
    ),
)

insight.py: the east-coast, west-coast and latitude numbers

"""Check the title claims: fast water along east coasts of continents and near the equator.

Needs the Natural Earth 1:50m countries in data/:
curl -o data/ne_50m_admin_0_countries.geojson https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_50m_admin_0_countries.geojson
"""
import json
import numpy as np
import polars as pl
import shapely
from shapely.geometry import shape

DEG, K = 0.25, 8  # 0.25 degree cells; coast strip = land within K cells (2 degrees) to the west/east
W, H = int(360 / DEG), int(120 / DEG)  # latitudes -60..60 (no sea ice / Antarctic shelf)
cx = -180 + (np.arange(W) + 0.5) * DEG
cy = -60 + (np.arange(H) + 0.5) * DEG
gx, gy = np.meshgrid(cx, cy)
land = np.zeros((H, W), bool)
for f in json.load(open("data/ne_50m_admin_0_countries.geojson"))["features"]:
    g = shape(f["geometry"])
    for poly in getattr(g, "geoms", [g]):
        if poly.area < 5:  # square degrees: drop small islands, keep Japan, Cuba, Madagascar, NZ, UK
            continue
        x0, y0, x1, y1 = poly.bounds
        sel = (gx >= x0) & (gx <= x1) & (gy >= y0) & (gy <= y1)
        land[sel] |= shapely.contains_xy(poly, gx[sel], gy[sel])
west = np.zeros_like(land); east = np.zeros_like(land)
for k in range(1, K + 1):
    west |= np.roll(land, k, axis=1)   # land k cells to the west
    east |= np.roll(land, -k, axis=1)  # land k cells to the east
zone = np.full((H, W), "open ocean", dtype=object)
zone[west & ~east] = "east-coast strip"
zone[east & ~west] = "west-coast strip"
zone[west & east] = "narrow sea"
zone[land] = "land"
zdf = pl.DataFrame({"ix": np.tile(np.arange(W), H), "iy": np.repeat(np.arange(H), W), "zone": zone.ravel().astype(str)}).with_columns(pl.col("ix", "iy").cast(pl.Int32))

lf = pl.scan_parquet("data/drifters.parquet").filter(pl.col("lat").is_between(-60, 60, closed="left"))
cells = (lf.select(ix=((pl.col("lon") + 180) / DEG).floor().clip(0, W - 1).cast(pl.Int32),
                   iy=((pl.col("lat") + 60) / DEG).floor().clip(0, H - 1).cast(pl.Int32), s="speed_d")
           .group_by("ix", "iy").agg(n=pl.len(), ssum=pl.col("s").sum(), fast=(pl.col("s") > 1).sum(),
                                     f05=(pl.col("s") > 0.5).sum()).collect()
           .join(zdf, on=["ix", "iy"]))
tot_fast = cells["fast"].sum()
z = (cells.group_by("zone").agg(hours=pl.col("n").sum(), mean=pl.col("ssum").sum() / pl.col("n").sum(),
                                gt05=pl.col("f05").sum() / pl.col("n").sum(), gt1=pl.col("fast").sum() / pl.col("n").sum(),
                                share_fast=pl.col("fast").sum() / tot_fast)
     .sort("mean", descending=True))
with pl.Config(tbl_rows=10, float_precision=4, tbl_width_chars=140):
    print(z)
bands = (lf.with_columns(band=(pl.col("lat") / 5).floor() * 5).group_by("band")
           .agg(mean=pl.col("speed_d").mean(), gt1=(pl.col("speed_d") > 1).mean(), hours=pl.len()).sort("band").collect())
with pl.Config(tbl_rows=30, float_precision=3):
    print(bands)