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

SWOT Hydrology Dataset Exploration on a local machine

Accessing and Visualizing SWOT Datasets

Requirement:

Local compute environment e.g. laptop, server: this tutorial can be run on your local machine.

Learning Objectives:

Access SWOT HR data products (archived in NASA Earthdata Cloud) by downloading to local machine
Visualize accessed data for a quick check

SWOT Level 2 KaRIn High Rate Version D Datasets:

River Vector Shapefile - SWOT_L2_HR_RIVERSP_D
Lake Vector Shapefile - SWOT_L2_HR_LAKESP_D
Water Mask Pixel Cloud NetCDF - SWOT_L2_HR_PIXC_D
Water Mask Pixel Cloud Vector Attribute NetCDF - SWOT_L2_HR_PIXCVec_D
Raster NetCDF - SWOT_L2_HR_Raster_D
Single Look Complex Data product - SWOT_L1B_HR_SLC_D

Notebook Author: Cassie Nickles, NASA PO.DAAC (Feb 2024) || Other Contributors: Zoe Walschots (PO.DAAC Summer Intern 2023), Catalina Taglialatela (NASA PO.DAAC), Luis Lopez (NASA NSIDC DAAC), Brent Williams (NASA JPL), Nicholas Tarpinian (NASA PO.DAAC) (Nov 2026)

Last updated: November 2026

Libraries Needed

import os
import glob
import h5netcdf
import xarray as xr
import pandas as pd
import geopandas as gpd
import contextily as cx
import numpy as np
import matplotlib.pyplot as plt
import hvplot.xarray
import zipfile
import earthaccess
import warnings

warnings.filterwarnings("ignore", category=FutureWarning, module="earthaccess")

Earthdata Login

An Earthdata Login account is required to access data, as well as discover restricted data, from the NASA Earthdata system. Thus, to access NASA data, you need Earthdata Login. If you don’t already have one, please visit https://urs.earthdata.nasa.gov to register and manage your Earthdata Login account. This account is free to create and only takes a moment to set up. We use earthaccess to authenticate your login credentials below.

auth = earthaccess.login()

Single File Access

1. River Vector Shapefiles

The https access link can be found using earthaccess data search. Since this collection consists of Reach and Node files, we need to extract only the granule for the Reach file. We do this by filtering for the ‘Reach’ title in the data link.

Alternatively, Earthdata Search (see tutorial) can be used to manually search in a GUI interface.

For additional tips on spatial searching of SWOT HR L2 data, see also PO.DAAC Cookbook - SWOT Chapter tips section.

Search for the data of interest

#Retrieves granule from the day we want, in this case by passing to `earthdata.search_data` function the data collection shortname, temporal bounds, and filter by wildcards
river_results = earthaccess.search_data(short_name = 'SWOT_L2_HR_RiverSP_D', 
                                        temporal = ('2024-02-01 00:00:00', '2024-07-15 23:59:59'), # can also specify by time
                                        granule_name = '*Reach*_287_NA*') # here we filter by Reach files (not node), pass=287, continent code=NA
print(f"Found {len(river_results)} granules")

Dowload, unzip, read the data

Let’s download the first data file! earthaccess.download has a list as the input format, so we need to put brackets around the single file we pass.

earthaccess.download([river_results[0]], "./data_downloads")

The native format for this data is a .zip file, and we want the .shp file within the .zip file, so we must first extract the data to open it. First, we’ll programmatically get the filename we just downloaded, and then extract all data to the data_downloads folder.

filename = earthaccess.results.DataGranule.data_links(river_results[0], access='external')
filename = filename[0].split("/")[-1]
filename

with zipfile.ZipFile(f'data_downloads/{filename}', 'r') as zip_ref:
    zip_ref.extractall('data_downloads')

Open the shapefile using geopandas

filename_shp = filename.replace('.zip','.shp')

SWOT_HR_shp1 = gpd.read_file(f'data_downloads/{filename_shp}') 

#view the attribute table
SWOT_HR_shp1

Quickly plot the SWOT river data

# Simple plot
fig, ax = plt.subplots(figsize=(10,7))
SWOT_HR_shp1.plot(ax=ax, color='black')
cx.add_basemap(ax, crs=SWOT_HR_shp1.crs, source=cx.providers.OpenTopoMap)

