From the PO.DAAC Cookbook, to access the GitHub version of the notebook, follow this link.

Note: This notebook uses the default SWOT RiverSP collection available through Hydrocron. At the time of writing, the default collection corresponds to Version D SWOT data.

Hydrocron API: Getting Started with SWOT Time Series

Authors: Nikki Tebaldi, Cassandra Nickles, and Brandi Downs, NASA PO.DAAC

Summary:

Hydrocron is an API that repackages the SWOT_L2_HR_RiverSP_D dataset into CSV or GeoJSON formats that make time series analysis easier. This notebook will highlight how to utilize Hydrocron and convert its output into a readable dataframe from multiple SWOT reaches identified from the SWORD Database.

Requirements:

Any compute environment, local or the cloud.

Learning Objectives:

• Obtain a list of SWORD IDs for a region of interest
• Access SWOT river vector product attributes for multiple reach IDs via the Hydrocron API
• Convert accessed time series data into readable dataframe
• Plot one time series variable from multiple reaches

Cite the Hydrocron API via the following:

Frank Greguska, Nikki Tebaldi, Victoria McDonald, Vanesa Gomez Gonzalez, & Cassie Nickles. (2024). podaac/hydrocron: 1.2.0 (1.2.0). Zenodo. https://doi.org/10.5281/zenodo.11176234

Import Packages

import pandas as pd
import requests
import time
import folium
import hvplot.pandas
import holoviews as hv
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt

Obtain SWORD reach IDs

In this section, we identify a list of river reach IDs from the SWOT River Database (SWORD). One river reach ID has the format CBBBBBRRRRT and a node ID has the format, CBBBBBRRRRNNNT where C stands for continent, B for basin, R for reach, N for node, and T for type. The first 6 digits of the id (CBBBBB) are the HydroBASINS Pfafstetter level 6 basin code, where the first digit represents one of nine continental regions (1 = Africa, 2 = Europe, 3 = Siberia, 4 = Asia, 5 = Oceania, 6 = South America, 7 = North America, 8 = Arctic, 9 = Greenland), and the remaining digits are nested basin levels 2–6.

In this example, we will use river reach IDs from the Skagit River in Washington. In December 2025, major flooding occurred in western Washington driven by mountain rainfall in early December and a widespread precipitation event on December 8-11.

To identify the river reach IDs for Skagit River, go to SWORD, select the ‘North America’ tab, and select the basin that covers Washington State (basin 78). After a few moments, the individual river reaches will load on the map, as shown in the screenshot below. Zoom in to the Skagit River in northern Washington, click on the individual reaches, and note the reach IDs.

image.png

# Define the reach IDs and the Hydrocron URL

# reach IDs for Skagit River in Washington, US
reach_ids = ['78310800191','78310800181','78310800171','78310800161','78310800081','78310800071',
             '78310800061','78310800051','78310800041','78310800031','78310800021','78310800015']

HYDROCRON_URL = "https://soto.podaac.earthdatacloud.nasa.gov/hydrocron/v1/timeseries"

Query Hydrocron for time series data

Depending on how many reaches you are querying for, this section could take between seconds to minutes.

All parameter options are listed in the Hydrocron documentation: https://podaac.github.io/hydrocron/timeseries/

Information about the fields (attributes that are available in the SWOT river shapefiles) can be found here: https://podaac.github.io/hydrocron/user-guide/fields/

# define a start and end time for retrieving SWOT riverSP data
start_time = "2025-10-01T00:00:00Z"
end_time = "2026-01-31T00:00:00Z"

# define the fields you want to retrieve
fields = "reach_id,time_str,wse,wse_u,geometry"


def query_hydrocron(query_url, reach_id, start_time, end_time, fields):
    """Query Hydrocron for reach-level time series data and geometry."""

    params = {
        "feature": "Reach",
        "feature_id": reach_id,
        "output": "geojson",
        "start_time": start_time,
        "end_time": end_time,
        "fields": fields
    }

    response = requests.get(query_url, params=params)
    response.raise_for_status()
    response = response.json()

    if "results" in response:
        geojson_data = response["results"]["geojson"]
        rows = []
        for feature in geojson_data["features"]:
            rows.append(feature["properties"])
        df = pd.DataFrame(rows)
    else:
        print(f"No results for reach {reach_id}: {response}")
        geojson_data = None
        df = pd.DataFrame()

    return df, geojson_data


