notebooks/iceflow/corrections.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying Coordinate Transformations to Facilitate Data Comparison\n",
    "\n",
    "## Applying ITRF and Plate Movement Models to Lidar Data from Pre-IceBridge, ICESat/GLAS, IceBridge, and ICESat-2 Data Sets\n",
    "\n",
    "Comparing measurements across altimetry data sets where latitude, longitude, and surface heights have been collected over long periods of time – a span of 25+ year in the case of ICESat/GLAS, IceBridge and ICESat-2 – requires us to account for evolving coordinate reference systems and tectonic plate movement. To demonstrate why plate movement matters over time, let’s look at the Australian plate. It is moving at a rate of around 7 cm / year, so over 25 years the continent will have moved nearly 1.8 m. All of Earth’s land masses are affected by plate movement to varying degrees, and over a multi-decade time series it’s important to recognize that the position recorded for a location in 1993 is going to be different from what would be recorded at that same location today. Plate movement and updates to coordinate reference systems over time must be accounted for when analyzing data over long time periods to avoid drawing spurious conclusions from the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Important Terms and Concepts\n",
    "\n",
    "* Coordinate system: The 2-dimensional (latitude, longitude) or 3-dimensional (latitude, longitude, height) system, and unit of measure used to define where a point is located in space.\n",
    "* Datum: A modeled version of the shape of the Earth which defines the origin used to place the coordinate system in space.\n",
    "* Coordinate Reference System (CRS): A CRS comprises a coordinate system and a datum.\n",
    "* Epoch (datum epoch): An epoch is used to identify the date at which a given reference frame was realized. They are also used to date when the coordinates for a particular observation (e.g., lat, lon, elevation) were obtained. Epochs are noted in decimal year format (YYYY.DD) where numbers to the right of the decimal are derived from the Julian day-in-year. The decimal year is equal to: year + day-in-year/365 (or 366 for leap years). For example, 1 January 2010 is written as 2010.00; 5 May 1991 is the 128th Julian day-in-year and written as 1991.35. To determine what month and day is represented by, say, 2018.5, multiply the numbers to the right of the decimal by the total number of days in 2018 (.5 * 365) to determine the Julian day-in-year is 183 (after rounding up to a whole number). Referencing a Julian date calendar (e.g., [here](https://www.timeanddate.com/calendar/custom.html?year=2020&country=1&doy=1&df=1) ) we see that 2018.5 corresponds to 2 July 2018.\n",
    "* International Terrestrial Reference Frame (ITRF): A realization (i.e., a major adjustment) of the International Terrestrial Reference System (ITRS). Its origin is at the center of mass of the whole earth including the oceans and atmosphere. New ITRF solutions are produced every few years as the Earth changes and to employ the latest mathematical and surveying techniques to attempt to realize the ITRS as precisely as possible. Examples include ITRF88, ITRF93, ITRF2000, etc., with ITRF2014 being the latest realization of the ITRS. The ITRS is a formal description of the coordinate reference system, along with the mathematical and experimental procedures used to create reference frames and is maintained by the International Earth Rotation Service (IERS). The IERS issues parameters for transforming to and from older ITRF realizations. This makes it convenient to compare locations at any point on the Earth, but only for the same time, or epoch. Comparing coordinates in a dynamic frame from different epochs presents problems since the geographic features are being moved by plate tectonics, resulting in different coordinates for the same feature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Data Set\n",
    "\n",
    "For the remainder of the notebook, we will use a fictitious data set comprising latitude, longitude, and elevation coordinates. We will illustrate the use of *PROJ* and *PyProj* libraries to perform transformations between CRSs, and coordinate propagation (the shift in coordinates for a given geographic feature as a result of tectonic plate movement) between epochs in the same CRS using a Plate Motion Model. We construct a set of coordinates roughly on the Greenland ice sheet, with fictitious elevations. The particular values are used to illustrate the effects of coordinate transformation and propagation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import axes3d, Axes3D\n",
    "import numpy as np\n",
    "from pyproj import Geod, Transformer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We're going to illustrate an example of fictitious data from 1993 being transformed to be comparable to data from 2018.\n",
    "epoch1 = 1993.5\n",
    "itrf1 = \"ITRF93\"\n",
    "epoch2 = 2018.5\n",
    "itrf2 = \"ITRF2014\"\n",
    "\n",
    "# Create lats and lons evenly spaced 1x1m grid (1 cm spacing) anchored to a point on the Greenland icesheet.\n",
    "lat0 = 70.0\n",
    "lon0 = -50.0\n",
    "lat_deg_m = (1 / 111000.0)  # (~degrees / m)\n",
    "lon_deg_m = (1 / 34300.0)  # (~degrees / m)\n",
    "delta_lat = lat_deg_m * 0.01  # 1 cm\n",
    "delta_lon = lon_deg_m * 0.01  # 1 cm\n",
    "num_points = 11\n",
    "lats = lat0 + np.arange(num_points) * delta_lat\n",
    "lons = lon0 + np.arange(num_points) * delta_lon\n",
    "lon_grid, lat_grid = np.meshgrid(lons, lats)\n",
    "\n",
    "# Create elevations for illustration of coordinate transformation and propagation effects. The resulting \"surface\" will just be plane at constant elevation.\n",
    "elev_grid = np.ones(lon_grid.shape)\n",
    "\n",
    "# We set the observation epoch of the points of our sample data set.\n",
    "epoch_grid = epoch1 * np.ones(lon_grid.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we show a 10 x 10 grid spaced at 1 cm intervals, and set the elevations of each point at 1 m. The absolute latitude longitude values are at 70 N, 50 W, a spot \"on\" the Greenland ice sheet. We use this grid of points as our fictitious data set that we are going to transform using two methods. In the examples that follow, we're going to assume this data was collected in 1993 and we'd like to compare it to data collected in 2018. The fictitious data set which we will transform will be in the ITRF93 reference frame, and assume the 2018 data we'd like to compare it to is in the ITRF2014 frame.\n",
    "\n",
    "Note the scale difference on the z-axis. All axes units are meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.view_init(40, 230)\n",
    "ax.set_xlim(-0.01, 0.11)\n",
    "ax.set_ylim(-0.01, 0.11)\n",
    "\n",
    "ax.plot_wireframe((lon_grid - lon0) / lon_deg_m, (lat_grid - lat0) / lat_deg_m, elev_grid, color='blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coordinate Transformations between Reference Systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As stated above, we'd like to transform the data collected in 1993 in the ITRF93 reference frame to allow comparison with points in the ITRF2014 frame. To do so, we create a *PROJ* pipeline that consists of a series of steps:\n",
    "\n",
    "* Declare that all the ellipsoids we're dealing with are WGS84\n",
    "* Convert from degrees to radians\n",
    "* Convert from lat/lon/elevations to Cartesian coordinates\n",
    "* Convert from the ITRF93 frame to ITRF2014\n",
    "* Convert back to latitude / longitude coordinates\n",
    "* Convert from radians to degrees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_transformer(from_crs, to_crs, to_epoch=None):\n",
    "    plate_model_step = \"\"\n",
    "    if to_epoch:\n",
    "        plate_model_step = f\"+step +inv +init={to_crs}:NOAM +t_epoch={to_epoch} \"\n",
    "\n",
    "    pipeline = (\n",
    "        f\"+proj=pipeline +ellps=WGS84 \"\n",
    "        f\"+step +proj=unitconvert +xy_in=deg +xy_out=rad \"\n",
    "        f\"+step +proj=latlon \"\n",
    "        f\"+step +proj=cart \"\n",
    "        f\"+step +inv +init={to_crs}:{from_crs} \"\n",
    "        + plate_model_step +                \n",
    "        f\"+step +proj=cart +inv \"\n",
    "        f\"+step +proj=unitconvert +xy_in=rad +xy_out=deg\"\n",
    "    )\n",
    "    return Transformer.from_pipeline(pipeline)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now calculate the magnitude of the coordinate shifts from ITRF93 to ITRF2014, and then plot the original grid (blue) and the newly shifted grid (green)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = create_transformer(itrf1, itrf2)\n",
    "\n",
    "(lons_2014, lats_2014, elevs_2014, epochs_2014) = t.transform(lon_grid.T.flatten(), \n",
    "                                                              lat_grid.T.flatten(), \n",
    "                                                              elev_grid, \n",
    "                                                              epoch1 * np.ones(lon_grid.shape))\n",
    "\n",
    "max_lon_delta = np.max(abs(lons_2014 - lon_grid.T.flatten()))\n",
    "print(f\"Longitude delta (m): {max_lon_delta / lon_deg_m:0.3f}\")\n",
    "max_lat_delta = np.max(abs(lats_2014 - lat_grid.T.flatten()))\n",
    "print(f\"Latitude delta (m): {max_lat_delta / lat_deg_m:0.3f}\")\n",
    "max_elev_delta = np.max(abs(elevs_2014 - elev_grid))\n",
    "print(f\"Elevation delta (m): {max_elev_delta:0.3f}\")\n",
    "\n",
    "lons_2014 = np.reshape(lons_2014, (num_points, num_points))\n",
    "lats_2014 = np.reshape(lats_2014, (num_points, num_points))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.view_init(40, 255)\n",
    "ax.set_xlim(-0.01, 0.11)\n",
    "ax.set_ylim(-0.01, 0.11)\n",
    "ax.set_zlim(1, 1.01)\n",
    "\n",
    "ax.plot_wireframe((lon_grid - lon0) / lon_deg_m, (lat_grid - lat0) / lat_deg_m, elev_grid, color='blue')\n",
    "ax.plot_wireframe((lons_2014 - lon0) / lon_deg_m, (lats_2014 - lat0) / lat_deg_m, elevs_2014, color='green')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coordinate Propagation between Epochs within One Coordinate Reference System\n",
    "\n",
    "Now that our example data set is in the same reference frame as the data with which we'd like to compare it, we'd like to account for the shift in coordinates for a given geographic feature as a result of tectonic plate movement. To do so, we can make use of the ITRF2014 Plate Motion Model (PMM) and apply both the ITRF93 => ITRF2014 coordinate transformation _*and*_ the epoch 1993.5 => 2018.5 shift using the ITRFs PMM model parameters. \n",
    "\n",
    "Here we show the original data set (blue), the data set with only the ITRF93 => ITRF2014 transformation applied (green), and the data set obtained by applying the PMM coordinate propagation from 1993.5 through 2018.5 in five-year increments (red). We see that the ITRF93 => ITRF2014 shift applied a relatively small vertical and horizontal shift, and the PMM results in a much larger shift roughly to the northwest. We see that the magnitude of this PMM shift is much more significant, and perhaps significant to the comparison of data between epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def propagated_data(from_crs, to_crs, obs_epoch, to_epoch, lon_grid, lat_grid, elev_grid):\n",
    "    t = create_transformer(from_crs, to_crs, to_epoch)\n",
    "    (lons_new, lats_new, elevs_new, epochs_new) = t.transform(lon_grid.T.flatten(), \n",
    "                                                              lat_grid.T.flatten(), \n",
    "                                                              elev_grid, \n",
    "                                                              obs_epoch * np.ones(lon_grid.shape))\n",
    "    lons_new = np.reshape(lons_new, (num_points, num_points))\n",
    "    lats_new = np.reshape(lats_new, (num_points, num_points))\n",
    "\n",
    "    return (lons_new, lats_new, elevs_new, epochs_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.view_init(40, 255)\n",
    "ax.set_zlim(1, 1.01)\n",
    "\n",
    "ax.plot_wireframe((lon_grid - lon0) / lon_deg_m, (lat_grid - lat0) / lat_deg_m, elev_grid, color='blue')\n",
    "\n",
    "count = 6\n",
    "propd_data = {y: propagated_data(itrf1, itrf2, epoch1, y, lon_grid, lat_grid, elev_grid) for y in np.linspace(epoch1, epoch2, count)}\n",
    "print(propd_data.keys())\n",
    "\n",
    "for y, (lons, lats, elevs, _) in propd_data.items():\n",
    "    ax.plot_wireframe((lons - lon0) / lon_deg_m, (lats - lat0) / lat_deg_m, elevs, color=mpl.colors.hsv_to_rgb([1, 0.5, 0.8]))\n",
    "    \n",
    "ax.plot_wireframe((lons_2014 - lon0) / lon_deg_m, (lats_2014 - lat0) / lat_deg_m, elevs_2014, color='green')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_lon_delta = np.max(abs(propd_data[2018.5][0] - propd_data[1993.5][0]))\n",
    "print(f\"Longitude delta (m): {max_lon_delta / lon_deg_m:0.2f}\")\n",
    "\n",
    "max_lat_delta = np.max(abs(propd_data[2018.5][1] - propd_data[1993.5][1]))\n",
    "print(f\"Latitude delta (m): {max_lat_delta / lat_deg_m:0.2f}\")\n",
    "\n",
    "max_elev_delta = np.max(abs(propd_data[2018.5][2] - propd_data[1993.5][2]))\n",
    "print(f\"Elevation delta (m): {max_elev_delta:0.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "### Overview of the Reference Frame and Plate Shift Problem\n",
    "1. [The many paths to a common ground: A comparison of transformations between GDA94 and ITRF](https://www.spatial.nsw.gov.au/__data/assets/pdf_file/0003/165702/2011_Haasdyk_and_Janssen_IGNSS2011_comparison_of_GDA94-ITRF_transformations.pdf)\n",
    "1. [GDA94, ITRF & WGS84Whatâ€™s the difference?](http://www.quickclose.com.au/stanaway07pres.pdf)\n",
    "1. [Misaligned maps? High-accuracy data must become time-dependent](https://www.opengeospatial.org/blog/3045)\n",
    "1. [Dynamic Datum Transformations in Australia and New Zealand](http://ceur-ws.org/Vol-1142/paper6.pdf)\n",
    "1. [Making Sense of Evolving Datums:WGS84 and NAD83](http://www.naref.org/transf/nad83_hydroscan2006.pdf)\n",
    "1. [Re: Dynamic/static WGS84 datum problem with web-mapping and map misalignment](http://osgeo-org.1560.x6.nabble.com/Re-Dynamic-static-WGS84-datum-problem-with-web-mapping-and-map-misalignment-td5412013.html)\n",
    "1. [From static to dynamic datums: 150 years of geodetic datums in New Zealand](https://doi.org/10.1080/00288306.2015.1128451)\n",
    "1. [The new world of Earth-fixed datums](https://www.spatialsource.com.au/company-industry/new-world-earth-fixed-datums)\n",
    "1. [FIG/IAG/UN-GGIM-AP/UN-ICG/NZISTechnical Seminar Reference Frame in Practice](https://www.fig.net/resources/proceedings/2016/2016_05_reference%20frame/1_Rizos.pdf)\n",
    "\n",
    "### ITRF\n",
    "1. [ITRF](https://www.iers.org/IERS/EN/DataProducts/ITRF/itrf.html;jsessionid=04D3589399F3433E18FD3739F941D96A.live1)\n",
    "1. [IERS CONVENTIONS (2003) / IERS Technical Note No. 32](https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn32.html)\n",
    "1. [IERS CONVENTIONS (2010) / IERS Technical Note No. 36](https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html)\n",
    "1. [Analysis and results of ITRF2014 / IERS Technical Note No. 38](https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn38.html)\n",
    "1. [ITRF2008: an improved solution of the internationalterrestrial reference frame](https://doi.org/10.1007/s00190-011-0444-4)\n",
    "1. [ITRF2008 plate motion model](https://doi.org/10.1029/2011JB008930)\n",
    "1. [ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions](https://doi.org/10.1002/2016JB013098)\n",
    "1. [ITRF2014 plate motion model](https://doi.org/10.1093/gji/ggx136)\n",
    "\n",
    "### GPS\n",
    "1. [GPS Reference Systems](http://www.nbmg.unr.edu/staff/pdfs/Blewitt_Encyclopedia_of_Geodesy.html)\n",
    "1. [GPS Geodetic Reference System  WGS 84](https://www.gps.gov/multimedia/presentations/2009/09/ICG/wiley.pdf)\n",
    "1. [Reference Frames in Practice Manual. FIG Publication No. 64](http://www.fig.net/resources/publications/figpub/pub64/Figpub64.pdf)\n",
    "1. [Data collection of WGS 84 information â€” or is it?](https://www.gpsworld.com/data-collection-of-wgs-84-information-or-is-it/)\n",
    "\n",
    "### Libraries\n",
    "1. [PROJ](https://proj.org/index.html)\n",
    "1. [PyProj](https://pyproj4.github.io/pyproj/stable/index.html)\n",
    "1. [GitHub: Add plate motion model parameters to ITRF init-files](https://github.com/OSGeo/PROJ/commit/403f930355926aced5caba5bfbcc230ad152cf86)\n",
    "1. [[PROJ] Static/Dynamic Webmapping Problem version 2.0](https://lists.osgeo.org/pipermail/proj/2019-July/008703.html)\n",
    "1. [[PROJ] Seeking clarification on PROJ support for temporal transformations](https://lists.osgeo.org/pipermail/proj/2019-August/008789.html)\n",
    "1. [[PROJ] OGC blog post summarising Web-mapping misalignment problem](https://lists.osgeo.org/pipermail/proj/2019-August/008776.html)\n",
    "1. [Re: [gdal-dev] Static/Dynamic datum problems](http://osgeo-org.1560.x6.nabble.com/Re-gdal-dev-Static-Dynamic-datum-problems-td5408589.html)\n",
    "\n",
    "### Plate Motion\n",
    "1. [UNAVCO Plate Motion Calculator](https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/plate-motion-calculator.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}