# Another way to plot geopandas dataframes is with `explore`, which also plots a basemap
#SWOT_HR_shp1.explore()

2. Lake Vector Shapefiles

The lake vector shapefiles can be accessed in the same way as the river shapefiles above.

For additional tips on spatial searching of SWOT HR L2 data, see also PO.DAAC Cookbook - SWOT Chapter tips section.

Search for data of interest

lake_results = earthaccess.search_data(short_name = 'SWOT_L2_HR_LakeSP_2.0', 
                                        temporal = ('2024-02-01 00:00:00', '2024-07-15 23:59:59'), # can also specify by time
                                        granule_name = '*Prior*_287_NA*') # here we filter by files with 'Prior' in the name (This collection has three options: Obs, Unassigned, and Prior), pass 287 and continent code=NA
print(f"Found {len(lake_results)} granules")

Let’s download the first data file! earthaccess.download has a list as the input format, so we need to put brackets around the single file we pass.

earthaccess.download([lake_results[0]], "./data_downloads")

filename2 = earthaccess.results.DataGranule.data_links(lake_results[0], access='external')
filename2 = filename2[0].split("/")[-1]
filename2

with zipfile.ZipFile(f'data_downloads/{filename2}', 'r') as zip_ref:
    zip_ref.extractall('data_downloads')

Open the shapefile using geopandas

filename_shp2 = filename2.replace('.zip','.shp')
filename_shp2

SWOT_HR_shp2 = gpd.read_file(f'data_downloads/{filename_shp2}') 

#view the attribute table
SWOT_HR_shp2

Quickly plot the SWOT lakes data

fig, ax = plt.subplots(figsize=(10,7))
SWOT_HR_shp2.plot(ax=ax, color='black')
cx.add_basemap(ax, crs=SWOT_HR_shp2.crs, source=cx.providers.OpenTopoMap)

Accessing the remaining files is different than the shp files above. We do not need to extract the shapefiles from a zip file because the following SWOT HR collections are stored in netCDF files in the cloud. For the rest of the products, we will open via xarray, not geopandas.

3. Water Mask Pixel Cloud NetCDF

Search for data collection and time of interest

pixc_results = earthaccess.search_data(short_name = 'SWOT_L2_HR_PIXC_D', 
                                        temporal = ('2024-01-25 00:00:00', '2024-05-29 23:59:59'), # can also specify by time
                                        bounding_box = (-106.62, 38.809, -106.54, 38.859)) # Lake Travis near Austin, TX
print(f"Found {len(pixc_results)} granules")

Let’s download one data file! earthaccess.download has a list as the input format, so we need to put brackets around the single file we pass.

earthaccess.download([pixc_results[0]], "./data_downloads")

Open data using xarray

The pixel cloud netCDF files are formatted with three groups titled, “pixel cloud”, “tvp”, or “noise” (more detail here). In order to access the coordinates and variables within the file, a group must be specified when calling xarray open_dataset.

ds_PIXC = xr.open_mfdataset("data_downloads/SWOT_L2_HR_PIXC_*.nc", group = 'pixel_cloud', engine='h5netcdf')
ds_PIXC

For plotting PIXC using classification and quality flags

# mask to get good water pixels
mask = np.where(np.logical_and(ds_PIXC.classification > 2, ds_PIXC.geolocation_qual <28800)) 
#For more conservative filtering, use ds_PIXC.geolocation_qual <4 to also remove suspect values 

plt.scatter(x=ds_PIXC.longitude[mask], y=ds_PIXC.latitude[mask], c=ds_PIXC.height[mask])
plt.colorbar().set_label('Height (m)')

4. Water Mask Pixel Cloud Vector Attribute NetCDF

Search for data of interest

#Let's plot the same pass and tile as the above
pixcvec_results = earthaccess.search_data(short_name = 'SWOT_L2_HR_PIXCVec_D', 
                                          granule_name = '*010_106_086R*') #The same cycle, pass and tile as previously downloaded
print(f"Found {len(pixcvec_results)} granules")