t0 = time.perf_counter()

results = []
geojson_results = {}

for reach in reach_ids:

    df_reach, geojson_data = query_hydrocron(
        HYDROCRON_URL,
        reach,
        start_time,
        end_time,
        fields
    )

    if not df_reach.empty:
        results.append(df_reach)

    if geojson_data is not None:
        geojson_results[reach] = geojson_data
        
if results:
    df = pd.concat(results, ignore_index=True)
    df["wse"] = pd.to_numeric(df["wse"], errors="coerce")    # errors="coerce" replaces values with NaN if they can't be converted to the indicated type rather than throwing an error
    df["wse_u"] = pd.to_numeric(df["wse_u"], errors="coerce")           
    df["time_str"] = pd.to_datetime(df["time_str"], errors="coerce")    
    df = df.dropna(subset=["time_str"])
else:
    df = pd.DataFrame()

print(f"Query time: {time.perf_counter()-t0}")

display(df)

Query time: 29.016685999988113

	reach_id	time_str	wse	wse_u	wse_units	wse_u_units
1	78310800191	2025-10-09 18:44:58+00:00	132.3752	0.09629	m	m
2	78310800191	2025-10-13 08:06:53+00:00	132.4446	0.09850	m	m
3	78310800191	2025-10-20 17:07:51+00:00	132.1659	0.17033	m	m
5	78310800191	2025-10-30 15:30:05+00:00	132.5509	0.09756	m	m
6	78310800191	2025-11-03 04:51:59+00:00	132.7871	0.10375	m	m
...	...	...	...	...	...	...
235	78310800015	2025-12-25 20:44:51+00:00	2.1087	0.09587	m	m
236	78310800015	2026-01-02 05:45:57+00:00	1.3743	0.25852	m	m
237	78310800015	2026-01-12 04:08:09+00:00	1.3693	0.09339	m	m
238	78310800015	2026-01-15 17:29:55+00:00	2.2662	0.09627	m	m
239	78310800015	2026-01-23 02:31:01+00:00	1.7769	0.21807	m	m

193 rows × 6 columns

Note: Interactive Folium maps and hvPlot visualizations may not render in GitHub’s notebook preview.
Run the notebook locally in JupyterLab to view these outputs.

# Map the river reaches

if not geojson_results:
    raise ValueError("No GeoJSON results were returned. Cannot create map.")
    
# generate colors for each reach
cmap = plt.colormaps.get_cmap("tab20").resampled(len(geojson_results))
colors = [
    mcolors.to_hex(cmap(i))
    for i in range(len(geojson_results))
]
color_dict = dict(zip(geojson_results.keys(), colors))

reach_map = folium.Map(tiles="cartodbpositron")

for reach_id, geojson_data in geojson_results.items():

    folium.GeoJson(
        geojson_data,
        name=f"Reach ID {reach_id}",
        style_function=lambda feature, reach_id=reach_id: {
            "color": color_dict[reach_id],
            "weight": 4
        },

        tooltip=folium.GeoJsonTooltip(
            fields=["reach_id"],
            aliases=["Reach ID:"]
        )

    ).add_to(reach_map)

folium.LayerControl().add_to(reach_map)

reach_map.fit_bounds(reach_map.get_bounds(), padding=(30, 30))

reach_map

Make this Notebook Trusted to load map: File -> Trust Notebook

In hvPlot, the asterisk (*) operator is used to overlay objects on the same axis. Here, we use this operator to combine multiple individual line and scatter plots into a single visual.

# Plot the time series results

# remove fill values
df_plot = df[df["wse"] > -9999].copy()

line_plots = []
scatter_plots = []

for reach_id in reach_ids:

    df_reach = df_plot[df_plot["reach_id"].astype(str) == str(reach_id)]

    color = color_dict[reach_id]

    line = df_reach.hvplot(
        x="time_str",
        y="wse",
        kind="line",
        color=color,
        label=f"Reach ID {reach_id}"
    ).opts(tools=[])  # tools=[] removes tooltips from the line segments

    scatter = df_reach.hvplot(
        x="time_str",
        y="wse",
        kind="scatter",
        color=color,
        hover_cols=["reach_id"]
    )

    line_plots.append(line)
    scatter_plots.append(scatter)

