Feature/coord fmt (#46)

* update .toml * update .toml * update README.md * update README.md * update python version of devcontainer to 3.11 * refactor * add tests for GNSSPositionEstimator * fix crash if navigation file does not contain a sys in rin obs * update .toml * add output crs as user argument * update requirements.txt * add example * fix typo * add example * update README.md
paarnes · Jan 12, 2025 · 26f8027 · 26f8027
1 parent 88aaf9f
commit 26f8027
Show file tree

Hide file tree

Showing 14 changed files with 312 additions and 199 deletions.
diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile
@@ -1,24 +1,16 @@
 # [Choice] Python version (use -bullseye variants on local arm64/Apple Silicon): 3, 3.10, 3.9, 3.8, 3.7, 3.6, 3-bullseye, 3.10-bullseye, 3.9-bullseye, 3.8-bullseye, 3.7-bullseye, 3.6-bullseye, 3-buster, 3.10-buster, 3.9-buster, 3.8-buster, 3.7-buster, 3.6-buster
-ARG VARIANT="3.10-bullseye"
-FROM mcr.microsoft.com/vscode/devcontainers/python:0-${VARIANT}
+ARG VARIANT="3.11-bullseye"
+FROM mcr.microsoft.com/vscode/devcontainers/python:${VARIANT}
 
-# [Choice] Node.js version: none, lts/*, 16, 14, 12, 10
-ARG NODE_VERSION="none"
-RUN if [ "${NODE_VERSION}" != "none" ]; then su vscode -c "umask 0002 && . /usr/local/share/nvm/nvm.sh && nvm install ${NODE_VERSION} 2>&1"; fi
-
-# [Optional] If your pip requirements rarely change, uncomment this section to add them to the image.
- COPY .devcontainer/requirements.txt /tmp/pip-tmp/
-  RUN pip3 --disable-pip-version-check --no-cache-dir install -r /tmp/pip-tmp/requirements.txt \
-    && rm -rf /tmp/pip-tmp
+COPY .devcontainer/requirements.txt /tmp/pip-tmp/
+RUN pip3 --disable-pip-version-check --no-cache-dir install -r /tmp/pip-tmp/requirements.txt \
+&& rm -rf /tmp/pip-tmp
 
 RUN python -m pip install --upgrade pip
 RUN python -m pip install notebook jupyterlab
 
-# [Optional] Uncomment this section to install additional OS packages.
 RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \
      && apt-get -y install --no-install-recommends libpq-dev gdal-bin libgdal-dev
-     # <your-package-list-here>
 
-# [Optional] Uncomment this line to install global node packages.
-# RUN su vscode -c "source /usr/local/share/nvm/nvm.sh && npm install -g <your-package-here>" 2>&1
+
 
diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json
@@ -9,7 +9,7 @@
 			// Update 'VARIANT' to pick a Python version: 3, 3.10, 3.9, 3.8, 3.7, 3.6
 			// Append -bullseye or -buster to pin to an OS version.
 			// Use -bullseye variants on local on arm64/Apple Silicon.
-			"VARIANT": "3.10-bullseye",
+			"VARIANT": "3.11-bullseye",
 			// Options
 			"NODE_VERSION": "lts/*"
 		}
@@ -38,7 +38,6 @@
 				"ms-python.vscode-pylance",
 				"42Crunch.vscode-openapi",
 				"ms-python.black-formatter",
-				"waderyan.gitblame",
 				"ms-toolsai.jupyter",
 				"ms-toolsai.jupyter-keymap",
 				"ms-toolsai.jupyter-renderers"

diff --git a/.devcontainer/requirements.txt b/.devcontainer/requirements.txt
@@ -3,4 +3,5 @@ numpy==1.24.4
 pandas==2.0.3
 tqdm==4.66.5
 zstandard==0.23.0
