-
Notifications
You must be signed in to change notification settings - Fork 42
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
added usgs transform to swot tutorial and dataset
- Loading branch information
Showing
2 changed files
with
768 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,342 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Compare SWOT water surface elevation with USGS gage heights\n", | ||
| "\n", | ||
| "#### *Author: 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", | ||
| "### 1. 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", | ||
| "\n", | ||
| "- GDAL must be installed (https://gdal.org/)\n", | ||
| "\n", | ||
| "### 2. Earthdata Login\n", | ||
| "\n", | ||
| "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. \n", | ||
| "\n", | ||
| "\n", | ||
| "### Import libraries" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 2, | ||
| "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 python module with the parameter code 62615. \n", | ||
| "\n", | ||
| "2. SWOT LakeSP dataset can be accessed using the EarthAccess python module or the PO.DAAC Data Downloader. \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", | ||
| "| 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" | ||
| ] | ||
| }, | ||
| { | ||
| "name": "stderr", | ||
| "output_type": "stream", | ||
| "text": [ | ||
| "c:\\Users\\kmcquil\\anaconda3\\envs\\swot-wind\\lib\\site-packages\\osgeo\\osr.py:410: FutureWarning: Neither osr.UseExceptions() nor osr.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.\n", | ||
| " warnings.warn(\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": 25, | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "name": "stdout", | ||
| "output_type": "stream", | ||
| "text": [ | ||
| "The number of cotemporal USGS and SWOT observations = 425\n" | ||
| ] | ||
| }, | ||
| { | ||
| "data": { | ||
| "text/plain": [ | ||
| "CompletedProcess(args='cd c:\\\\Users\\\\kmcquil\\\\Documents\\\\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": 25, | ||
| "metadata": {}, | ||
| "output_type": "execute_result" | ||
| } | ||
| ], | ||
| "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(\"notebooks/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(\"notebooks/resources/gdal_in.txt\", header=False, index=False)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "#### Transform using gdal" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Write the gdal command and run in the shell \n", | ||
| "cd_command = \"cd \" + os.getcwd() + \"\\\\notebooks\\\\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" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 28, | ||
| "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(\"notebooks/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": "swot-wind", | ||
| "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.9.18" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 2 | ||
| } |
Oops, something went wrong.