Let’s download the first data file! earthaccess.download has a list as the input format, so we need to put brackets around the single file we pass.

earthaccess.download([pixcvec_results[0]], "./data_downloads")

Open data using xarray

First, we’ll programmatically get the filename we just downloaded and then view the file via xarray.

ds_PIXCVEC = xr.open_mfdataset("data_downloads/SWOT_L2_HR_PIXCVec_*.nc", decode_cf=False,  engine='h5netcdf')
ds_PIXCVEC

Simple plot

pixcvec_htvals = ds_PIXCVEC.height_vectorproc.compute()
pixcvec_latvals = ds_PIXCVEC.latitude_vectorproc.compute()
pixcvec_lonvals = ds_PIXCVEC.longitude_vectorproc.compute()

#Before plotting, we set all fill values to nan so that the graph shows up better spatially
pixcvec_htvals[pixcvec_htvals > 15000] = np.nan
pixcvec_latvals[pixcvec_latvals == 0] = np.nan
pixcvec_lonvals[pixcvec_lonvals == 0] = np.nan

plt.scatter(x=pixcvec_lonvals, y=pixcvec_latvals, c=pixcvec_htvals)
plt.colorbar().set_label('Height (m)')

5. Raster NetCDF

Search for data of interest

raster_results = earthaccess.search_data(short_name = 'SWOT_L2_HR_Raster_D', 
                                        temporal = ('2024-01-25 00:00:00', '2024-07-15 23:59:59'), # can also specify by time
                                        granule_name = '*100m*', # here we filter by files with '100m' in the name (This collection has two resolution options: 100m & 250m)
                                        bounding_box = (-106.62, 38.809, -106.54, 38.859)) # Lake Travis near Austin, TX
print(f"Found {len(raster_results)} granules")

Let’s download one data file.

earthaccess.download([raster_results[0]], "./data_downloads")

Open data with xarray

First, we’ll programmatically get the filename we just downloaded and then view the file via xarray.

ds_raster = xr.open_mfdataset(f'data_downloads/SWOT_L2_HR_Raster*', engine='h5netcdf')
ds_raster

Quick interactive plot with `hvplot`

Note: this is not filtered by quality

ds_raster.wse.hvplot.quadmesh(y='y', x='x')

6. SLC

Search for data collection and time of interest

slc_results = earthaccess.search_data(short_name = 'SWOT_L1B_HR_SLC_D', 
                                        #temporal = ('2024-01-01 00:00:00', '2024-01-30 23:59:59'), # can also specify by time
                                        bounding_box = (-106.62, 38.809, -106.54, 38.859)) # Lake Travis near Austin, TX
print(f"Found {len(slc_results)} granules")

Let’s download one data file! earthaccess.download has a list as the input format, so we need to put brackets around the single file we pass.

earthaccess.download([slc_results[0]], "./data_downloads")

Open data using xarray

The L1B_HR_SLC product file contains five NetCDF data group called the slc, xfactor, noise, tvp, and grdem groups. More info can be found in the product description document within the dataset table for each group.

# Find all SLC granules and open them without loading the data into memory.
# chunks={} makes xarray read the data in small pieces later, when a result is requested.
paths = sorted(glob.glob("data_downloads/SWOT_L1B_HR_SLC*.nc"))
dsets = [xr.open_dataset(p, group='slc', engine='h5netcdf', chunks={}) for p in paths]

# Why not open_mfdataset? The files have different sizes (e.g. 22956 vs 22804 lines),
# and SLC data has no shared coordinates to align on, so we stack them manually.

# The range dimension (num_pixels) differs by one pixel between files (5622 vs 5623).
# Trim every file to the smallest width so they can be stacked.
# This drops one range pixel from the wider files. That is fine for viewing, but not
# for phase-sensitive processing that relies on the original pixel grid.
n = min(d.sizes['num_pixels'] for d in dsets)
trimmed = [d.isel(num_pixels=slice(0, n)) for d in dsets]

# Stack the granules end to end along-track. The 'minimal'/'override' options keep
# variables without a num_lines dimension (scalars, metadata) from being duplicated,
# using the first file's copy.
ds_SLC = xr.concat(trimmed, dim='num_lines',
                   data_vars='minimal', coords='minimal', compat='override')