-pytest
+pytest
+pyproj
diff --git a/README.md b/README.md
@@ -3,6 +3,11 @@
 [![Python application](https://github.com/paarnes/GNSS_Multipath_Analysis_Software/actions/workflows/run-tests.yml/badge.svg)](https://github.com/paarnes/GNSS_Multipath_Analysis_Software/actions/workflows/run-tests.yml)
 [![PyPI version](https://badge.fury.io/py/gnssmultipath.svg)](https://badge.fury.io/py/gnssmultipath)
 [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
+[![Python Versions](https://img.shields.io/pypi/pyversions/gnssmultipath.svg)](https://pypi.org/project/gnssmultipath/)
+[![Downloads](https://pepy.tech/badge/gnssmultipath)](https://pepy.tech/project/gnssmultipath)
+
+
+
 
 ## Table of Contents
 - [Introduction](#introduction)
@@ -59,6 +64,8 @@ GNSS Multipath Analysis is a software tool for analyzing the multipath effect on
 - Allows selection of specific navigation systems and signal bands for analysis.
 - Estimate the approximate position of the receiver using pseudoranges from the RINEX observation file.
   - Supports both SP3 and RINEX navigation files.
+  - The software will estimate the receiver's position if it is not provided in the header of the RINEX observation file.
+  - Supports user-defined Coordinate Reference System (CRS). The estimated coordinates can be delivered in the desired CRS.
   - Calculates statistical measures for the estimated position, including:
     - Residuals
     - Sum of Squared Errors (SSE)
@@ -110,7 +117,7 @@ If you have a SP3 file, and not a RINEX navigation file, you just replace the ke
 1. Reads in the RINEX observation file
 2. Reads the RINEX navigation file or the precise satellite coordinates in SP3-format (depends on what’s provided)
 3. If a navigation file is provided, the satellite coordinates will be transformed from Kepler-elements to ECEF for GPS, Galileo and BeiDou. For GLONASS the navigation file is containing a state vector. The coordinates then get interpolated to the current epoch by solving the differential equation using a 4th order Runge-Kutta. If a SP3 file is provided, the interpolation is done using ``Neville's algorithm``.
-4. Satellites elevation and azimuth angles get computed.
+4. Compute the satellites' elevation and azimuth angles. If the receiver's approximate position is not provided in the header of the RINEX observation file, the software automatically estimates it based on pseudoranges using the ``GNSSPositionEstimator`` class.
 5. Cycle slip detection by using both ionospheric residuals and a code-phase combination. These linear combinations are given as
 
 $$
@@ -443,6 +450,34 @@ print(f'\nDOP values:\n' + '\n'.join([f'{k} = {v}' for k, v in stats["DOPs"].ite
 ```
 
 
+### Estimate the receiver's position based on pseudoranges in the desired CRS
+Define a specific Coordinate Reference System (CRS) to output the estimated receiver's coordinates. In this case the
+coordinates will be given in WGS84 UTM zone 32N (EPSG:32632) and ellipsoidal heights.
+
+Note: You can use the [EPSG GeoRepository](https://epsg.org/home.html) to find the EPSG code for the desired CRS.
+
+```python
+from gnssmultipath import GNSSPositionEstimator
+import numpy as np
+
+
+rinObs = 'OPEC00NOR_S_20220010000_01D_30S_MO_3.04_croped.rnx'
+rinNav = 'BRDC00IGS_R_20220010000_01D_MN.rnx'
+
+# Set desired time for when to estimate position and which system to use
+desired_time = np.array([2022, 1, 1, 1, 5, 30.0000000])
+desired_system = "E"  # GPS
+desired_crs = "EPSG:32632"  # Desired CRS for the estimated receiver coordinates (WGS84 UTM zone 32N)
+gnsspos, stats = GNSSPositionEstimator(rinObs,
+                                    rinex_nav_file=rinNav,
+                                    desired_time = desired_time,
+                                    desired_system = desired_system,
+                                    elevation_cut_off_angle = 10,
+                                    crs=desired_crs).estimate_position()
+
+print('Estimated coordinates in ECEF (m):\n' + '\n'.join([f'{axis} = {coord}' for axis, coord in zip(['Eastin', 'Northing', 'Height (ellipoidal)'], np.round(gnsspos[:-1], 3))]))
+```
+
 
 ## Some background information on implementation
 
@@ -532,7 +567,7 @@ $$
 \cos(\nu) = \frac{\cos(E) - e}{1 - e \cos(E)}, \quad \sin(\nu) = \frac{\sqrt{1 - e^2} \sin(E)}{1 - e \cos(E)}
 $$
 
-Use the arctangent to find $ \nu $:
+Use the arctangent to find $\nu$:
 
 $$
 \nu = \arctan2(\sin(\nu), \cos(\nu))
@@ -617,7 +652,7 @@ The Earth's rotation during the signal's travel introduces a positional error if
 
 **Iterative Process**:
 
-Update the longitude of the ascending node ($ \Omega_k $) to account for the Earth's rotation during the signal travel time:
+Update the longitude of the ascending node ($\Omega_k$) to account for the Earth's rotation during the signal travel time:
 
 $$
 \begin{equation*}
@@ -834,7 +869,7 @@ $$
 t_{n+1} = t_n + t_\text{step}
 $$
 
-**Reduce $ \Delta t $**:
+**Reduce**$\Delta t$:
 Adjust the remaining time difference:
 
 $$

diff --git a/pyproject.toml b/pyproject.toml
@@ -5,7 +5,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "gnssmultipath"
-version = "1.4.4.6"
+version = "1.5.0"
 authors = [
   { name="Per Helge Aarnes", email="[email protected]" },
 ]
@@ -15,6 +15,20 @@ maintainers  = [
 description = "A python software that provides comprehensive solutions for GNSS multipath analysis."
 readme = "README.md"
 requires-python = ">=3.8"
+keywords = [
+    "GNSS",
+    "multipath",
+    "analytics",
+    "gps",
+    "glonass",
+    "galileo",
+    "beidou",
+    "rinex",
+    "sp3",
+    "snr",
+    "position-estimation",
+]
+license = { text = "MIT" }
 classifiers = [
     "Programming Language :: Python :: 3",
     "Programming Language :: Python :: 3.8",

diff --git a/src/Examples_on_how_to_run_it.ipynb b/src/Examples_on_how_to_run_it.ipynb
@@ -12,13 +12,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
     "import os\n",
     "from gnssmultipath import GNSS_MultipathAnalysis\n",
     "\n",
+    "\n",
     "base_path = os.getcwd()\n",
     "parent_dir = os.path.abspath(os.path.join(base_path, os.pardir))\n",
     "\n",
@@ -230,49 +231,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 41,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "INFO(rinexReadObsFileHeader304): Rinex header has been read\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "Rinex observations are being read: 100%|██████████| (100/100)\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "INFO(readRinexObs304): The following GNSS systems have been read into the data: GPS, GLONASS, Galileo, Beidou\n",
-      "INFO(readRinexObs304): The following GPS observation types have been registered: C1C, L1C, C1P, C2W, L2W, C2X, L2X, C5X, L5X\n",
-      "INFO(readRinexObs304): The following GLONASS observation types have been registered: C1C, L1C, C1P, L1P, C2P, L2P, C2C, L2C\n",
-      "INFO(readRinexObs304): The following Galileo observation types have been registered: C1X, L1X, C7X, L7X, C5X, L5X, C8X, L8X\n",
-      "INFO(readRinexObs304): The following Beidou observation types have been registered: C2X, L2X, C7X, L7X, C6X, L6X\n",
-      "INFO(readRinexObs304): LLI have been read (if present in observation file)\n",
-      "INFO(readRinexObs304): SS have been read (if present in observation file)\n",
-      "INFO(readRinexObs304): Total processing time: 0.359375 seconds\n",
-      "\n",
-      "\n",
-      "Estimated coordinates in ECEF (m):\n",
-      "X = 3149786.364\n",
-      "Y = 598260.72\n",
-      "Z = 5495355.141\n",
-      "\n",
-      "Standard deviation of the estimated coordinates (m):\n",
-      "Sx = 4.233\n",
-      "Sy = 2.409\n",
-      "Sz = 7.096\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from gnssmultipath import GNSSPositionEstimator\n",
     "import numpy as np\n",
@@ -303,61 +264,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 57,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "Rinex navigation file is being read: 100%|██████████| (100/100)\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "INFO(rinexReadObsFileHeader304): Rinex header has been read\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "Rinex observations are being read: 100%|██████████| (100/100)\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "INFO(readRinexObs304): The following GNSS systems have been read into the data: GPS, GLONASS, Galileo, Beidou\n",
-      "INFO(readRinexObs304): The following GPS observation types have been registered: C1C, L1C, C1P, C2W, L2W, C2X, L2X, C5X, L5X\n",
-      "INFO(readRinexObs304): The following GLONASS observation types have been registered: C1C, L1C, C1P, L1P, C2P, L2P, C2C, L2C\n",
-      "INFO(readRinexObs304): The following Galileo observation types have been registered: C1X, L1X, C7X, L7X, C5X, L5X, C8X, L8X\n",
-      "INFO(readRinexObs304): The following Beidou observation types have been registered: C2X, L2X, C7X, L7X, C6X, L6X\n",
-      "INFO(readRinexObs304): LLI have been read (if present in observation file)\n",
-      "INFO(readRinexObs304): SS have been read (if present in observation file)\n",
-      "INFO(readRinexObs304): Total processing time: 0.359375 seconds\n",
-      "\n",
-      "\n",
-      "Estimated coordinates in ECEF (m):\n",
-      "X = 3149785.771\n",
-      "Y = 598281.355\n",
-      "Z = 5495359.165\n",
-      "\n",
-      "Standard deviation of the estimated coordinates (m):\n",
-      "Sx = 1.483\n",
-      "Sy = 1.634\n",
-      "Sz = 2.343\n",
-      "\n",
-      "DOP values:\n",
-      "PDOP = 1.66\n",
-      "TDOP = 0.833\n",
-      "GDOP = 1.857\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from gnssmultipath import GNSSPositionEstimator\n",
     "import numpy as np\n",
@@ -380,6 +289,40 @@
     "print('\\nStandard deviation of the estimated coordinates (m):\\n' + '\\n'.join([f'{k} = {v}' for k, v in stats[\"Standard Deviations\"].items() if k in ['Sx', 'Sy', 'Sz']]))\n",
     "print(f'\\nDOP values:\\n' + '\\n'.join([f'{k} = {v}' for k, v in stats[\"DOPs\"].items()]))"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Get approximate recevier coordinates in user defined CRS"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from gnssmultipath import GNSSPositionEstimator\n",
+    "import numpy as np\n",
+    "\n",
+    "\n",
+    "rinObs = os.path.join(path_to_testdata, 'ObservationFiles/OPEC00NOR_S_20220010000_01D_30S_MO_3.04_croped.rnx')\n",
+    "rinNav = os.path.join(path_to_testdata, 'NavigationFiles/BRDC00IGS_R_20220010000_01D_MN.rnx')\n",
+    "\n",
+    "# Set desired time for when to estimate position and which system to use\n",
+    "desired_time = np.array([2022, 1, 1, 1, 5, 30.0000000])\n",
+    "desired_system = \"E\"  # GPS\n",
+    "desired_crs = \"EPSG:32632\"  # Desired CRS for the estimated receiver coordinates (WGS84 UTM zone 32N)\n",
+    "gnsspos, stats = GNSSPositionEstimator(rinObs,\n",
+    "                                    rinex_nav_file=rinNav,\n",
+    "                                    desired_time = desired_time,\n",
+    "                                    desired_system = desired_system,\n",
+    "                                    elevation_cut_off_angle = 10,\n",
+    "                                    crs=desired_crs).estimate_position()\n",
+    "\n",
+    "print('Estimated coordinates in ECEF (m):\\n' + '\\n'.join([f'{axis} = {coord}' for axis, coord in zip(['Eastin', 'Northing', 'Height (ellipoidal)'], np.round(gnsspos[:-1], 3))]))"
+   ]
   }
  ],
  "metadata": {