wip: cleaning up notebook

cleanup notebook

Author: Sammy Guo

python/gtsam/examples/DifferentialPseudorangeExample.ipynb

@@ -11,7 +11,7 @@
     "## Overview\n",
     "Differential GNSS (DGNSS) improves positioning accuracy by using a reference receiver at a precisely known location to correct errors affecting nearby GNSS users. The reference station receives the same satellite signals as the user and compares the measured ranges to the true ranges based on its known position. The differences represent combined errors from satellite clocks, orbits, and atmospheric delays. These corrections are transmitted to the user receiver, which applies them to its own measurements. Because both receivers experience nearly the same errors, especially when close together, many of these errors cancel out, reducing position errors from several meters to sub-meter or better accuracy.\n",
     "\n",
-    "GNSS positioning problems can be encoded as factor graphs in GTSAM; this notebook extends the groundwork from [SinglePointPositioningExample.ipynb](https://github.com/borglab/gtsam/blob/develop/python/gtsam/examples/SinglePointPositioningExample.ipynb) to improve receiver positioning accuracy using corrections from a nearby reference receiver."
+    "GNSS positioning problems can be encoded as factor graphs in GTSAM; this notebook extends the groundwork from [SinglePointPositioningExample.ipynb](https://borglab.github.io/gtsam/singlepointpositioningexample/) to improve receiver positioning accuracy using corrections from a nearby reference receiver."
    ]
   },
   {
@@ -43,10 +43,12 @@
    "source": [
     "try:\n",
     "    import google.colab\n",
-    "    %pip install --quiet pyrtklib numpy gtsam-develop pyproj tabulate\n",
+    "    %pip install --quiet pyrtklib numpy gtsam-develop pyproj tabulate folium\n",
+    "    !wget https://raw.githubusercontent.com/borglab/gtsam/refs/heads/develop/python/gtsam/examples/gnss_utils.py\n",
     "except ImportError:\n",
     "    pass\n",
     "\n",
+    "import folium\n",
     "import gnss_utils\n",
     "import gtsam\n",
     "from gtsam.symbol_shorthand import B, C, X\n",
@@ -165,7 +167,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "So far the solution passes the eye-ball test: It's in the right US state, in roughly the correct part of the SF Bay Area. But we can quantify error by comparing the [actual position of P222](https://www.ngs.noaa.gov/cgi-cors/CorsSidebarSelect.prl?site=P222&option=Coordinates20):\n",
+    "So far this solution passes the eye-ball test: It's in the right US state, in roughly the correct part of the SF Bay Area. But we can quantify error by comparing the [actual position of P222](https://www.ngs.noaa.gov/cgi-cors/CorsSidebarSelect.prl?site=P222&option=Coordinates20):\n",
     "```\n",
     "| ITRF2020 POSITION (EPOCH 2020.0)                                            |\n",
     "| Computed in Apr 2025 using data through gpswk 2237.                         |\n",
@@ -204,43 +206,9 @@
     "```\n",
     "In a factor graph framework, these survey coordinates correspond to a stiff `PriorFactor` constraint on ZOA1's position node to keep it anchored in the global Earth frame. The rest of the factor graph is constructed as follows:\n",
     "\n",
-    "```mermaid\n",
-    "graph TD\n",
-    "    %% Variable nodes\n",
-    "    c1((Differential Correction))\n",
-    "    \n",
-    "    subgraph P222\n",
-    "    direction LR\n",
-    "    x1((P222 Position))\n",
-    "    b1((P222 Clock Bias))\n",
-    "    end\n",
-    "\n",
-    "    subgraph ZOA1\n",
-    "    direction LR\n",
-    "    x2((ZOA1 Position))\n",
-    "    b2((ZOA1 Clock Bias))\n",
-    "    end\n",
-    "\n",
-    "    %% Factor nodes\n",
-    "    p1[PriorFactor]\n",
-    "    dpr1[DifferentialPseudorangeFactor]\n",
-    "    dpr2[DifferentialPseudorangeFactor]\n",
-    "\n",
-    "    %% Connect ZOA1's prior\n",
-    "    p1 --- x2\n",
+    "![illustrative factor graph example showcasing differential correction factors](https://raw.githubusercontent.com/masoug/gtsam-gnss/refs/heads/feature/differential-pseudorange/python/gtsam/examples/images/dpr_fg.jpg)\n",
     "\n",
-    "    %% Connect ZOA1's pseudorange\n",
-    "    dpr1 --- c1\n",
-    "    dpr1 --- x2\n",
-    "    dpr1 --- b2\n",
-    "\n",
-    "    %% Connect P222's pseudorange\n",
-    "    dpr2 --- c1\n",
-    "    dpr2 --- x1\n",
-    "    dpr2 --- b1\n",
-    "```\n",
-    "\n",
-    "The bottom row are variable nodes, and the top row are factors. A `DifferentialPseudorangeFactor` is added for each receiver observation, and their correction nodes are indexed by satellite ID. For example, ZOA1's pseudorange factor of satellite 25 links to the same correction node for P222's pseudorange observation of satellite 25. Since different satellites are in different positions in the sky, each satellite has their own correction node to account for spatial variation in atmospheric effects. Furthermore, for long-time-series observations, the correction nodes must also account for temporal variation in the atmosphere. In that case, a linear chain of correction nodes linked with between-factors is recommended. The above structure is encoded into GTSAM below:"
+    "A `DifferentialPseudorangeFactor` is added for each receiver observation, and their correction nodes are indexed by satellite ID. For example, ZOA1's pseudorange factor of satellite 25 links to the same correction node for P222's pseudorange observation of satellite 25. Since different satellites are in different positions in the sky, each satellite has their own correction node to account for spatial variation in atmospheric effects. Furthermore, for long-time-series observations, the correction nodes must also account for temporal variation in the atmosphere. In that case, a linear chain of correction nodes linked with between-factors is recommended. The above structure is encoded into GTSAM below:"
    ]
   },
   {
@@ -305,19 +273,76 @@
     "\n",
     "# Process the results:\n",
     "clock_bias = result.atDouble(p222_bias_key)\n",
-    "receiver_position = result.atVector(p222_pos_key)\n",
-    "\n",
+    "corrected_receiver_position = result.atVector(p222_pos_key)\n",
+    "print(f\"Corrected P222 position: {corrected_receiver_position} ECEF meters\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's plot these coordinates on a map, including the reference stations and single-point positions, to get a big-picture overview of how the positioning results stack up:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "# Convenience plot:\n",
+    "lon, lat, alt = ecef2lla.transform(corrected_receiver_position[0], corrected_receiver_position[1], corrected_receiver_position[2])\n",
+    "print(f\"Corrected geodetic coordinates: {lat}, {lon}, {alt} m\")\n",
+    "\n",
+    "# Center the map on P222\n",
+    "m = folium.Map(\n",
+    "    location=[lat, lon],\n",
+    "    zoom_start=17\n",
+    ")\n",
+    "\n",
+    "# Plot surveyed P222 position:\n",
+    "folium.Circle(\n",
+    "    location=[37.5392407028, -122.0832682528],\n",
+    "    radius=10,\n",
+    "    tooltip=\"P222 Truth\",\n",
+    "    popup=\"Actual P222 position\"\n",
+    ").add_to(m)\n",
+    "\n",
+    "# Plot the corrected receiver position:\n",
+    "folium.Marker(\n",
+    "    location=[lat, lon],\n",
+    "    tooltip=\"P222 Corrected\",\n",
+    "    popup=\"Estimated P222 position using corrections from ZOA1\",\n",
+    "    icon=folium.Icon(color=\"green\", icon=\"plus\")\n",
+    ").add_to(m)\n",
+    "\n",
+    "# Plot the single-point-positioning solution:\n",
     "lon, lat, alt = ecef2lla.transform(receiver_position[0], receiver_position[1], receiver_position[2])\n",
-    "print(f\"geodetic: {lat}, {lon}, {alt} m\")\n",
-    "gnss_utils.plotMap(lat, lon)"
+    "folium.Marker(\n",
+    "    location=[lat, lon],\n",
+    "    tooltip=\"P222 SPP\",\n",
+    "    popup=\"Estimated P222 position without differential corrections\",\n",
+    "    icon=folium.Icon(color=\"orange\", icon=\"minus\")\n",
+    ").add_to(m)\n",
+    "\n",
+    "# Plot ZOA1's position:\n",
+    "folium.Marker(\n",
+    "    location=[37.5430544917, -122.0159500139],\n",
+    "    tooltip=\"ZOA1\",\n",
+    "    popup=\"ZOA1 reference station surveyed position\",\n",
+    "    icon=folium.Icon(color=\"blue\")\n",
+    ").add_to(m)\n",
+    "\n",
+    "m"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We got an updated position; so let's see if it's closer to ground truth?"
+    "Scroll and zoom the map to get an idea for the relative positioning error between the single-point and differential solutions. The blue circle indicates the ground-truth surveyed postion of P222. The green \"plus\" marker shows the differentially-corrected solution. The orange \"minus\" marker indicates the uncorrected position. Zoom out further on the map to make the blue ZOA1 receiver position visible.\n",
+    "\n",
+    "Let's see if it's closer to ground truth?"
    ]
   },
   {
@@ -327,7 +352,7 @@
    "outputs": [],
    "source": [
     "print(f\"ZOA1 position adjustment: {1e3*np.linalg.norm(result.atVector(zoa1_pos_key) - zoa1_ref_pos)} millimeters\")\n",
-    "error = ground_truth - receiver_position\n",
+    "error = ground_truth - corrected_receiver_position\n",
     "print(f\"Total positioning error: {np.linalg.norm(error)} meters\")"
    ]
   },
@@ -347,7 +372,7 @@
     "- [PseudorangeFactor.h](https://github.com/borglab/gtsam/blob/develop/gtsam/navigation/PseudorangeFactor.h)\n",
     "- [PseudorangeFactor.cpp](https://github.com/borglab/gtsam/blob/develop/gtsam/navigation/PseudorangeFactor.cpp)\n",
     "- [PseudorangeFactor.ipynb](https://github.com/borglab/gtsam/blob/develop/gtsam/navigation/doc/PseudorangeFactor.ipynb)\n",
-    "- [SinglePointPositioningExample.ipynb](https://github.com/masoug/gtsam-gnss/blob/develop/python/gtsam/examples/SinglePointPositioningExample.ipynb)"
+    "- [SinglePointPositioningExample.ipynb](https://borglab.github.io/gtsam/singlepointpositioningexample/)"
    ]
   }
  ],