Feat/mapping functions

AGENTS.md

@@ -0,0 +1,62 @@
+# AGENTS.md
+
+## Project Overview
+
+This project, `PyGNSS-TEC`, is a Python package designed for processing and analyzing Total Electron Content (TEC) data from Global Navigation Satellite System (GNSS) observations. It is a hybrid Python and Rust project, with the performance-intensive RINEX file parsing implemented in Rust and exposed to Python using `pyo3` and `maturin`.
+
+The core functionalities of the library include:
+- Reading RINEX observation and navigation files (versions 2.x and 3.x), including Hatanaka compressed files.
+- Calculating Slant and Vertical Total Electron Content (STEC and VTEC).
+- Correcting for satellite and receiver differential code biases (DCBs).
+- Supporting multiple GNSS constellations.
+
+The library uses the [Polars](https://pola.rs/) DataFrame library for efficient data manipulation.
+
+## Building and Running
+
+### Dependencies
+
+- Python >= 3.10
+- Rust
+- `uv` (recommended for Python environment management)
+
+### Building the project
+
+1.  Clone the repository:
+    ```bash
+    git clone https://github.com/Eureka-0/pygnss-tec.git
+    cd pygnss-tec
+    ```
+
+2.  Build the project using `maturin`:
+    ```bash
+    uv run maturin build --release
+    ```
+    The compiled wheel will be in the `target/wheels` directory.
+
+### Running tests
+
+The project uses `pytest` for testing. To run the tests, execute the following command from the root of the project:
+
+```bash
+uv run pytest
+```
+
+## Development Conventions
+
+- **Linting**: The project uses `ruff` for Python code linting.
+- **Code Formatting**: The project uses `ruff` to format the Python code.
+
+You can run the linter and formatter with the following commands:
+
+```bash
+# Lint
+uv run ruff check .
+
+# Format
+uv run ruff format .
+```
+
+- **Rust**: The Rust code follows standard Rust conventions and is formatted using `cargo fmt`.
+- **Testing**: Tests are located in the `tests` directory and use `pytest`. The tests cover both the Python and the Rust parts of the library.
+- **Commits**: There are no explicit commit message conventions mentioned, but the commit history shows a preference for descriptive and concise messages.

gnss_tec/rinex/__init__.py

@@ -1,8 +1,15 @@
 from .read_rinex import (
     ALL_CONSTELLATIONS,
+    LEAP_SECONDS,
     RinexObsHeader,
     get_leap_seconds,
     read_rinex_obs,
 )
 
-__all__ = ["RinexObsHeader", "ALL_CONSTELLATIONS", "get_leap_seconds", "read_rinex_obs"]
+__all__ = [
+    "ALL_CONSTELLATIONS",
+    "LEAP_SECONDS",
+    "RinexObsHeader",
+    "get_leap_seconds",
+    "read_rinex_obs",
+]

gnss_tec/rinex/read_rinex.py

@@ -21,6 +21,47 @@
 }
 """All supported GNSS constellations for RINEX file reading."""
 
+LEAP_SECONDS = [
+    (pl.datetime(1980, 1, 1), pl.duration(seconds=0)),
+    (pl.datetime(1981, 7, 1), pl.duration(seconds=1)),
+    (pl.datetime(1982, 7, 1), pl.duration(seconds=2)),
+    (pl.datetime(1983, 7, 1), pl.duration(seconds=3)),
+    (pl.datetime(1985, 7, 1), pl.duration(seconds=4)),
+    (pl.datetime(1988, 1, 1), pl.duration(seconds=5)),
+    (pl.datetime(1990, 1, 1), pl.duration(seconds=6)),
+    (pl.datetime(1991, 1, 1), pl.duration(seconds=7)),
+    (pl.datetime(1992, 7, 1), pl.duration(seconds=8)),
+    (pl.datetime(1993, 7, 1), pl.duration(seconds=9)),
+    (pl.datetime(1994, 7, 1), pl.duration(seconds=10)),
+    (pl.datetime(1996, 1, 1), pl.duration(seconds=11)),
+    (pl.datetime(1997, 7, 1), pl.duration(seconds=12)),
+    (pl.datetime(1999, 1, 1), pl.duration(seconds=13)),
+    (pl.datetime(2006, 1, 1), pl.duration(seconds=14)),
+    (pl.datetime(2009, 1, 1), pl.duration(seconds=15)),
+    (pl.datetime(2012, 7, 1), pl.duration(seconds=16)),
+    (pl.datetime(2015, 7, 1), pl.duration(seconds=17)),
+    (pl.datetime(2016, 12, 31), pl.duration(seconds=18)),
+]
+"""List of leap seconds as (datetime, duration) tuples."""
+
+
+def get_leap_seconds(time_col: pl.Expr | str) -> pl.Expr:
+    if isinstance(time_col, str):
+        time_col = pl.col(time_col)
+    time_col = time_col.dt.replace_time_zone(None)
+
+    expr = None
+    for dt, duration in LEAP_SECONDS[::-1]:
+        if expr is None:
+            expr = pl.when(time_col.ge(dt)).then(duration)
+        else:
+            expr = expr.when(time_col.ge(dt)).then(duration)
+
+    if expr is None:
+        raise ValueError("LEAP_SECONDS list is empty.")
+
+    return expr
+
 
 @dataclass
 class RinexObsHeader:
