For an updated notebook using the latest data, see this notebook in the PO.DAAC Cookbook

SWOT Hydrology Science Application Tutorial on the Cloud

Retrieving SWOT attributes (WSE, width, slope) and plotting a longitudinal profile along a river or over a basin

Requirement

This tutorial can only be run in an AWS cloud instance running in us-west-2: NASA Earthdata Cloud data in S3 can be directly accessed via earthaccess python library; this access is limited to requests made within the US West (Oregon) (code: us-west-2) AWS region.

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. 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.

This code runs using SWOT Level 2 Data Products (Version 1.1).

Notebook Authors: Arnaud Cerbelaud, Jeffrey Wade, NASA Jet Propulsion Laboratory - California Institute of Technology (Jan 2024)

Learning Objectives

  1. Retrieve SWOT hydrological attributes on river reaches within the AWS cloud (Cal/Val data). Query reaches by:
    • River name
    • Spatial bounding box
    • Downstream tracing from reach id (e.g. headwater to outlet) for river longitudinal profiles
    • Upstream tracing from reach id (e.g. outlet to full river network) for watershed analysis
  2. Plot a time series of WSE, width, slope data on the filtered data
  3. Visualize an interactive map of WSE, width, slope data on the filtered data

Last updated: 2 Feb 2024

Import Packages

import glob
import os
import requests
import s3fs
import fiona
import netCDF4 as nc
import h5netcdf
import xarray as xr
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import hvplot.xarray
import earthaccess
from earthaccess import Auth, DataCollections, DataGranules, Store
from pathlib import Path
import zipfile

Authenticate

Authenticate your Earthdata Login (EDL) information using the earthaccess python package as follows:

earthaccess.login() # Login with your EDL credentials if asked
Enter your Earthdata Login username:  cerbelaud
Enter your Earthdata password:  ········
<earthaccess.auth.Auth at 0x7f114ba7bd30>

1. Retrieve SWOT hydrological attributes on river reaches within the AWS cloud (Cal/Val data)

What data should we download?

  • Optional step: Get the .kmz file of SWOT passes/swaths (Version 1.1) and import it into Google Earth for visualization
  • Determine which pass number corresponds to the river/basin you want to look at! #### Search for multiple days of data
# Enter pass number
pass_number    = "009"   # e.g. 009 for Connecticut in NA, 003 for Rhine in EU
# Enter continent code
continent_code = "NA"     # e.g. "AF", "NA", "EU", "SI", "AS", "AU", "SA", "AR", "GR"

# Retrieves granule from the day we want, in this case by passing to `earthdata.search_data` function the data collection shortname and temporal bounds
river_results = earthaccess.search_data(short_name = 'SWOT_L2_HR_RIVERSP_1.1', 
                                        temporal = ('2023-04-08 00:00:00', '2023-04-12 23:59:59'),
                                        granule_name = "*Reach*_" + pass_number + "_" + continent_code + "*")

# Create fiona session to read data
fs_s3 = earthaccess.get_s3fs_session(results=river_results)
fiona_session=fiona.session.AWSSession(
        aws_access_key_id=fs_s3.storage_options["key"],
        aws_secret_access_key=fs_s3.storage_options["secret"],
        aws_session_token=fs_s3.storage_options["token"]
    )
Granules found: 5

Unzip selected files in Fiona session

# Initialize list of shapefiles containing all dates
SWOT_HR_shps = []

# Loop through queried granules to stack all acquisition dates
for j in range(len(river_results)):
    
    # Get the link for each zip file
    river_link = earthaccess.results.DataGranule.data_links(river_results[j], access='direct')[0]
    
    # We use the zip+ prefix so fiona knows that we are operating on a zip file
    river_shp_url = f"zip+{river_link}"
    
    # Read shapefile
    with fiona.Env(session=fiona_session):
        SWOT_HR_shps.append(gpd.read_file(river_shp_url)) 

Aggregate unzipped files into dataframe

# Combine granules from all acquisition dates into one dataframe
SWOT_HR_df = gpd.GeoDataFrame(pd.concat(SWOT_HR_shps, ignore_index=True))

# Sort dataframe by reach_id and time
SWOT_HR_df = SWOT_HR_df.sort_values(['reach_id', 'time'])

