Skip to content

Commit

Permalink
Merge pull request #244 from podaac/swot_kmcquillan_tutorials
Browse files Browse the repository at this point in the history 
added swot tutorial
  • Loading branch information
@cassienickles
cassienickles authored Aug 22, 2024
2 parents 2de0b32 + 887b7c5 commit 6248383
Show file tree
Hide file tree
Showing 5 changed files with 796 additions and 0 deletions.
2 changes: 2 additions & 0 deletions _quarto.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,8 @@ website:
href: notebooks/GIS/SWOT_datetime_GIS.ipynb
- section: "Transform Data"
contents:
- text: "Datum Conversion"
href: quarto_text/SWOT_Datum.qmd
- text: "River Time Series"
href: notebooks/datasets/Hydrocron_SWOT_timeseries_examples.ipynb
- text: "Lake Time Series"
Expand Down
345 changes: 345 additions & 0 deletions notebooks/datasets/SWOT_USGS_Comparison.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,345 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](https://img.shields.io/badge/Community-Contribution-%20?color=grey&labelColor=yellow)\n",
"\n",
"> From the PO.DAAC Cookbook, to access the GitHub version of the notebook, follow [this link](https://github.com/podaac/tutorials/blob/master/notebooks/datasets/SWOT_USGS_Comparison.ipynb). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Compare SWOT water surface elevation with USGS gage heights\n",
"### Datum Transformation Tutorial\n",
"\n",
"#### *Authors: Katie McQuillan and George Allen, Virginia Tech*\n",
"\n",
"## Summary \n",
"This notebook showcases how to transform the horizontal and vertical coordinates of USGS gage heights to match SWOT LakeSP water surface elevation using GDAL. \n",
"\n",
"## Requirements\n",
"\n",
"### Compute environment:\n",
"\n",
"This tutorial is written to run in the following environment:\n",
"- **Local compute environment** e.g. laptop, server: this tutorial can be run on your local machine\n",
"- **GDAL** must be installed (https://gdal.org/)\n",
"\n",
"### Import libraries"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"from osgeo import osr\n",
"import pandas as pd\n",
"import numpy as np\n",
"import math\n",
"import subprocess"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Datasets"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**1. USGS Lake/Reservoir Water Surface Elevation dataset** can be acccessed using the [DataRetrieval](https://github.com/DOI-USGS/dataretrieval-python) python module with the parameter code 62615. \n",
"\n",
"**2. SWOT LakeSP dataset** can be accessed using the [Earthaccess](https://earthaccess.readthedocs.io/en/latest/) python module or the [PO.DAAC Data Downloader](https://podaac.github.io/tutorials/quarto_text/DataSubscriberDownloader.html). \n",
"\n",
"The above data were combined outside of this tutorial into a csv file called 'usgs_swot_merged_example.csv'.\n",
"\n",
"### Datum Transformation to Compare with SWOT Data\n",
"\n",
"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. \n",
"\n",
"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. \n",
"\n",
"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. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Table 1. SWOT and USGS datum information\n",
"\n",
"| Data source | Horizontal Datum | Reference Ellipsoid | Vertical Datum | EPSG Code |\n",
"| --- | --- | --- | --- | --- |\n",
"| SWOT | ITRF2014 | WGS84 | EGM2008 | EPSG:9057+EPSG:3855 |\n",
"| USGS | NAD83 (2011) | GRS80 | NAVD88 | EPSG:6349 |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"USGS data is represented using EPSG:6349. EPSG:6349 is a compound CRS that represents NAD83 (2011) + NAVD88 height. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"COMPD_CS[\"NAD83(2011) + NAVD88 height\",\n",
" GEOGCS[\"NAD83(2011)\",\n",
" DATUM[\"NAD83_National_Spatial_Reference_System_2011\",\n",
" SPHEROID[\"GRS 1980\",6378137,298.257222101,\n",
" AUTHORITY[\"EPSG\",\"7019\"]],\n",
" AUTHORITY[\"EPSG\",\"1116\"]],\n",
" PRIMEM[\"Greenwich\",0,\n",
" AUTHORITY[\"EPSG\",\"8901\"]],\n",
" UNIT[\"degree\",0.0174532925199433,\n",
" AUTHORITY[\"EPSG\",\"9122\"]],\n",
" AXIS[\"Latitude\",NORTH],\n",
" AXIS[\"Longitude\",EAST],\n",
" AUTHORITY[\"EPSG\",\"6318\"]],\n",
" VERT_CS[\"NAVD88 height\",\n",
" VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n",
" AUTHORITY[\"EPSG\",\"5103\"]],\n",
" UNIT[\"metre\",1,\n",
" AUTHORITY[\"EPSG\",\"9001\"]],\n",
" AXIS[\"Gravity-related height\",UP],\n",
" AUTHORITY[\"EPSG\",\"5703\"]],\n",
" AUTHORITY[\"EPSG\",\"6349\"]]\n"
]
}
],
"source": [
"# Details of the the compound NAD83(2011) + NAVD88 (EPSG:6349)\n",
"src = osr.SpatialReference()\n",
"src.ImportFromEPSG(6349)\n",
"print(src)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"SWOT data is represented using EPSG:9057 + EPSG:3855. EPSG:9057 represents WGS84 (G1762) and EPSG: 3855 represents EGM 2008. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"GEOGCS[\"WGS 84 (G1762)\",\n",
" DATUM[\"World_Geodetic_System_1984_G1762\",\n",
" SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
" AUTHORITY[\"EPSG\",\"7030\"]],\n",
" AUTHORITY[\"EPSG\",\"1156\"]],\n",
" PRIMEM[\"Greenwich\",0,\n",
" AUTHORITY[\"EPSG\",\"8901\"]],\n",
" UNIT[\"degree\",0.0174532925199433,\n",
" AUTHORITY[\"EPSG\",\"9122\"]],\n",
" AXIS[\"Latitude\",NORTH],\n",
" AXIS[\"Longitude\",EAST],\n",
" AUTHORITY[\"EPSG\",\"9057\"]]\n"
]
}
],
"source": [
"# Details for EPSG:9057 WGS84 (G1762)\n",
"dst = osr.SpatialReference()\n",
"dst.ImportFromEPSG(9057)\n",
"print(dst)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"VERT_CS[\"EGM2008 height\",\n",
" VERT_DATUM[\"EGM2008 geoid\",2005,\n",
" AUTHORITY[\"EPSG\",\"1027\"]],\n",
" UNIT[\"metre\",1,\n",
" AUTHORITY[\"EPSG\",\"9001\"]],\n",
" AXIS[\"Gravity-related height\",UP],\n",
" AUTHORITY[\"EPSG\",\"3855\"]]\n"
]
}
],
"source": [
"# Details of EPSG:3855 EGM2008\n",
"v_dst = osr.SpatialReference()\n",
"v_dst.ImportFromEPSG(3855)\n",
"print(v_dst)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Transform USGS coordinates for comparison with SWOT LakeSP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Prep the data"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The number of cotemporal USGS and SWOT observations = 425\n"
]
}
],
"source": [
"# Change working directory to tutorials folder\n",
"os.chdir('.')\n",
"\n",
"# Open the dataset of cotemporal SWOT and USGS observations\n",
"swot_usgs_df = pd.read_csv(\"../resources/usgs_swot_merged_example.csv\", index_col=0)\n",
"\n",
"# How many cotemporal observations? \n",
"print('The number of cotemporal USGS and SWOT observations = ' + str(swot_usgs_df.shape[0]))\n",
"\n",
"# Get data into correct format to pass to gdal including longitude, latitude, and gage height in feet\n",
"in_gdal = swot_usgs_df[[\"usgs_long\", \"usgs_lat\", \"usgs_X_62615_00000\"]].copy()\n",
"\n",
"# Since the USGS heights are in feet but the projection we have assigned are in meters, convert heights to meters\n",
"in_gdal.loc[:,'usgs_X_62615_00000'] = in_gdal.loc[:,'usgs_X_62615_00000'] * 0.3048\n",
"\n",
"# Create a column with long, lat, height combined \n",
"in_gdal.loc[:,\"out\"] = [\n",
" str(i) + \" \" + str(j) + \" \" + str(k)\n",
" for i, j, k in zip(\n",
" in_gdal[\"usgs_long\"],\n",
" in_gdal[\"usgs_lat\"],\n",
" in_gdal[\"usgs_X_62615_00000\"],\n",
" )\n",
"]\n",
"\n",
"# Save the combined column to a text file\n",
"in_gdal[\"out\"].to_csv(\"../resources/gdal_in.txt\", header=False, index=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Transform using gdal"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"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)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Write the gdal command and run in the shell \n",
"os.chdir('../')\n",
"cd_command = \"cd \" + os.getcwd() + \"//resources && \"\n",
"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\"'\n",
"subprocess.run(cd_command + gdal_command, shell=True)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Calculate error statistics between USGS datum-transformed data and SWOT data"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" MBE (m) MAE (m) RMSE (m) One-Sigma (m)\n",
"0 -0.144761 0.299489 0.974874 0.139806\n"
]
}
],
"source": [
"# Merge back with the original data \n",
"out_gdal = pd.read_csv(\"resources/gdal_out.txt\", header=None)\n",
"out_gdal = out_gdal.rename(columns={0: \"result\"})\n",
"out_gdal[[\"usgs_long\", \"usgs_lat\", \"usgs_X_62615_00000_egm0_meters\"]] = out_gdal[\"result\"].str.split(\" \", expand=True)\n",
"swot_usgs_df[\"usgs_X_62615_00000_egm0_meters\"] = out_gdal[\"usgs_X_62615_00000_egm0_meters\"].astype(float)\n",
"\n",
"# Error stats\n",
"mae = np.mean(np.abs(np.subtract(swot_usgs_df[\"usgs_X_62615_00000_egm0_meters\"], swot_usgs_df[\"swot_wse\"])))\n",
"bias = np.mean(np.subtract(swot_usgs_df[\"usgs_X_62615_00000_egm0_meters\"], swot_usgs_df[\"swot_wse\"]))\n",
"rmse = math.sqrt(np.square(np.subtract(swot_usgs_df[\"usgs_X_62615_00000_egm0_meters\"], swot_usgs_df[\"swot_wse\"])).mean())\n",
"one_sigma = np.quantile(np.abs(np.subtract(swot_usgs_df[\"usgs_X_62615_00000_egm0_meters\"], swot_usgs_df[\"swot_wse\"])),0.68)\n",
"results = pd.DataFrame([[bias, mae, rmse, one_sigma]], columns=[\"MBE (m)\", \"MAE (m)\", \"RMSE (m)\", \"One-Sigma (m)\"])\n",
"print(results)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.13"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Loading

0 comments on commit 6248383

Please sign in to comment.