ds_SLC

Plot SLC data

# function to get interferogram from slc images

# Multilook size: number of pixels averaged into one output pixel (azimuth x range).
# Larger values use less memory but give a coarser image.
looks_az, looks_rg = 30, 5

def interfere_multilook(slc, looks_az, looks_rg):
    # Combine the real and imaginary layers into complex arrays.
    # These lines only define the calculation; no data is read yet.
    slc_plus_y = slc.slc_plus_y[:, :, 0] + 1j * slc.slc_plus_y[:, :, 1]
    slc_minus_y = slc.slc_minus_y[:, :, 0] + 1j * slc.slc_minus_y[:, :, 1]

    # Interferogram: plus_y times the complex conjugate of minus_y
    ifgram = slc_plus_y * np.conj(slc_minus_y)

    # Average over blocks BEFORE plotting so only a small array is held in memory.
    # Average the complex values, not the phase angles, so the phase wrap at +/- pi
    # doesn't corrupt the result. boundary='trim' drops leftover edge pixels.
    ifgram_ml = ifgram.coarsen(num_lines=looks_az, num_pixels=looks_rg, boundary='trim').mean()

    # .compute() runs the calculation now: the files are read in pieces and only
    # the small averaged result is kept in memory.
    ifgram_ml = ifgram_ml.compute()

    mag = np.abs(ifgram_ml)
    phase = np.arctan2(np.imag(ifgram_ml), np.real(ifgram_ml))
    return ifgram_ml, mag, phase

ifgram, mag, phase = interfere_multilook(ds_SLC, looks_az, looks_rg)
print(ifgram.shape)

# function to get interferogram from slc images

# ignore the divide-by-zero warning in log10
warnings.filterwarnings("ignore")

# make a plot of slc interferogram
# Magnitude in dB (grayscale) with the phase overlaid semi-transparent
plt.figure(figsize=(8, 10))
plt.imshow(10*np.log10(mag), interpolation='none', cmap='gray', aspect='auto')
plt.clim((0, 90))
plt.colorbar().set_label('interferogram magnitude (dB)')
plt.imshow(phase, interpolation='none', cmap='hsv', alpha=0.5, aspect='auto')
plt.colorbar().set_label('interferometric phase')

# Amplitude image for the plus_y channel, using the same multilook size as the interferogram

# Combine the real and imaginary layers into a complex array.
# This line only defines the calculation; no data is read yet.
slc_plus_y = ds_SLC.slc_plus_y[:, :, 0] + 1j * ds_SLC.slc_plus_y[:, :, 1]

# Average over blocks (looks_az x looks_rg pixels) first, then run the calculation with
# .compute(). The files are read in pieces and only the small averaged result stays in memory.
power_plus_y = (np.abs(slc_plus_y)
                .coarsen(num_lines=looks_az, num_pixels=looks_rg, boundary='trim')
                .mean()
                .compute())

plt.figure(figsize=(8, 10))
plt.imshow(10*np.log10(power_plus_y), interpolation='none', cmap='gray', aspect='auto')
plt.clim((0, 45))
plt.colorbar().set_label('interferometric power (plus_y channel)')

SWOT Hydrology Dataset Exploration on a local machine

Accessing and Visualizing SWOT Datasets

Requirement:

Learning Objectives:

SWOT Level 2 KaRIn High Rate Version D Datasets:

Libraries Needed

Earthdata Login

Single File Access

1. River Vector Shapefiles

Search for the data of interest

Dowload, unzip, read the data

Quickly plot the SWOT river data

2. Lake Vector Shapefiles

Search for data of interest

Quickly plot the SWOT lakes data

3. Water Mask Pixel Cloud NetCDF

Search for data collection and time of interest

Open data using xarray

For plotting PIXC using classification and quality flags

4. Water Mask Pixel Cloud Vector Attribute NetCDF

Search for data of interest

Open data using xarray

Simple plot

5. Raster NetCDF

Search for data of interest

Open data with xarray

Quick interactive plot with hvplot

6. SLC

Search for data collection and time of interest

Open data using xarray

Plot SLC data

Quick interactive plot with `hvplot`