@@ -270,122 +311,3 @@ def calc_az_el(df: pl.DataFrame) -> pl.DataFrame:
         )
 
     return header, lf
-
-
-def get_leap_seconds(time_col: pl.Expr | str) -> pl.Expr:
-    if isinstance(time_col, str):
-        time_col = pl.col(time_col)
-
-    return (
-        pl.when(
-            time_col.is_between(
-                pl.datetime(1980, 1, 1), pl.datetime(1981, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=0))
-        .when(
-            time_col.is_between(
-                pl.datetime(1981, 7, 1), pl.datetime(1982, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=1))
-        .when(
-            time_col.is_between(
-                pl.datetime(1982, 7, 1), pl.datetime(1983, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=2))
-        .when(
-            time_col.is_between(
-                pl.datetime(1983, 7, 1), pl.datetime(1985, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=3))
-        .when(
-            time_col.is_between(
-                pl.datetime(1985, 7, 1), pl.datetime(1988, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=4))
-        .when(
-            time_col.is_between(
-                pl.datetime(1988, 1, 1), pl.datetime(1990, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=5))
-        .when(
-            time_col.is_between(
-                pl.datetime(1990, 1, 1), pl.datetime(1991, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=6))
-        .when(
-            time_col.is_between(
-                pl.datetime(1991, 1, 1), pl.datetime(1992, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=7))
-        .when(
-            time_col.is_between(
-                pl.datetime(1992, 7, 1), pl.datetime(1993, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=8))
-        .when(
-            time_col.is_between(
-                pl.datetime(1993, 7, 1), pl.datetime(1994, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=9))
-        .when(
-            time_col.is_between(
-                pl.datetime(1994, 7, 1), pl.datetime(1996, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=10))
-        .when(
-            time_col.is_between(
-                pl.datetime(1996, 1, 1), pl.datetime(1997, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=11))
-        .when(
-            time_col.is_between(
-                pl.datetime(1997, 7, 1), pl.datetime(1999, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=12))
-        .when(
-            time_col.is_between(
-                pl.datetime(1999, 1, 1), pl.datetime(2006, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=13))
-        .when(
-            time_col.is_between(
-                pl.datetime(2006, 1, 1), pl.datetime(2009, 1, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=14))
-        .when(
-            time_col.is_between(
-                pl.datetime(2009, 1, 1), pl.datetime(2012, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=15))
-        .when(
-            time_col.is_between(
-                pl.datetime(2012, 7, 1), pl.datetime(2015, 7, 1), "left"
-            )
-        )
-        .then(pl.duration(seconds=16))
-        .when(
-            time_col.is_between(
-                pl.datetime(2015, 7, 1), pl.datetime(2016, 12, 31), "left"
-            )
-        )
-        .then(pl.duration(seconds=17))
-        .when(time_col.ge(pl.datetime(2016, 12, 31)))
-        .then(pl.duration(seconds=18))
-        .otherwise(None)
-    )

gnss_tec/tec/constants.py

@@ -18,7 +18,7 @@
 SUPPORTED_CONSTELLATIONS = {"C": "BDS", "G": "GPS"}
 """Supported GNSS constellations for TEC calculation."""
 
-SIGNAL_FREQ = {
+SIGNAL_FREQ: dict[str, dict[str, float]] = {
     "C": {
         "B1-2": 1561.098e6,
         "B1": 1575.42e6,
@@ -31,14 +31,20 @@
 }
 """Signal frequencies for supported constellations and signals in Hz."""
 
-DEFAULT_C1_CODES = {
-    "C": ["C2I", "C2D", "C2X", "C1I", "C1D", "C1X", "C2W", "C1C"],
-    "G": ["C1W", "C1C", "C1X"],
+DEFAULT_C1_CODES: dict[str, dict[str, list[str]]] = {
+    "2": {"G": ["C1"]},
+    "3": {
+        "C": ["C2I", "C2D", "C2X", "C1I", "C1D", "C1X", "C2W", "C1C"],
+        "G": ["C1W", "C1C", "C1X"],
+    },
 }
 
-DEFAULT_C2_CODES = {
-    "C": ["C6I", "C6D", "C6X", "C7I", "C7D", "C7X", "C5I", "C5D", "C5X"],
-    "G": ["C2W", "C2C", "C2X", "C5W", "C5C", "C5X"],
+DEFAULT_C2_CODES: dict[str, dict[str, list[str]]] = {
+    "2": {"G": ["C2", "C5"]},
+    "3": {
+        "C": ["C6I", "C6D", "C6X", "C7I", "C7D", "C7X", "C5I", "C5D", "C5X"],
+        "G": ["C2W", "C2C", "C2X", "C5W", "C5C", "C5X"],
+    },
 }
 
 
@@ -58,10 +64,10 @@ class TECConfig:
     min_snr: float = 30.0
     """Minimum signal-to-noise ratio in dB-Hz."""
 
-    c1_codes: Mapping[str, list[str]] = field(default_factory=lambda: {})
+    c1_codes: Mapping[str, Mapping[str, list[str]]] = field(default_factory=lambda: {})
     """Observation codes priority list for C1 measurements."""
 
-    c2_codes: Mapping[str, list[str]] = field(default_factory=lambda: {})
+    c2_codes: Mapping[str, Mapping[str, list[str]]] = field(default_factory=lambda: {})
     """Observation codes priority list for C2 measurements."""
 
     rx_bias: Literal["external", "mstd", "lsq"] | None = "external"
@@ -73,6 +79,12 @@ class TECConfig:
         - None: Do not correct receiver bias.
     """
 
+    mapping_function: Literal["slm", "mslm"] = "slm"
+    """Mapping function to use:
+        - "slm": Single Layer Model
+        - "mslm": Modified Single Layer Model
+    """
+
     retain_intermediate: str | Iterable[str] | None | Literal["all"] = None
     """Names of intermediate columns to retain in the output DataFrame."""
 
@@ -82,38 +94,49 @@ def ipp_height_m(self) -> float:
         return self.ipp_height * 1e3
 
     @property
-    def code2band(self) -> Mapping[str, int]:
+    def mslm_height_m(self) -> float:
+        """Ionospheric pierce point height for Modified Single Layer Model in meters."""
+        return 506.7e3
+
+    @property
+    def alpha(self) -> float:
+        """Correction factor for Modified Single Layer Model."""
+        return 0.9782
+
+    def iter_c1_codes(
+        self, version: Literal["2", "3"]
+    ) -> Iterator[tuple[str, str, int]]:
+        for constellation, codes in self.c1_codes.get(version, {}).items():
+            for i, code in enumerate(codes):
+                yield constellation, code, i
+
+    def iter_c2_codes(
+        self, version: Literal["2", "3"]
+    ) -> Iterator[tuple[str, str, int]]:
+        for constellation, codes in self.c2_codes.get(version, {}).items():
+            for i, code in enumerate(codes):
+                yield constellation, code, i
+
+    def code2band(self, version: Literal["2", "3"]) -> Mapping[str, int]:
         code_band: dict[str, int] = {}
-        for const, code, _ in self.iter_c1_codes():
+        for const, code, _ in self.iter_c1_codes(version):
             code_band[f"{const}_{code}"] = 1  # C1 band
-        for const, code, _ in self.iter_c2_codes():
+        for const, code, _ in self.iter_c2_codes(version):
             code_band[f"{const}_{code}"] = 2  # C2 band
         return code_band
 
-    @property
-    def c1_priority(self) -> Mapping[str, int]:
+    def c1_priority(self, version: Literal["2", "3"]) -> Mapping[str, int]:
         priority: dict[str, int] = {}
-        for const, code, i in self.iter_c1_codes():
+        for const, code, i in self.iter_c1_codes(version):
             priority[f"{const}_{code}"] = i
         return priority
 
-    @property
-    def c2_priority(self) -> Mapping[str, int]:
+    def c2_priority(self, version: Literal["2", "3"]) -> Mapping[str, int]:
         priority: dict[str, int] = {}
-        for const, code, i in self.iter_c2_codes():
+        for const, code, i in self.iter_c2_codes(version):
             priority[f"{const}_{code}"] = i
         return priority
 
-    def iter_c1_codes(self) -> Iterator[tuple[str, str, int]]:
-        for constellation, codes in self.c1_codes.items():
-            for i, code in enumerate(codes):
-                yield constellation, code, i
-
-    def iter_c2_codes(self) -> Iterator[tuple[str, str, int]]:
-        for constellation, codes in self.c2_codes.items():
-            for i, code in enumerate(codes):
-                yield constellation, code, i
-
     def __post_init__(self):
         # Validate constellations
         allowed = set(SUPPORTED_CONSTELLATIONS.keys())