# build up one large overlay plot from individual line and scatter plots
if not line_plots:
    raise ValueError("No valid data available to plot.")
plot = line_plots[0]
for p in line_plots[1:]:
    plot *= p
for p in scatter_plots:
    plot *= p

plot.opts(
    width=900,
    height=500,
    xlabel="Date",
    ylabel="Water Surface Elevation (m)",
    title="WSE by Reach ID and Observation Date",
    xrotation=90,
    legend_position="right",
    show_grid=True
)

The plotted results show an increase in water surface elevation on December 11, 2025 for the reach IDs with measurements on that date, coinciding with the time period of flooding.

However, for reach IDs 78310800031 and 78310800021, there appear to be some anomalous data points, in particular for 2025-12-04 and 2026-01-02. We can investigate these data points further using the water surface elevation uncertainty (wse_u) parameter. wse_u is the total uncertainty (random and systematic) in the reach WSE and has units of meters. High wse_u values are not automatically invalid, but they can help identify measurements that may be less reliable and worth reviewing before interpretation.

# Plot the water surface elevation uncertainty

# remove fill values
df_plot = df[
    (df["wse"] > -9999) &
    (df["wse_u"] > -9999)
].copy()

scatter_plots = []

for reach_id in reach_ids:

    df_reach = df_plot[df_plot["reach_id"].astype(str) == str(reach_id)]

    scatter = df_reach.hvplot(
        x="time_str",
        y="wse_u",
        kind="scatter",
        color=color_dict[reach_id],
        label=f"Reach ID {reach_id}",
        hover_cols=["reach_id"]
    )

    scatter_plots.append(scatter)

# build up one large overlay plot from individual scatter plots
if not scatter_plots:
    raise ValueError("No valid uncertainty data available to plot.")
plot = scatter_plots[0]
for p in scatter_plots[1:]:
    plot *= p
    
# add a horizontal line corresponding to an uncertainty level of 0.5 m
threshold = hv.HLine(0.5).opts(line_dash="dashed", color="gray")

# place annotation above horizontal line
x_pos = df_plot["time_str"].min()
y_pos = 0.52
label = hv.Text(x_pos, y_pos, "0.5 m").opts(text_align="left", text_color="gray", text_font_size="10pt")

# add horizontal line and annotation to the plot
plot = plot * threshold * label

plot.opts(
    width=900,
    height=500,
    xlabel="Date",
    ylabel="WSE Uncertainty (m)",
    title="WSE Uncertainty per Measurement",
    xrotation=90,
    legend_position="right",
    yformatter="%.2f",
    yticks=5,
    ylim=(0,1.0),
    show_grid=True
)

# Plot the uncertainty-filtered time series results


# set water surface elevation uncertainty threshold
WSE_U_THRESHOLD = 0.5


# remove fill values and filter out data points with high uncertainty
df_plot = df[
    (df["wse"] > -9999) &
    (df["wse_u"] > -9999) &
    (df["wse_u"] <= WSE_U_THRESHOLD)
].copy()


line_plots = []
scatter_plots = []

for reach_id in reach_ids:

    df_reach = df_plot[df_plot["reach_id"].astype(str) == str(reach_id)]

    color = color_dict[reach_id]

    line = df_reach.hvplot(
        x="time_str",
        y="wse",
        kind="line",
        color=color,
        label=f"Reach ID {reach_id}"
    ).opts(tools=[])

    scatter = df_reach.hvplot(
        x="time_str",
        y="wse",
        kind="scatter",
        color=color,
        hover_cols=["wse_u","reach_id"]
    )

    line_plots.append(line)
    scatter_plots.append(scatter)

# build up one large overlay plot from individual line and scatter plots
if not line_plots:
    raise ValueError("No valid data available to plot.")
plot = line_plots[0]
for p in line_plots[1:]:
    plot *= p
for p in scatter_plots:
    plot *= p

plot.opts(
    width=900,
    height=500,
    xlabel="Date",
    ylabel="Water Surface Elevation (m)",
    title="WSE After Removing High Uncertainty Data Points",
    xrotation=90,
    legend_position="right",
    show_grid=True
)