SWOT_HR_df
reach_id time time_tai time_str p_lat p_lon river_name wse wse_u wse_r_u ... p_wid_var p_n_nodes p_dist_out p_length p_maf p_dam_id p_n_ch_max p_n_ch_mod p_low_slp geometry
1958 72270400231 -1.000000e+12 -1.000000e+12 no_data 56.794630 -66.461122 no_data -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 1618.079 53 230079.472 10630.774958 -1.000000e+12 0 1 1 0 LINESTRING (-66.48682 56.83442, -66.48645 56.8...
2641 72270400231 -1.000000e+12 -1.000000e+12 no_data 56.794630 -66.461122 no_data -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 1618.079 53 230079.472 10630.774958 -1.000000e+12 0 1 1 0 LINESTRING (-66.48682 56.83442, -66.48645 56.8...
1959 72270400241 -1.000000e+12 -1.000000e+12 no_data 56.744361 -66.341905 no_data -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 4831.128 53 240728.430 10648.958289 -1.000000e+12 0 2 1 0 LINESTRING (-66.41602 56.75969, -66.41553 56.7...
2642 72270400241 -1.000000e+12 -1.000000e+12 no_data 56.744361 -66.341905 no_data -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 4831.128 53 240728.430 10648.958289 -1.000000e+12 0 2 1 0 LINESTRING (-66.41602 56.75969, -66.41553 56.7...
1960 72270400251 7.344991e+08 7.344991e+08 2023-04-11T03:31:02Z 56.676087 -66.219811 no_data 2.855340e+02 1.008930e+00 1.004910e+00 ... 412.975 58 252286.769 11558.339052 -1.000000e+12 0 2 1 0 LINESTRING (-66.27812 56.71437, -66.27775 56.7...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
615 73130000061 -1.000000e+12 -1.000000e+12 no_data 41.424945 -73.230960 Housatonic River -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 10469.706 82 47161.191 16303.793847 -1.000000e+12 0 2 1 0 LINESTRING (-73.17436 41.38295, -73.17465 41.3...
1283 73130000061 -1.000000e+12 -1.000000e+12 no_data 41.424945 -73.230960 Housatonic River -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 10469.706 82 47161.191 16303.793847 -1.000000e+12 0 2 1 0 LINESTRING (-73.17436 41.38295, -73.17465 41.3...
1957 73130000061 -1.000000e+12 -1.000000e+12 no_data 41.424945 -73.230960 Housatonic River -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 10469.706 82 47161.191 16303.793847 -1.000000e+12 0 2 1 0 LINESTRING (-73.17436 41.38295, -73.17465 41.3...
2640 73130000061 -1.000000e+12 -1.000000e+12 no_data 41.424945 -73.230960 Housatonic River -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 10469.706 82 47161.191 16303.793847 -1.000000e+12 0 2 1 0 LINESTRING (-73.17436 41.38295, -73.17465 41.3...
3312 73130000061 -1.000000e+12 -1.000000e+12 no_data 41.424945 -73.230960 Housatonic River -1.000000e+12 -1.000000e+12 -1.000000e+12 ... 10469.706 82 47161.191 16303.793847 -1.000000e+12 0 2 1 0 LINESTRING (-73.17436 41.38295, -73.17465 41.3...

3313 rows × 127 columns

Exploring the dataset

What acquisition dates and rivers do our downloaded files cover?

print('Available dates are:')
print(np.unique([i[:10] for i in SWOT_HR_df['time_str']]))
print('Available rivers are:')
print(np.unique([i for i in SWOT_HR_df['river_name']]))
Available dates are:
['2023-04-08' '2023-04-09' '2023-04-10' '2023-04-11' '2023-04-12'
 'no_data']
Available rivers are:
['Allagash River' 'Androscoggin River' 'Big Black River' 'Canal de fuite'
 'Carrabassett River' 'Concord River' 'Concord River; Sudbury River'
 'Connecticut River' 'Connecticut River; Westfield River'
 'Connecticut River; White River' 'Dead River'
 'Dead River (Kennebec River)' 'Deerfield River' 'Farmington River'
 'Housatonic River' 'Ikarut River' 'Kennebec River' 'Komaktorvik River'
 'Magalloway River' 'Merrimack River' 'Moose River'
 'North Branch Penobscot River' 'Passumsic River' 'Pemigewasset River'
 'Penobscot River West Branch' 'Quinebaug River' 'Saco River'
 'Saguenay River' 'Saint Francis River' 'Saint John River'
 'Saint Lawrence River' 'Sandy River' 'Shetucket River' 'Sudbury River'
 'Thames River' 'West River' 'White River' 'no_data']

Filter dataframe by river name of interest and plot selected reaches:

Note: Some rivers have multiple names, hence using the contains function

# Enter river name
river = "Connecticut River"  # e.g. "Rhine", "Connecticut River"

## Filter dataframe
SWOT_HR_df_river = SWOT_HR_df[(SWOT_HR_df.river_name.str.contains(river))]

# Plot geopandas dataframe with 'explore' by reach id
SWOT_HR_df_river[['reach_id','river_name','geometry']].explore('reach_id', style_kwds=dict(weight=6))
Make this Notebook Trusted to load map: File -> Trust Notebook