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

Compare SWOT water surface elevation with USGS gage heights

Datum Transformation Tutorial

Authors: Katie McQuillan and George Allen, Virginia Tech

Summary

This notebook showcases how to transform the horizontal and vertical coordinates of USGS gage heights to match SWOT LakeSP water surface elevation using GDAL.

Requirements

Compute environment:

This tutorial is written to run in the following environment: - Local compute environment e.g. laptop, server: this tutorial can be run on your local machine - GDAL must be installed (https://gdal.org/)

Import libraries

import os
from osgeo import osr
import pandas as pd
import numpy as np
import math
import subprocess

Datasets

1. USGS Lake/Reservoir Water Surface Elevation dataset can be acccessed using the DataRetrieval python module with the parameter code 62615.

2. SWOT LakeSP dataset can be accessed using the Earthaccess python module or the PO.DAAC Data Downloader.

The above data were combined outside of this tutorial into a csv file called ‘usgs_swot_merged_example.csv’.

Datum Transformation to Compare with SWOT Data

Cotemporal SWOT LakeSP and USGS observations were used to form the dataset used for comparisons including X lakes with gages. To directly compare SWOT and USGS datasets, the USGS horizontal and vertical coordiantes must be transformed to match the SWOT datums. Datums for each dataset are noted in Table 1.

It is important to note that using the generic WGS84 (EPSG:4326) is not recommended because it is based on a datum ensemble whose positional accuracy is approximately two meters. Instead, a realization such as WGS84 (G1762) is recommended. WGS84 (G1762) and ITRF 2014 are equivalent for all practical purposes when their epochs are the same.

Epochs are used to define a reference date for positions esablished using the datum ellipsoid and reference frame. Due to tectonic plate movement, landmasses are constantly moving in relationship to each other and in relation to the reference frame. Therefore, accounting for the epoch is necessary for accurate coordinate transformations. The NAD83 (2011) epoch is 2010.0. The standard epoch of WGS84 (G1762) is 2005.0 and the standard epoch of ITRF2014 is 2010.0. Since SWOT is based on ITRF2014, we set the target epoch to 2010.0.

Table 1. SWOT and USGS datum information

Data source	Horizontal Datum	Reference Ellipsoid	Vertical Datum	EPSG Code
SWOT	ITRF2014	WGS84	EGM2008	EPSG:9057+EPSG:3855
USGS	NAD83 (2011)	GRS80	NAVD88	EPSG:6349

USGS data is represented using EPSG:6349. EPSG:6349 is a compound CRS that represents NAD83 (2011) + NAVD88 height.

# Details of the the compound NAD83(2011) + NAVD88 (EPSG:6349)
src = osr.SpatialReference()
src.ImportFromEPSG(6349)
print(src)

COMPD_CS["NAD83(2011) + NAVD88 height",
    GEOGCS["NAD83(2011)",
        DATUM["NAD83_National_Spatial_Reference_System_2011",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","1116"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AXIS["Latitude",NORTH],
        AXIS["Longitude",EAST],
        AUTHORITY["EPSG","6318"]],
    VERT_CS["NAVD88 height",
        VERT_DATUM["North American Vertical Datum 1988",2005,
            AUTHORITY["EPSG","5103"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Gravity-related height",UP],
        AUTHORITY["EPSG","5703"]],
    AUTHORITY["EPSG","6349"]]

SWOT data is represented using EPSG:9057 + EPSG:3855. EPSG:9057 represents WGS84 (G1762) and EPSG: 3855 represents EGM 2008.

# Details for EPSG:9057 WGS84 (G1762)
dst = osr.SpatialReference()
dst.ImportFromEPSG(9057)
print(dst)

GEOGCS["WGS 84 (G1762)",
    DATUM["World_Geodetic_System_1984_G1762",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","1156"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","9057"]]

# Details of EPSG:3855 EGM2008
v_dst = osr.SpatialReference()
v_dst.ImportFromEPSG(3855)
print(v_dst)

VERT_CS["EGM2008 height",
    VERT_DATUM["EGM2008 geoid",2005,
        AUTHORITY["EPSG","1027"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Gravity-related height",UP],
    AUTHORITY["EPSG","3855"]]

Transform USGS coordinates for comparison with SWOT LakeSP

Prep the data

# Change working directory to tutorials folder
os.chdir('.')

# Open the dataset of cotemporal SWOT and USGS observations
swot_usgs_df = pd.read_csv("../resources/usgs_swot_merged_example.csv", index_col=0)

# How many cotemporal observations? 
print('The number of cotemporal USGS and SWOT observations = ' + str(swot_usgs_df.shape[0]))

# Get data into correct format to pass to gdal including longitude, latitude, and gage height in feet
in_gdal = swot_usgs_df[["usgs_long", "usgs_lat", "usgs_X_62615_00000"]].copy()

# Since the USGS heights are in feet but the projection we have assigned are in meters, convert heights to meters
in_gdal.loc[:,'usgs_X_62615_00000'] = in_gdal.loc[:,'usgs_X_62615_00000'] * 0.3048

# Create a column with long, lat, height combined 
in_gdal.loc[:,"out"] = [
    str(i) + " " + str(j) + " " + str(k)
    for i, j, k in zip(
        in_gdal["usgs_long"],
        in_gdal["usgs_lat"],
        in_gdal["usgs_X_62615_00000"],
    )
]

# Save the combined column to a text file
in_gdal["out"].to_csv("../resources/gdal_in.txt", header=False, index=False)

The number of cotemporal USGS and SWOT observations = 425

Transform using gdal

# Write the gdal command and run in the shell 
os.chdir('../')
cd_command = "cd " + os.getcwd() + "//resources && "
gdal_command = 'gdaltransform -s_coord_epoch "2010.0" -t_coord_epoch "2010.0" -s_srs "EPSG:6349" -t_srs "EPSG:9057+EPSG:3855" < "gdal_in.txt" > "gdal_out.txt"'
subprocess.run(cd_command + gdal_command, shell=True)

CompletedProcess(args='cd /home/jovyan/tutorials/notebooks//resources && gdaltransform -s_coord_epoch "2010.0" -t_coord_epoch "2010.0" -s_srs "EPSG:6349" -t_srs "EPSG:9057+EPSG:3855" < "gdal_in.txt" > "gdal_out.txt"', returncode=0)

Calculate error statistics between USGS datum-transformed data and SWOT data

# Merge back with the original data 
out_gdal = pd.read_csv("resources/gdal_out.txt", header=None)
out_gdal = out_gdal.rename(columns={0: "result"})
out_gdal[["usgs_long", "usgs_lat", "usgs_X_62615_00000_egm0_meters"]] = out_gdal["result"].str.split(" ", expand=True)
swot_usgs_df["usgs_X_62615_00000_egm0_meters"] = out_gdal["usgs_X_62615_00000_egm0_meters"].astype(float)

# Error stats
mae = np.mean(np.abs(np.subtract(swot_usgs_df["usgs_X_62615_00000_egm0_meters"], swot_usgs_df["swot_wse"])))
bias = np.mean(np.subtract(swot_usgs_df["usgs_X_62615_00000_egm0_meters"], swot_usgs_df["swot_wse"]))
rmse = math.sqrt(np.square(np.subtract(swot_usgs_df["usgs_X_62615_00000_egm0_meters"], swot_usgs_df["swot_wse"])).mean())
one_sigma = np.quantile(np.abs(np.subtract(swot_usgs_df["usgs_X_62615_00000_egm0_meters"], swot_usgs_df["swot_wse"])),0.68)
results = pd.DataFrame([[bias, mae, rmse, one_sigma]], columns=["MBE (m)", "MAE (m)", "RMSE (m)", "One-Sigma (m)"])
print(results)

    MBE (m)   MAE (m)  RMSE (m)  One-Sigma (m)
0 -0.144761  0.299489  0.974874       0.139806