import requests
import json
import pprint
import xarray as xr
import geojson
import matplotlib.pyplot as plt
import contextily as cx
import plotly.graph_objects as go
from IPython.display import JSON, Image
import earthaccessFrom the PO.DAAC Cookbook, to access the GitHub version of the notebook, follow this link.
SWORD River Demo
Summary
This notebook will query the SWOT River Database (SWORD) for river reaches (segments) or nodes (points) and visualize results using PO.DAAC’s Feature Translation Service. We use geospatial coordinates of the queried features (here river reaches along the Kasai River, a tributary of the Congo River in Africa) to then query and download or access via the cloud Pre-SWOT Hydrology data along the specified regions. This is a programmatic approach to using the ‘Advanced Search -> River Reach’ query in the Earthdata Search GUI.
Requirements
1. Compute environment
This tutorial can be run in the following environments: - Local compute environment e.g. laptop, server: this tutorial can be run on your local machine - AWS instance running in us-west-2: NASA Earthdata Cloud data in S3 can be directly accessed via temporary credentials; this access is limited to requests made within the US West (Oregon) (code: us-west-2) AWS region.
2. 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.
Import libraries
Authenticate Earthdata Login via earthaccess
auth = earthaccess.login()Define a Function to Query Easier
First, we define a function to query a single reach, multiple reaches, or river nodes. This function returns the properties of river reaches including name, length, coordinates as well as individual node properties.
def response_to_FeatureCollection(response):
"""
This function will return a geojson.FeatureCollection representation of the features found in the provided response.
Parameters
----------
response : requests.Response - a Response object returned from a GET request on the rivers or nodes endpoint.
Returns
-------
geojson.FeatureCollection - FeatureCollection containing all features extracted from the response.
"""
featureList = []
for reach_id, reach_json in response.json()['results'].items():
reach_feature = geojson.loads(json.dumps(reach_json['geojson']))
reach_feature['properties']={k:v for k,v in reach_json.items() if k not in ['geojson', 'geometry']}
featureList.append(reach_feature)
featureCollection = geojson.FeatureCollection(featureList)
return featureCollectionSingle River Reach from the SWORD Database
In this section, we query the Feature Translation Service (FTS) SWORD service using a single Reach ID (from SWORD) to give us all of the metadata regarding the particular reach. In this example, we use the river Reach ID 13227000061. This ID represents a specific reach along the Kasai River, a tributary of the Congo River in Africa.
#change the reach ID in the link below for a different location.
response_reach = requests.get("https://fts.podaac.earthdata.nasa.gov/rivers/reach/13227000061")
featureCollection_reach = response_to_FeatureCollection(response_reach)
pprint.pprint(response_reach.json(), compact=True, width=60, depth=2){'hits': 1,
'results': {'13227000061': {...}},
'search on': {'exact': False,
'page_number': 1,
'page_size': 100,
'parameter': 'reach'},
'status': '200 OK',
'time': '6.2 ms.'}Note that we haven’t looked at any data from a collection yet, we’ve simply found the geospatial coordinates of our river reach of interest within the metadata.
Query Data by Coordinates
We can use results obtained from the FTS query to then directly and automatically query data using the earthaccess python library. We will use the coordinate information of a single reach to search for granules (files) available through the SWOT L2 raster files, whose data has the short-name SWOT_L2_HR_Raster_2.0.
#obtain the associated lat/lons for the feature
lats = [xy[1] for feature in featureCollection_reach['features'] for xy in feature['coordinates']]
lons = [xy[0] for feature in featureCollection_reach['features'] for xy in feature['coordinates']]
# find max and min of lat and lon derived when visualizing the reach for the bounding box input to earthaccess
maxlat, maxlon, minlat, minlon = max(lats), max(lons), min(lats), min(lons)#earthaccess data search
results = earthaccess.search_data(short_name = 'SWOT_L2_HR_Raster_2.0', bounding_box=(minlon,minlat,maxlon,maxlat))Granules found: 30Results contains a link to the data file (granule) from the SWOT_L2_HR_RiverSP_2.0 data collection that overlaps the geospatial search from FTS-SWORD for the river reaches of interest, which we can then download or use within the cloud. Earthaccess will return the link required for the environment (https if local and s3 if in the cloud).
Accessing Data in the Cloud (skip this chunk if local)
This code chunk will only work if you are running this script in the cloud, AWS us-west-2 region.
#access data within the cloud and open into an xarray dataset
ds_SWOT_raster = xr.open_mfdataset(earthaccess.open([results[0]]), engine="h5netcdf")
ds_SWOT_rasterOpening 1 granules, approx size: 0.01 GB<xarray.Dataset>
Dimensions: (x: 581, y: 582)
Coordinates:
* x (x) float64 3.585e+05 3.588e+05 ... 5.035e+05
* y (y) float64 9.358e+06 9.358e+06 ... 9.503e+06
Data variables: (12/39)
crs object ...
longitude (y, x) float64 dask.array<chunksize=(582, 581), meta=np.ndarray>
latitude (y, x) float64 dask.array<chunksize=(582, 581), meta=np.ndarray>
wse (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
wse_qual (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
wse_qual_bitwise (y, x) float64 dask.array<chunksize=(582, 581), meta=np.ndarray>
... ...
load_tide_fes (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
load_tide_got (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
pole_tide (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
model_dry_tropo_cor (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
model_wet_tropo_cor (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
iono_cor_gim_ka (y, x) float32 dask.array<chunksize=(582, 581), meta=np.ndarray>
Attributes: (12/49)
Conventions: CF-1.7
title: Level 2 KaRIn High Rate Raster Data Product
source: Ka-band radar interferometer
history: 2023-11-30T07:08:40Z : Creation
platform: SWOT
references: V1.1.1
... ...
x_min: 358500.0
x_max: 503500.0
y_min: 9357750.0
y_max: 9503000.0
institution: CNES
product_version: 01Downloading Data on a Local Machine
#download data into folder on local machine
earthaccess.download([results[0]], "datasets/data_downloads") Getting 1 granules, approx download size: 0.01 GB['datasets\\data_downloads\\SWOT_L2_HR_Raster_250m_UTM34M_N_x_x_x_007_055_073F_20231125T084630_20231125T084651_PIC0_01.nc']#open dataset for visualization
ds_SWOT_raster = xr.open_mfdataset(f'datasets/data_downloads/SWOT_L2_HR_Raster*', engine='h5netcdf')
ds_SWOT_raster<xarray.Dataset>
Dimensions: (x: 1920, y: 2101)
Coordinates:
* x (x) float64 2.969e+05 2.97e+05 ... 5.035e+05
* y (y) float64 4.274e+06 4.274e+06 ... 9.503e+06
Data variables: (12/39)
crs (y, x) object b'1' b'1' b'1' ... b'1' b'1' b'1'
longitude (y, x) float64 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
latitude (y, x) float64 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
wse (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
wse_qual (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
wse_qual_bitwise (y, x) float64 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
... ...
load_tide_fes (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
load_tide_got (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
pole_tide (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
model_dry_tropo_cor (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
model_wet_tropo_cor (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
iono_cor_gim_ka (y, x) float32 dask.array<chunksize=(1519, 1920), meta=np.ndarray>
Attributes: (12/49)
Conventions: CF-1.7
title: Level 2 KaRIn High Rate Raster Data Product
source: Ka-band radar interferometer
history: 2023-12-03T08:26:57Z : Creation
platform: SWOT
references: V1.1.1
... ...
x_min: 296900.0
x_max: 448800.0
y_min: 4274000.0
y_max: 4425800.0
institution: CNES
product_version: 04… Continue to Science!
Other Options to query the SWORD Database:
Multiple River Reaches
We can query the FTS SWORD service over multiple river reaches too. In this example, we use ID 132270000. This ID represents multiple reaches along the Kasai River, a tributary of the Congo River in Africa. Note that this reach includes the reach ID of 13227000061 we plotted earlier. The response includes 9 individual reaches.
response_multi = requests.get("https://fts.podaac.earthdata.nasa.gov/rivers/reach/132270000")
featureCollection_multi = response_to_FeatureCollection(response_multi)
pprint.pprint(response_multi.json(), compact=True, width=60, depth=2){'hits': 9,
'results': {'13227000011': {...},
'13227000021': {...},
'13227000031': {...},
'13227000041': {...},
'13227000051': {...},
'13227000061': {...},
'13227000071': {...},
'13227000081': {...},
'13227000091': {...}},
'search on': {'exact': False,
'page_number': 1,
'page_size': 100,
'parameter': 'reach'},
'status': '200 OK',
'time': '58.247 ms.'}River Nodes within a Reach
We can also query the FTS SWORD service for river nodes, points along the river reaches. In this example, we use the ID of 13227000060. Note that this ID corresponds to the same reach ID 13227000061 we used earlier. The only difference is that the last digit of 0 corresponds to all the individual nodes along the reach. If the last digit is 1 in the reach ID, it corresponds to the overall properties of the reach itself. The response returns 52 nodes along the reach.
response = requests.get("https://fts.podaac.earthdata.nasa.gov/rivers/node/13227000060")
featureCollection = response_to_FeatureCollection(response)
pprint.pprint(response.json(), compact=True, width=60, depth=2){'hits': 52,
'results': {'13227000060011': {...},
'13227000060021': {...},
'13227000060031': {...},
'13227000060041': {...},
'13227000060051': {...},
'13227000060061': {...},
'13227000060071': {...},
'13227000060081': {...},
'13227000060091': {...},
'13227000060101': {...},
'13227000060111': {...},
'13227000060121': {...},
'13227000060131': {...},
'13227000060141': {...},
'13227000060151': {...},
'13227000060161': {...},
'13227000060171': {...},
'13227000060181': {...},
'13227000060191': {...},
'13227000060201': {...},
'13227000060211': {...},
'13227000060221': {...},
'13227000060231': {...},
'13227000060241': {...},
'13227000060251': {...},
'13227000060261': {...},
'13227000060271': {...},
'13227000060281': {...},
'13227000060291': {...},
'13227000060301': {...},
'13227000060311': {...},
'13227000060321': {...},
'13227000060331': {...},
'13227000060341': {...},
'13227000060351': {...},
'13227000060361': {...},
'13227000060371': {...},
'13227000060381': {...},
'13227000060391': {...},
'13227000060401': {...},
'13227000060411': {...},
'13227000060421': {...},
'13227000060431': {...},
'13227000060441': {...},
'13227000060451': {...},
'13227000060461': {...},
'13227000060471': {...},
'13227000060481': {...},
'13227000060491': {...},
'13227000060501': {...},
'13227000060511': {...},
'13227000060521': {...}},
'search on': {'exact': False,
'page_number': 1,
'page_size': 100,
'parameter': 'node'},
'status': '200 OK',
'time': '119.271 ms.'}