{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "87fbf7ef",
   "metadata": {},
   "source": [
    "## Geolocation Error Budget for NISAR\n",
    "\n",
    "Incidence angle = 42 degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0675459e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from mintpy.utils import readfile, utils as ut\n",
    "from mintpy.simulation import iono\n",
    "\n",
    "inc_angle = 42 / 180 * np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80e02e01",
   "metadata": {},
   "source": [
    "### Table IV - Range Geolocation Error Budget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "18f2beef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iono: residual GIM: 0.223 +/- 1.206 m\n",
      "Iono: topside  TEC: 1.739 +/- 0.517 m\n"
     ]
    }
   ],
   "source": [
    "## ionosphere\n",
    "# residual ionosphere based on the bias/uncertainty of GIM\n",
    "# reference: Hernandez et al. (2009, Table I)\n",
    "# note that we are using JHR GIM rather than JLR GIM, which from our study here at least, JHR is significantly better than JLR.\n",
    "tec_jlr = np.array([0.72, 4.49], dtype=np.float32)     # mean and STD in TECUm\n",
    "iono_delay_L = iono.vtec2range_delay(tec_jlr, inc_angle=42, freq=iono.SAR_BAND['L'])\n",
    "iono_delay_S = iono.vtec2range_delay(tec_jlr, inc_angle=42, freq=iono.SAR_BAND['S'])\n",
    "print('Iono: residual GIM: {:.3f} +/- {:.3f} m'.format(iono_delay_L[0], iono_delay_L[1]))\n",
    "\n",
    "# worse case scenario spatially based on total TEC and the map of snapshot of topTECs\n",
    "top_tec = np.array([6.6, 1.8], dtype=np.float32)\n",
    "top_tec_delay_L = iono.vtec2range_delay(top_tec, inc_angle=42, freq=iono.SAR_BAND['L'])\n",
    "top_tec_delay_S = iono.vtec2range_delay(top_tec, inc_angle=42, freq=iono.SAR_BAND['S'])\n",
    "print('Iono: topside  TEC: {:.3f} +/- {:.3f} m'.format(top_tec_delay_L[0], top_tec_delay_L[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "381aee0b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tropo: 0.040 +/- 0.050 m\n"
     ]
    }
   ],
   "source": [
    "## residual troposphere\n",
    "# mean: +/- 3 cm in zenith direction (Yu et al., 2021, Fig. 1a1)\n",
    "# STD:  +/- 5 cm in LOS    direction (Fattahi & Amelung, 2015, sec 4.2.2)\n",
    "tropo = np.array([0.03 / np.cos(inc_angle), 0.05], dtype=np.float32)\n",
    "print('Tropo: {:.3f} +/- {:.3f} m'.format(tropo[0], tropo[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a17ecc27",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Earth motion: 0.000 +/- 0.044 m\n"
     ]
    }
   ],
   "source": [
    "## uncompensated tidal effects\n",
    "OTL_ENU = np.array([0.0100, 0.0100, 0.0500], dtype=np.float32)    # reference: Martens et al. (2016); Gisinger (0.01, 0.01, 0.05)\n",
    "PT_ENU  = np.array([0.0035, 0.0035, 0.0125], dtype=np.float32)    # reference: Petit & Luzum (2010)\n",
    "ATM_ENU = np.array([0.0005, 0.0005, 0.0050], dtype=np.float32)    # reference: Dong et al. (2002)\n",
    "OTL_LOS = ut.enu2los(OTL_ENU[0], OTL_ENU[1], OTL_ENU[2], inc_angle=inc_angle*180/np.pi, az_angle=-102)\n",
    "PT_LOS  = ut.enu2los( PT_ENU[0],  PT_ENU[1],  PT_ENU[2], inc_angle=inc_angle*180/np.pi, az_angle=-102)\n",
    "ATM_LOS = ut.enu2los(ATM_ENU[0], ATM_ENU[1], ATM_ENU[2], inc_angle=inc_angle*180/np.pi, az_angle=-102)\n",
    "EM_LOS = np.array([0, np.sqrt(OTL_LOS**2 + PT_LOS**2 + ATM_LOS**2)], dtype=np.float32)\n",
    "print('Earth motion: {:.3f} +/- {:.3f} m'.format(EM_LOS[0], EM_LOS[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4a53f346",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Orbital error: 0.005 +/- 0.030 m\n"
     ]
    }
   ],
   "source": [
    "## orbital error\n",
    "# reference: Peter et al. (2017)\n",
    "orb_s1 = np.array([0.005, 0.030], dtype=np.float32)\n",
    "print('Orbital error: {:.3f} +/- {:.3f} m'.format(orb_s1[0], orb_s1[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2fac8768",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DEM error: 4.442 +/- 0.002 m\n"
     ]
    }
   ],
   "source": [
    "## DEM error\n",
    "# reference: 1. NISAR handbook;  2. Sansosti et al. (2006, Eq. 30);  3. Jung et al. (2019, Eq. 8);  4. Fahrland et al. (2020)\n",
    "# Copernicus global DEM accuracy: https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198\n",
    "rg = 980e3       # m, NISAR\n",
    "bperp_max = 350  # m\n",
    "dem_err = 4      # m, absolute vertical accuracy\n",
    "dem_LOS_rel = dem_err * bperp_max / (rg * np.sin(inc_angle)) * 1\n",
    "dem_LOS_abs = dem_err / np.tan(inc_angle)\n",
    "dem_LOS = np.array([dem_LOS_abs, dem_LOS_rel], dtype=np.float32)\n",
    "print('DEM error: {:.3f} +/- {:.3f} m'.format(dem_LOS[0], dem_LOS[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2a9191cb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NISAR L-band:\n",
      "  Iono residual GIM: 0.223 +/- 1.206 m\n",
      "  Iono topside  TEC: 1.739 +/- 0.517 m\n",
      "  Tropo            : 0.040 +/- 0.050 m\n",
      "  Tidal residual   : 0.000 +/- 0.044 m\n",
      "  Orbital error    : 0.005 +/- 0.030 m\n",
      "  DEM error        : 4.442 +/- 0.002 m\n",
      "  Overall (L-band) : 6.449 +/- 1.314 m\n",
      "NISAR S-band:\n",
      "  Iono residual GIM: 0.037 +/- 0.215 m\n",
      "  Iono topside  TEC: 0.307 +/- 0.091 m\n",
      "  ...\n",
      "  Overall (S-band) : 4.832 +/- 0.244 m\n"
     ]
    }
   ],
   "source": [
    "geoloc_L = np.sqrt(iono_delay_L**2 + top_tec_delay_L**2 + tropo**2 + EM_LOS**2 + orb_s1**2 + dem_LOS**2)\n",
    "geoloc_L[0] = iono_delay_L[0] + top_tec_delay_L[0] + tropo[0] + EM_LOS[0] + orb_s1[0] + dem_LOS[0]\n",
    "print('NISAR L-band:')\n",
    "print('  Iono residual GIM: {:.3f} +/- {:.3f} m'.format(iono_delay_L[0], iono_delay_L[1]))\n",
    "print('  Iono topside  TEC: {:.3f} +/- {:.3f} m'.format(top_tec_delay_L[0], top_tec_delay_L[1]))\n",
    "print('  Tropo            : {:.3f} +/- {:.3f} m'.format(tropo[0], tropo[1]))\n",
    "print('  Tidal residual   : {:.3f} +/- {:.3f} m'.format(EM_LOS[0], EM_LOS[1]))\n",
    "print('  Orbital error    : {:.3f} +/- {:.3f} m'.format(orb_s1[0], orb_s1[1]))\n",
    "print('  DEM error        : {:.3f} +/- {:.3f} m'.format(dem_LOS[0], dem_LOS[1]))\n",
    "print('  Overall (L-band) : {:.3f} +/- {:.3f} m'.format(geoloc_L[0], geoloc_L[1]))\n",
    "\n",
    "geoloc_S = np.sqrt(iono_delay_S**2 + top_tec_delay_S**2 + tropo**2 + EM_LOS**2 + orb_s1**2 + dem_LOS**2)\n",
    "geoloc_S[0] = iono_delay_S[0] + top_tec_delay_S[0] + tropo[0] + EM_LOS[0] + orb_s1[0] + dem_LOS[0]\n",
    "print('NISAR S-band:')\n",
    "print('  Iono residual GIM: {:.3f} +/- {:.3f} m'.format(iono_delay_S[0], iono_delay_S[1]))\n",
    "print('  Iono topside  TEC: {:.3f} +/- {:.3f} m'.format(top_tec_delay_S[0], top_tec_delay_S[1]))\n",
    "print('  ...')\n",
    "print('  Overall (S-band) : {:.3f} +/- {:.3f} m'.format(geoloc_S[0], geoloc_S[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f97a5156",
   "metadata": {},
   "source": [
    "### OPERA NI-CSLC: Relative Geolocation Accuracy (using ALOS-2 as a proxy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3fa2c782",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range  : 0.560 m --> 0.09 pixel\n",
      "Azimuth: 0.810 m --> 0.12 pixel\n"
     ]
    }
   ],
   "source": [
    "## range\n",
    "rg_reso = 6.25   # m for NISAR-L 24 MHz\n",
    "rg_rmse = 0.56   # m from ALOS-2 desc track 23 over Kyushu, Japan\n",
    "\n",
    "## azimuth\n",
    "az_reso = 7  # m for NISAR-L (NISAR handbook)\n",
    "mintpy_dir = os.path.expanduser('~/data/geolocation/KyushuAlos2DT23/mintpy_offset')\n",
    "ts_file   = os.path.join(mintpy_dir, 'timeseriesAz.h5')  # timeseriesRg.h5 or timeseriesAz.h5\n",
    "\n",
    "# calculate az geoloc rmse\n",
    "mask_file = os.path.join(mintpy_dir, 'maskResInv.h5')\n",
    "ts_data = readfile.read(ts_file)[0];  ts_data = ts_data.reshape(ts_data.shape[0], -1)\n",
    "mask = readfile.read(mask_file)[0].flatten()\n",
    "ts_data[:, mask == 0] = np.nan\n",
    "ts_med = np.nanmedian(ts_data, axis=-1)\n",
    "ts_med -= np.nanmedian(ts_med)\n",
    "az_rmse = ut.root_mean_sq_error(ts_med)\n",
    "\n",
    "print('Range  : {:.3f} m --> {:.2f} pixel'.format(rg_rmse, rg_rmse / rg_reso))\n",
    "print('Azimuth: {:.3f} m --> {:.2f} pixel'.format(az_rmse, az_rmse / az_reso))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18820d37",
   "metadata": {},
   "source": [
    "### OPERA NI-CSLC: Absolute Geolocation Accuracy in Range Direction (w/o DEM error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "544b7bc3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range w/o Topo: 1.8 m\n"
     ]
    }
   ],
   "source": [
    "geoloc = np.sqrt(iono_delay_L**2 + top_tec_delay_L**2 + tropo**2 + EM_LOS**2 + orb_s1**2)\n",
    "print('Range w/o Topo: {:.1f} m'.format(geoloc[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc20bb09",
   "metadata": {},
   "source": [
    "### OPERA NI-CSLC: Absolute Geolocation Accuracy in Azimuth Direction (w/o DEM error)\n",
    "\n",
    "We do NOT know due to the lack of TEC gradient/slope knowledge!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c719734",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}