diff --git a/README.md b/README.md index f691859..f93d5de 100644 --- a/README.md +++ b/README.md @@ -501,9 +501,9 @@ A command line utility and Python class `RinexConverter` to convert binary GNSS **NB: This ALPHA preview release is limited to the following experimental functionality:** -1. [RINEX version 3.05](https://files.igs.org/pub/data/format/rinex305.pdf). +1. RINEX versions 3.04 and 4.02. 1. Convert binary UBX RXM-RAW or RXM-RAWX (raw observation) data from u-blox receivers (e.g. ZED-F9P) to RINEX Observation file format. -1. Convert binary RXM-SFRBX (navigation subframe) data from u-blox receivers to RINEX Navigation file format. **Currently only GPS LNAV data is supported**, but the underlying `RinexConverterNavigation` class is readily extensible. +1. Convert binary RXM-SFRBX (navigation subframe) data from u-blox receivers to RINEX Navigation file format. **Currently only GPS LNAV and CNAV data is supported**, but the underlying `RinexConverterNavigation` class is readily extensible. 1. Convert RTCM3 Ephemerides messages (1019, 1020, 1041-1046) from any source (including NTRIP caster or RTK base station receiver) to RINEX Navigation file format. 1. Convert NMEA MWD (wind speed and direction) and XDR (temperature and pressure) sensor data to RINEX Meteorology file format. diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index c265864..5853194 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -1,5 +1,16 @@ # pygnssutils +### RELEASE 1.2.1 + +FIXES: + +1. Restore (nested) pynmeagps dependency in pyproject.toml. Fixes #136 +1. Fix blank RINEX meteorology output file error. Fixes #138 + +CHANGES: + +1. Updates to experimental RINEX conversion to add support for GPS CNAV and RINEX version 4.02. + ### RELEASE 1.2.0 CHANGES: diff --git a/examples/process_rxmsfrbx_frames.py b/examples/process_rxmsfrbx_frames.py index cf5a915..28eeefd 100644 --- a/examples/process_rxmsfrbx_frames.py +++ b/examples/process_rxmsfrbx_frames.py @@ -1,129 +1,71 @@ """ process_rxmsfrbx_frames.py -Demo script to illustrate acquisition and collation of -raw NAV subframes from a UBX RXM-SFRBX data log using -the pygnssutils.RawNav utility class. +Illustration of parsing raw GPS LNAV and CNAV subframes from +UBX RXM-SFRBX binary datalog using RawNav utility class and +subframe definitions derived from the relevant GPS ICD. -This collates the following subframes for each SV4 -into a single RawNav object. +https://www.gps.gov/sites/default/files/2025-07/IS-GPS-200N.pdf - - subframe 1: clock corrections, sv health, etc. - - subframes 2 & 3: ephemerides - - subframe 4 page 18: (optional) ionospheric and time corrections - -(NB: subframe data dictionaries only currently defined -for GPS LNAV, but others can be added) - -python3 process_rxmsfrbx_frames.py infile=pygpsdata-rxmsfrbx.log - -Created on 21 Apr 2026 +Created on 14 May 2026 :author: semuadmin (Steve Smith) :copyright: semuadmin © 2026 :license: BSD 3-Clause """ -# pylint: disable=no-member - -from sys import argv - -from pyubx2 import UBXMessage - -from pygnssutils import GPS, GNSSReader, RawNav, RINEXProcessingError -from pygnssutils.rinex_subframes_gps import ( - GPS_LNAV_SUBFRAME_1, - GPS_LNAV_SUBFRAME_2, - GPS_LNAV_SUBFRAME_3, - GPS_LNAV_SUBFRAME_4_P18, -) - -SFR1 = 1 -SFR2 = 2 -SFR3 = 4 -SFR4 = 8 -TARGET_SFR = SFR1 | SFR2 | SFR3 # | SFR4 - - -def main(**kwargs): - """ - Main routine. - """ - - infile = kwargs["infile"] - - tot = 0 - sfr = 0 - err = 0 - lnavs = {} - navs = {} - with open(infile, "rb") as stream: - gnr = GNSSReader(stream) - for _, parsed in gnr: - if not isinstance(parsed, UBXMessage): - continue - if parsed.identity != "RXM-SFRBX": # UBX RXM-SFRBX (raw subframes) only - continue - tot += 1 - - sfrdata = RawNav.process_rxm_sfrbx(parsed) - if sfrdata.get("gnss", "") != GPS or sfrdata.get("sigid", "") not in ( - "1C", - ): - continue - - gnss = sfrdata["gnss"] - svid = sfrdata["svid"] - sigid = sfrdata["sigid"] - subframeid = sfrdata["subframeid"] - svcode = sfrdata.get("svcode", 0) - subframe = sfrdata["subframe"] - - try: - navs[(gnss, svid)] = navs.get((gnss, svid), RawNav(gnss, svid, sigid)) - nav = navs[(gnss, svid)] - if subframeid == 1: # clock parameters, sv health, etc. - wn = subframe >> 230 & 0b1111111111 - toc = (subframe >> 66 & 0b1111111111111111) * 16 - tow = subframe >> 253 & 0b11111111111111111 - # print(f"{tow=}, {wn=}, {toc=}") - nav.parse(subframe, GPS_LNAV_SUBFRAME_1) - elif subframeid == 2: # ephemerides - nav.parse(subframe, GPS_LNAV_SUBFRAME_2) - elif subframeid == 3: # ephemerides - nav.parse(subframe, GPS_LNAV_SUBFRAME_3) - elif subframeid == 4: - if svcode == 56: # page 18, ionospheric corrections - nav.parse(subframe, GPS_LNAV_SUBFRAME_4_P18) - if nav.sfracq & 0b111 == TARGET_SFR: - frame = navs.pop((gnss, svid)) - print(f"{str(frame)}\n") - # if frame.svid == 19: - # print(f"{frame.gnss}{frame.svid=:02d}, {frame.tow=}, {frame.iodc=}") - lnavs[svid] = lnavs.get(svid, 0) + 1 - sfr += 1 - except RINEXProcessingError: - err += 1 - - if tot: - print(f"Total RXM-SFRBX records: {tot}, errors: {err} ({err*100/tot:.0f}%)") - print(f"Total GPS LNAV subrames processed: {sfr}") - print( - ( - f"Total GPS LNAV full frames acquired (sfracq={TARGET_SFR}): " - f"{sum(lnavs.values())}" - ) - ) - print( - ( - f"Breakdown of GPS LNAV frames by SV: {sorted(lnavs.items())} " - f"Sum: {sum(lnavs.values())}" - ) - ) - else: - print("No RXM-SFRBX records in file") - - -if __name__ == "__main__": - - main(**dict(arg.split("=") for arg in argv[1:])) +from pygnssutils import GNSSReader +import pygnssutils.rinex_subframes_gps as rsg + +from pygnssutils.rawnav import RawNav + +SFRMAP = { + 1: rsg.GPS_LNAV_SUBFRAME_1, + 2: rsg.GPS_LNAV_SUBFRAME_2, + 3: rsg.GPS_LNAV_SUBFRAME_3, + 4: rsg.GPS_LNAV_GENERIC, + 5: rsg.GPS_LNAV_GENERIC, + 10: rsg.GPS_CNAV_SUBFRAME_10, + 11: rsg.GPS_CNAV_SUBFRAME_11, + 12: rsg.GPS_CNAV_SUBFRAME_12, + 13: rsg.GPS_CNAV_SUBFRAME_13, + 14: rsg.GPS_CNAV_SUBFRAME_14, + 15: rsg.GPS_CNAV_SUBFRAME_15, + 30: rsg.GPS_CNAV_SUBFRAME_30, + 31: rsg.GPS_CNAV_SUBFRAME_31, + 32: rsg.GPS_CNAV_SUBFRAME_32, + 33: rsg.GPS_CNAV_SUBFRAME_33, + 34: rsg.GPS_CNAV_SUBFRAME_34, + 35: rsg.GPS_CNAV_SUBFRAME_35, + 36: rsg.GPS_CNAV_SUBFRAME_36, + 37: rsg.GPS_CNAV_SUBFRAME_37, + 40: rsg.GPS_CNAV_SUBFRAME_40, +} + +INFILE = "pygpsdata-rxmsfrbx.log" + +rxm = 0 +gps = 0 +subframes = {} +with open(INFILE, "rb") as stream: + gnr = GNSSReader(stream) + for raw, parsed in gnr: + if parsed.identity == "RXM-SFRBX": + rxm += 1 + if parsed.gnssId == 0: # GPS + gps += 1 + output = RawNav.process_rxm_sfrbx(parsed) + sv = f"{output["gnss"]}{output["svid"]:02d}_{output["sigid"]}" + subframes[sv] = subframes.get(sv, {}) + msg = RawNav(output["gnss"], output["svid"], output["sigid"]) + msg.parse(output["subframe"], SFRMAP[output["subframeid"]], rsg.GPS_SFRACQ_MAP) + print(msg) + subframes[sv][output["subframeid"]] = ( + subframes[sv].get(output["subframeid"], 0) + 1 + ) + + # print summary of subframes captured + print(f"\nTotal records processed - RXM-SFRBX: {rxm:,}, GPS: {gps:,}") + print("Summary of subframes processed by PRN/Signal:") + for key, val in sorted(subframes.items()): + print(f"{key} {sorted(val.items())}") diff --git a/pyproject.toml b/pyproject.toml index ea03b96..4de4137 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ dependencies = [ "pyserial>=3.5", "pyspartn>=1.0.8", "pyubx2>=1.3.0", + "pynmeagps>=1.1.4", "pysbf2>=1.0.4", "pyubxutils>=1.0.6", "pyqgc>=0.2.2", diff --git a/src/pygnssutils/_version.py b/src/pygnssutils/_version.py index 3fd3ad2..7db9738 100644 --- a/src/pygnssutils/_version.py +++ b/src/pygnssutils/_version.py @@ -8,4 +8,4 @@ :license: BSD 3-Clause """ -__version__ = "1.2.0" +__version__ = "1.2.1" diff --git a/src/pygnssutils/rawnav.py b/src/pygnssutils/rawnav.py index 58aa9e8..9670abe 100644 --- a/src/pygnssutils/rawnav.py +++ b/src/pygnssutils/rawnav.py @@ -8,6 +8,23 @@ etc.) of one or more raw GNSS NAV subframe messages, as a precursor to RINEX conversion. +Once a RawNav object is instantiated, the `parse` function can +be invoked repeatedly to collate data from separate sequential +subframes e.g. for GPS LNAV, subframe 1 contains clock corrections, +subframes 2 & 3 contain ephemerides and subframe 4 page 18 contains +ionospheric corrections. + +An `sfracq` bitfield signifies which subframe IDs have been acquired, +and hence whether or not the RawNav object contains sufficient +information to be converted to a NAV record. + +A boolean `sequence` argument determines whether subframes are +processed as a sequence e.g. for GNSS where MSB and LSB attributes +are in separate, sequential (but not necessarily contiguous) subframes. + +Static methods are provided to facilitate acquisition of NAV subframe data +from UBX RXM-SFRBX messages. + The objective is to handle any GNSS subframe format for which: - data is available as a raw, unpadded little-endian integer. @@ -18,13 +35,6 @@ signed int and IEEE 754 float) have been defined and implemented in `bits2val`. -Once a RawNav object is instantiated, the `parse` function can -be invoked repeatedly to collate data from separate sequential -subframes e.g. for GPS LNAV, subframe 1 contains clock corrections, -subframes 2 & 3 contain ephemerides and subframe 4 page 18 contains -ionospheric corrections. An `sfracq` bitfield signifies which -subframe IDs have been acquired. - Format of data definition dictionary:: dict[field_name, tuple[offset, length, encoding, scaling] @@ -35,13 +45,6 @@ MSB and LSB field names MUST be suffixed "_msb" and "_lsb" respectively - the `parse` function will automatically combine them. -A boolean `sequence` argument determines whether subframes are -processed as a sequence e.g. for GNSS where MSB and LSB attributes -are in separate, sequential (but not necessarily contiguous) subframes. - -Static methods are provided to facilitate acquisition of NAV subframe data -from UBX RXM-SFRBX messages (currently for GPS LNAV only). - Created on 20 Apr 2026 :author: semuadmin (Steve Smith) @@ -50,15 +53,17 @@ """ import struct -from datetime import datetime, timezone +from datetime import datetime +from logging import getLogger from types import NoneType from typing import Literal -from pynmeagps import utc2wnotow, wnotow2utc +from pynmeagps import EPOCH0_GPS from pyubx2 import UBXMessage from pygnssutils.exceptions import RINEXProcessingError from pygnssutils.rinex_globals import GPS, UBXRINEXGNSS, UBXRINEXOBSCODE +from pygnssutils.rinex_helpers import get_epoch LSB = "_lsb" """MSB field name suffix""" @@ -104,13 +109,13 @@ def __init__( :param dict kwargs: optional keyword arguments """ + self.logger = getLogger(__name__) self._gnss = gnss self._svid = svid self._sigid = sigid - self._epoch = datetime.now(timezone.utc) - wno, tow, _ = utc2wnotow(self._epoch, gnss) - self.wn = wno - self.toc = int(tow / 1000) + self._epoch = EPOCH0_GPS + self.wn = 0 + self.toc = 0 self._sfracq = 0 self._msb = {} self._firsttoc = 999999999 @@ -122,6 +127,7 @@ def parse( self, data: int, definition: dict[str, tuple[int, int, str, int]], + sfrmap: dict | NoneType = None, sequence: bool = False, ): """ @@ -130,6 +136,7 @@ def parse( :param int data: raw, unpadded input data :param dict[str, tuple[int, int, str, int]] def: subframe definition dictionary \ (from GNSS ICD) + :param dict | NoneType sfrmap: map of subframe ids to sfracq bitmask (None) :param bool sequence: process subframe as part of a contiguous sequence (False) :raises: RINEXProcessingError """ @@ -137,6 +144,8 @@ def parse( try: + if sfrmap is None: + sfrmap = {} # get exemplary preamble value if one is available valpre = definition.pop(VALPREAMBLE, 0) @@ -180,23 +189,18 @@ def parse( setattr(self, att, val) if att == SFR: # update subframe acquisition status - self._sfracq |= int(2 ** (val - 1)) - - # update epoch with last acquisition timestamp - # TODO is toc the correction value to use here? - # doesn't seem to reflect actual UTC time and - # doesn't appear to increment by seconds like - # the tow in HOW??? - self._epoch = wnotow2utc( + self._sfracq |= sfrmap.get(val, 0) + + # update epoch with latest acquisition timestamp + # TODO check this results in the correct epoch? + epoch, wn = get_epoch( wno=int(getattr(self, WN)), tow=int(getattr(self, TOC) * 1000), - ls=None, gnss=self._gnss, - autoroll=True, - modwno=True, ) - wno, tow, _ = utc2wnotow(utc=self._epoch, gnss=self._gnss, modwno=False) - self.wn = wno + if epoch > self._epoch: + self._epoch = epoch + self.wn = wn if not sequence: self._store_orphaned_msb() @@ -308,6 +312,17 @@ def svid(self) -> int: return self._svid + @property + def svcode(self) -> str: + """ + Getter for SV code (gnss & prn). + + :return: svcode e.g. "G14" + :rtype: str + """ + + return f"{self._gnss}{self._svid:02d}" + @property def sigid(self) -> str: """ @@ -352,6 +367,9 @@ def process_rxm_sfrbx(data: UBXMessage) -> dict[str, str | int | float | NoneTyp """ Reassemble individual subframe from UBX RXM-SFRBX dwrds. + CURRENTLY ONLY GPS LNAV & CNAV IMPLEMENTED BUT METHOD + READILY EXTENSIBLE. + :param UBXMessage data: parsed UBX RXM-SFRBX message :return: dict of subframe attributes :rtype: dict[str, str | int | float | NoneType] @@ -382,15 +400,12 @@ def process_rxm_sfrbx(data: UBXMessage) -> dict[str, str | int | float | NoneTyp tow = None subframeid = 0 - # for GPS LNAV, words are 30 bits each, padded with 2 bits at end - if gnss == GPS and sigid == "1C": + # for GPS LNAV, subframe = 10 * 30 bits, with each 32-bit dwrd padded with 2 bits at end + if gnss == GPS and sigid == "1C": # GPS LNAV for i in range(numw): wrd = getattr(data, f"dwrd_{i+1:02d}") & 0xFFFFFFFC >> 2 subframe += wrd << (30 * (numw - 1 - i)) - # tow in HOW is 17 MSB of Z-count; full tow is 19 bits. - # rolls over at 100,799 - # tow = (subframe >> 251) & 0x7FFFC # << 2 to make 19 bits - tow = (subframe >> 253) & 0b11111111111111111 # << 2 to make 19 bits + tow = (subframe >> 253) & 0b11111111111111111 subframeid = (subframe >> 248) & 0x7 output = { "gnss": gnss, @@ -402,7 +417,23 @@ def process_rxm_sfrbx(data: UBXMessage) -> dict[str, str | int | float | NoneTyp } if subframeid in (4, 5): output["dataid"] = subframe >> 234 & 0x3 - output["svcode"] = subframe >> 232 & 0x3F + output["pageid"] = subframe >> 232 & 0x3F + # for GPS CNAV, subframe = 3 * 100 bits, final 20 bits of 320 bit dwrd is padding + elif gnss == GPS and sigid in ("2L", "2S", "5I", "5Q"): # GPS CNAV + for i in range(numw): + wrd = getattr(data, f"dwrd_{i+1:02d}") + subframe += wrd << (32 * (numw - 1 - i)) + subframe = (subframe >> 20) & (2**300 - 1) + subframeid = (subframe >> 280) & 0b111111 + tow = (subframe >> 263) & 0b11111111111111111 + output = { + "gnss": gnss, + "svid": svid, + "sigid": sigid, + "subframeid": subframeid, + "tow": tow, + "subframe": subframe, + } # elif gnss == x and sigid == x: # TODO add other subframe processing algorithms here... # (would need 'tow equivalent' for GLONASS time system diff --git a/src/pygnssutils/rinex_conv.py b/src/pygnssutils/rinex_conv.py index e24c464..eed5bc0 100644 --- a/src/pygnssutils/rinex_conv.py +++ b/src/pygnssutils/rinex_conv.py @@ -43,6 +43,7 @@ MET, NAV, OBS, + RINEX4, RINEX_CANCELLED, RINEX_ERROR, RINEX_NORECS, @@ -51,8 +52,11 @@ ) from pygnssutils.rinex_helpers import ( format_comments, + format_doi, format_filename, + format_licenseofuse, format_runby, + format_stationinfo, format_version, ) @@ -154,6 +158,9 @@ def __init__( NAV: {MAX: EPOCHMIN, MIN: EPOCHMAX, CUR: EPOCHMIN, FRQ: 0}, MET: {MAX: EPOCHMIN, MIN: EPOCHMAX, CUR: EPOCHMIN, FRQ: 0}, } + self._doi = kwargs.get("doi", "") + self._license = kwargs.get("license", "") + self._station = kwargs.get("station", "") self.verbosity = int(verbosity) self.logtofile = logtofile @@ -266,7 +273,7 @@ def process_input_data( for raw, parsed in gnr: if stopevent is not None: if stopevent.is_set(): - raise RINEXProcessingError("Cancelled") + raise KeyboardInterrupt("Terminated by user") if raw is not None: self._tot += 1 self._progress = int(round(100 * self._tot / self._msgcount, 0)) @@ -294,10 +301,10 @@ def process_input_data( self.process_output_data(rinextypes) res = RINEX_OK - # except (TypeError, ValueError, AttributeError) as err: - # self.logger.error(err) - # res = RINEX_ERROR - except (RINEXProcessingError, KeyboardInterrupt): + except RINEXProcessingError as err: + self.logger.error(f"Processing error {err}") + res = RINEX_ERROR + except KeyboardInterrupt: self.logger.warning("Terminated by user") res = RINEX_CANCELLED @@ -332,11 +339,18 @@ def format_header_common(self, rinextype: Literal["O", "N", "M"]) -> str: :rtype: str """ - return ( + hdr = ( format_version(self._rinex_version, rinextype, self._gnssfilter) + format_runby() + format_comments(self.user_comments) ) + if self._rinex_version >= RINEX4: + hdr += ( + format_doi(self._doi) + + format_licenseofuse(self._license) + + format_stationinfo(self._station) + ) + return hdr def output(self, data: str | Any | NoneType, rinextype: Literal["O", "N", "M"]): """ diff --git a/src/pygnssutils/rinex_conv_cli.py b/src/pygnssutils/rinex_conv_cli.py index e02928a..83e8c13 100644 --- a/src/pygnssutils/rinex_conv_cli.py +++ b/src/pygnssutils/rinex_conv_cli.py @@ -33,6 +33,7 @@ RINEX_OK, RINEXTYPE, RINEXVER_DEFAULT, + RINEXVERSIONS, ) from pygnssutils.rinex_helpers import listify @@ -44,10 +45,7 @@ def main(): ap = ArgumentParser( epilog=EPILOG, - description=( - "NB: %(prog)s is currently an experimental facility with " - "limited functionality. NOT FOR PRODUCTION USE." - ), + description=("NB: %(prog)s is currently an experimental facility."), formatter_class=ArgumentDefaultsHelpFormatter, ) ap.add_argument( @@ -64,8 +62,9 @@ def main(): "-R", "--rinexver", required=False, - help="RINEX Version e.g. '3.05'", + help="RINEX Version", type=str, + choices=RINEXVERSIONS, default=RINEXVER_DEFAULT, ) ap.add_argument( @@ -93,7 +92,7 @@ def main(): ap.add_argument( "--obsfilter", required=False, - help="Comma-separated list of observation codes to process, or blank for all e.g '1C,2B'", + help="Comma-separated list of observation codes to process, or blank for all e.g '1C,2S'", type=str, default="", ) @@ -153,6 +152,27 @@ def main(): type=str, default="", ) + ap.add_argument( + "--doi", + required=False, + help="Digital Object Identifier (RINEX 4 only)", + type=str, + default="", + ) + ap.add_argument( + "--license", + required=False, + help="License of Use (RINEX 4 only)", + type=str, + default="", + ) + ap.add_argument( + "--station", + required=False, + help="Station Information (RINEX 4 only)", + type=str, + default="", + ) ap.add_argument( "--comments", required=False, @@ -191,6 +211,9 @@ def main(): protfilter=int( kwargs.pop("protfilter", NMEA_PROTOCOL | UBX_PROTOCOL | RTCM3_PROTOCOL) ), + doi=kwargs.pop("doi", ""), + license=kwargs.pop("license", ""), + station=kwargs.pop("station", ""), **kwargs, ) res = rc.process_input( diff --git a/src/pygnssutils/rinex_conv_met.py b/src/pygnssutils/rinex_conv_met.py index cc503c7..5b380e9 100644 --- a/src/pygnssutils/rinex_conv_met.py +++ b/src/pygnssutils/rinex_conv_met.py @@ -27,6 +27,7 @@ from pyrtcm import RTCMMessage from pyubx2 import UBXMessage +from pygnssutils.exceptions import RINEXProcessingError from pygnssutils.globals import VERBOSITY_MEDIUM from pygnssutils.rinex_globals import BDS, COLWIDTH, EPOCHMIN, MET from pygnssutils.rinex_helpers import ( @@ -205,7 +206,7 @@ def convert_nmea_mwd(self, data: NMEAMessage): try: # NB: NMEA MWD sentence has no timestamp, so epoch must be # obtained from another NMEA RMC message in the same data stream - epoch = self.__app.current_epoch + epoch = self.__app.get_current_epoch(MET) if epoch == EPOCHMIN: # epoch not yet established return winddir = data.dirM @@ -228,7 +229,8 @@ def convert_nmea_mwd(self, data: NMEAMessage): self._sensortypes[obscode].get("count", 0) + 1 ) except (AttributeError, TypeError) as err: - print(f"something went wrong {err}") + raise RINEXProcessingError(err) from err + # print(f"something went wrong {err}") def convert_nmea_xdr(self, data: NMEAMessage): """ @@ -248,7 +250,7 @@ def geta(att: str, i: int): try: # NB: NMEA XDR sentence has no timestamp, so epoch must be # obtained from another NMEA RMC message in the same data stream - epoch = self.__app.current_epoch + epoch = self.__app.get_current_epoch(MET) if epoch == EPOCHMIN: # epoch not yet established return self._metdata[epoch] = self._metdata.get(epoch, {}) @@ -282,7 +284,8 @@ def geta(att: str, i: int): break except (AttributeError, TypeError) as err: - print(f"something went wrong {err}") + raise RINEXProcessingError(err) from err + # print(f"something went wrong {err}") def get_nmea_epoch(self, data: NMEAMessage): """ diff --git a/src/pygnssutils/rinex_conv_nav.py b/src/pygnssutils/rinex_conv_nav.py index 7b962e0..8158577 100644 --- a/src/pygnssutils/rinex_conv_nav.py +++ b/src/pygnssutils/rinex_conv_nav.py @@ -27,11 +27,11 @@ from datetime import datetime, timezone from logging import getLogger -from math import pi +from math import pi, sqrt from types import NoneType from typing import Any, Literal -from pynmeagps import NMEAMessage, wnotow2utc +from pynmeagps import NMEAMessage from pyrtcm import RTCMMessage from pyubx2 import UBXMessage @@ -40,6 +40,7 @@ from pygnssutils.rawnav import RawNav from pygnssutils.rinex_globals import ( BDS, + EPHNAVTYPES, EPOCHMIN, GAL, GLO, @@ -47,31 +48,50 @@ IRN, NAV, QZS, + RINEX4, RINEX_URA, + SBA, ) from pygnssutils.rinex_helpers import ( DRNX, + format_eop, format_fileend, format_headerend, + format_ion, format_iono_corr, format_leapseconds, + format_nav_typesvmssg, + format_sto, format_time_corr, + format_timefirstlast, + get_epoch, get_fithours, get_svcode_rtcm, ) from pygnssutils.rinex_subframes_gps import ( + GPS_CNAV_SUBFRAME_10, + GPS_CNAV_SUBFRAME_11, + GPS_CNAV_SUBFRAME_30, + GPS_CNAV_SUBFRAME_32, + GPS_CNAV_SUBFRAME_33, GPS_LNAV_SUBFRAME_1, GPS_LNAV_SUBFRAME_2, GPS_LNAV_SUBFRAME_3, GPS_LNAV_SUBFRAME_4_P18, + GPS_SFRACQ_MAP, ) +AREF = 26559710 BOD = "bod" CLKBIAS = "clkbias" CLKDRIFT = "clkdrift" CLKRATE = "clkrate" +EPH = "EPH" EPOCH = "epoch" -TARGET_SFR = 0b111 +OMEGADOTREF = -2.6e-9 +RECTYPE = "rectype" +TARGET_SFR_GPS_LNAV = 0b111 # SFR 1(1), 2(2), 4p18(4) +TARGET_SFR_GPS_CNAV = 0b111 # SFR 10(1), 11(2), 30(4) class RinexConverterNavigation: @@ -120,6 +140,7 @@ def __init__( self._rinex_version = rinex_version self._ionocorr = {} self._timecorr = {} + self._eopcorr = {} self._leapseconds = "" self._navdata = {} # holder for converted RINEX nav data self._navframes = {} # holder for acquired partial NAV frames @@ -142,11 +163,8 @@ def process_input_data( res = 0 if isinstance(parsed, UBXMessage): if parsed.identity == "RXM-SFRBX": - self._convert_rxmsfrbx(parsed) + self._collate_rxmsfrbx(parsed) res += 1 - elif isinstance(parsed, RawNav) and self._datasource in ("R", "S", "U"): - self._convert_rawnav(parsed) - res += 1 elif isinstance(parsed, RTCMMessage) and self._datasource in ("N"): if parsed.identity in ("1005", "1006"): # station ID self._convert_rtcm1005(parsed) @@ -204,11 +222,12 @@ def _format_nav_data(self, navdata: dict[tuple, dict] | str = ""): """ Format navigation data for each svcode/iodc from navdata dict. - Format of navdata dict:: + Format of navdata dict: navdata = { (svcode (str), iodc (int)): { "epoch": epoch (datetime), + "rectype: rectype (str), "clkbias": clkbias (float), "clkdrift": clkdrift (float), "clkrate": clkrate (float), @@ -230,12 +249,25 @@ def _format_nav_data(self, navdata: dict[tuple, dict] | str = ""): if navdata == "": navdata = {} + if self._rinex_version >= RINEX4: + for _, sto in self._timecorr.items(): # RINEX 4 system time offsets + self.__app.output(sto, NAV) + for _, ion in self._ionocorr.items(): # RINEX 4 ionospheric corrections + self.__app.output(ion, NAV) + for _, eop in self._eopcorr.items(): # RINEX 4 earth orient corrections + self.__app.output(eop, NAV) + for (svcode, _), data in navdata.items(): # TODO sort??? + if self._rinex_version >= RINEX4: + rectyp = format_nav_typesvmssg(EPH, svcode, data[RECTYPE]) + else: + rectyp = "" epoch = data[EPOCH].strftime("%Y %m %d %H %M %S") clkbias = data[CLKBIAS] clkdrift = data[CLKDRIFT] clkrate = data[CLKRATE] nav = ( + f"{rectyp}" f"{svcode} {epoch}{DRNX(clkbias,19,12)}" f"{DRNX(clkdrift,19,12)}{DRNX(clkrate,19,12)}\n" ) @@ -251,16 +283,21 @@ def _format_header(self): Format navigation header lines. """ - hdr = ( - self.__app.format_header_common(NAV) - + format_iono_corr(self._ionocorr) - + format_time_corr(self._timecorr) - + format_leapseconds( - self.__app.get_start_epoch(NAV), # TODO check this is correct date - self._gnss_filter, - ) - + format_headerend() + hdr = self.__app.format_header_common(NAV) + if self._rinex_version < RINEX4: + for _, ion in self._ionocorr.items(): # RINEX 3 ionospheric corrections + hdr += ion + for _, sto in self._timecorr.items(): # RINEX 3 system time offsets + hdr += sto + hdr += format_leapseconds( + self.__app.get_start_epoch(NAV), # TODO check this is correct date + self._gnss_filter, ) + # debug only vvv + # hdr += format_timefirstlast(self.__app.get_start_epoch(NAV), "FIRST") + # hdr += format_timefirstlast(self.__app.get_end_epoch(NAV), "LAST") + # debug only ^^^ + hdr += format_headerend() self.__app.output(hdr, NAV) def _convert_rtcm1005(self, data: RTCMMessage): @@ -289,7 +326,8 @@ def _convert_rtcm1019(self, data: RTCMMessage): """ svcode = get_svcode_rtcm(GPS, data.DF009) - epoch = self._get_epoch(wno=data.DF076, toc=data.DF081, gnss=GPS) + epoch, _ = get_epoch(wno=data.DF076, tow=data.DF081, gnss=GPS) + self.__app.set_current_epoch(epoch, NAV) tom = data.DF093 if epoch == EPOCHMIN: return @@ -298,6 +336,7 @@ def _convert_rtcm1019(self, data: RTCMMessage): nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[GPS, "L1 C/A"] nvd[CLKBIAS] = data.DF084 # clock bias nvd[CLKDRIFT] = data.DF083 # clock drift nvd[CLKRATE] = data.DF082 # clock drift rate @@ -354,7 +393,8 @@ def _convert_rtcm1020(self, data: RTCMMessage): svcode = get_svcode_rtcm(GLO, data.DF038) # TODO which GLONASS fields represent wno and toc? - # epoch = self._get_epoch(wno=data.DF???, toc=data.DF???, gnss=GLO) + # epoch = get_epoch(wno=data.DF???, tow=data.DF???, gnss=GLO) + # self.__app.set_current_epoch(epoch, NAV) epoch = self._get_external_epoch(NAV) if epoch == EPOCHMIN: return @@ -362,6 +402,7 @@ def _convert_rtcm1020(self, data: RTCMMessage): self._navdata[(svcode, iodc)] = {} nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[GLO, "L1 C/A"] nvd[CLKBIAS] = ( data.DF124 ) # SV clock bias (sec) (-TauN) TODO should val be negated? @@ -395,7 +436,7 @@ def _convert_rtcm1020(self, data: RTCMMessage): def _convert_rtcm1045(self, data: RTCMMessage): """ - Format Galileo 1045 F/NAV, 1046 I/NAV broadcast orbit blocks. + Format Galileo 1045 F/NAV broadcast orbit blocks. https://github.com/semuconsulting/pyrtcm/blob/3c763bd0016698864ccccbb9988403d2b3a5ba7f/src/pyrtcm/rtcmtypes_get.py#L868 @@ -403,7 +444,8 @@ def _convert_rtcm1045(self, data: RTCMMessage): """ svcode = get_svcode_rtcm(GAL, data.DF252) - epoch = self._get_epoch(wno=data.DF289, toc=data.DF304, gnss=GAL) + epoch, _ = get_epoch(wno=data.DF289, tow=data.DF304, gnss=GAL) + self.__app.set_current_epoch(epoch, NAV) tom = data.DF304 if epoch == EPOCHMIN: return @@ -411,6 +453,7 @@ def _convert_rtcm1045(self, data: RTCMMessage): self._navdata[(svcode, iodc)] = {} nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[GAL, "E5a"] nvd[CLKBIAS] = data.DF296 # clock bias nvd[CLKDRIFT] = data.DF295 # clock drift nvd[CLKRATE] = data.DF294 # clock drift rate @@ -465,7 +508,8 @@ def _convert_rtcm1046(self, data: RTCMMessage): """ svcode = get_svcode_rtcm(GAL, data.DF252) - epoch = self._get_epoch(wno=data.DF289, toc=data.DF304, gnss=GAL) + epoch, _ = get_epoch(wno=data.DF289, tow=data.DF304, gnss=GAL) + self.__app.set_current_epoch(epoch, NAV) tom = data.DF304 if epoch == EPOCHMIN: return @@ -473,6 +517,7 @@ def _convert_rtcm1046(self, data: RTCMMessage): self._navdata[(svcode, iodc)] = {} nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[GAL, "E5b"] nvd[CLKBIAS] = data.DF296 # clock bias nvd[CLKDRIFT] = data.DF295 # clock drift nvd[CLKRATE] = data.DF294 # clock drift rate @@ -527,7 +572,8 @@ def _convert_rtcm1042(self, data: RTCMMessage): """ svcode = get_svcode_rtcm(BDS, data.DF488) - epoch = self._get_epoch(wno=data.DF489, toc=data.DF493, gnss=BDS) + epoch, _ = get_epoch(wno=data.DF489, tow=data.DF493, gnss=BDS) + self.__app.set_current_epoch(epoch, NAV) tom = data.DF505 if epoch == EPOCHMIN: return @@ -536,6 +582,7 @@ def _convert_rtcm1042(self, data: RTCMMessage): nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[BDS, "B1C"] nvd[CLKBIAS] = data.DF496 # clock bias nvd[CLKDRIFT] = data.DF495 # clock drift nvd[CLKRATE] = data.DF494 # clock drift rate @@ -589,7 +636,8 @@ def _convert_rtcm1044(self, data: RTCMMessage): """ svcode = get_svcode_rtcm(QZS, data.DF429) - epoch = self._get_epoch(wno=data.DF452, toc=data.DF442, gnss=QZS) + epoch, _ = get_epoch(wno=data.DF452, tow=data.DF442, gnss=QZS) + self.__app.set_current_epoch(epoch, NAV) tom = data.DF442 if epoch == EPOCHMIN: return @@ -598,6 +646,7 @@ def _convert_rtcm1044(self, data: RTCMMessage): nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[QZS, "L1 C/A"] nvd[CLKBIAS] = data.DF433 # clock bias nvd[CLKDRIFT] = data.DF432 # clock drift nvd[CLKRATE] = data.DF431 # clock drift rate @@ -641,48 +690,6 @@ def _convert_rtcm1044(self, data: RTCMMessage): nvb[6][2] = "" # - Spare nvb[6][3] = "" # - Spare - # def _convert_navsource_sbas(self, data: RTCMMessage): - # """ - # Format SBAS broadcast orbit blocks. - - # :param RTCMMessage data: parsed SBAS Ephemerides message - # """ - - # svcode = get_svcode_rtcm(SBA, data.DF009) - # if self._useextepoch: - # epoch = self._get_external_epoch(NAV) - # else: - # epoch = self._get_local_epoch(data.wno, data.tow, SBA) - # if epoch == EPOCHMIN: - # return - # iodc = data.iodc - # self._navdata[(svcode, iodc)] = {} - # nvd = self._navdata[(svcode, iodc)] - - # nvd[EPOCH] = epoch - # nvd[CLKBIAS] = data.a0 # clock bias - # nvd[CLKDRIFT] = data.a1 # clock drift - # nvd[CLKRATE] = data.a2 # clock drift rate - # nvd[BOD] = [] - # nvb = nvd[BOD] - # for _ in range(3): # broadcast orbit data blocks * 3 - # nvb.append(["", "", "", ""]) # 4X,4D19.12 - # # BROADCAST ORBIT - 1 - # nvb[0][0] = data.posX # - Satellite position X (km) - # nvb[0][1] = data.velX # - velocity X dot (km/sec) - # nvb[0][2] = data.accX # - X acceleration (km/sec2) - # nvb[0][3] = data.svhealth # - health (0=healthy, 1=unhealthy) (MSB of 3-bit Bn) - # # BROADCAST ORBIT - 2 - # nvb[1][0] = data.posY # - Satellite position Y (km) - # nvb[1][1] = data.velY # - velocity Y dot (km/sec) - # nvb[1][2] = data.accY # - Y acceleration (km/sec2) - # nvb[1][3] = data.ura # - Accuracy code (URA, meters) - # # BROADCAST ORBIT - 3 - # nvb[2][0] = data.posZ # - Satellite position Z (km) - # nvb[2][1] = data.velZ # - velocity Z dot (km/sec) - # nvb[2][2] = data.accZ # - Z acceleration (km/sec2) - # nvb[2][3] = data.iodn # - IODN (Issue of Data Navigation - def _convert_rtcm1041(self, data: RTCMMessage): """ Format NavIC broadcast orbit blocks. @@ -693,7 +700,8 @@ def _convert_rtcm1041(self, data: RTCMMessage): """ svcode = get_svcode_rtcm(IRN, data.DF516) - epoch = self._get_epoch(wno=data.DF517, toc=data.DF522, gnss=IRN) + epoch, _ = get_epoch(wno=data.DF517, tow=data.DF522, gnss=IRN) + self.__app.set_current_epoch(epoch, NAV) tom = data.DF537 if epoch == EPOCHMIN: return @@ -702,6 +710,7 @@ def _convert_rtcm1041(self, data: RTCMMessage): nvd = self._navdata[(svcode, iodc)] nvd[EPOCH] = epoch + nvd[RECTYPE] = EPHNAVTYPES[IRN, "L1"] nvd[CLKBIAS] = data.DF518 # clock bias nvd[CLKDRIFT] = data.DF519 # clock drift nvd[CLKRATE] = data.DF520 # clock drift rate @@ -745,22 +754,21 @@ def _convert_rtcm1041(self, data: RTCMMessage): nvb[6][2] = "" # - Spare nvb[6][3] = "" # - Spare - def _convert_rawnav(self, data: RawNav): + def _convert_rawnav_gps_lnav(self, data: RawNav): """ - Format RawNav broadcast orbit blocks. + Format RawNav GPS LNAV broadcast orbit blocks. :param RawNav data: RawNav object containing data \ collated from UBX RXM-SFRBX messages or other \ raw NAV subframe sources. """ - svcode = get_svcode_rtcm(data.gnss, data.svid) - epoch = self._get_epoch(wno=data.wn, toc=data.toc, gnss=data.gnss) - iodc = data.iodc - self._navdata[(svcode, iodc)] = {} - nvd = self._navdata[(svcode, iodc)] + self.__app.set_current_epoch(data.epoch, NAV) + self._navdata[(data.svcode, data.iodc)] = {} + nvd = self._navdata[(data.svcode, data.iodc)] - nvd[EPOCH] = epoch + nvd[EPOCH] = data.epoch + nvd[RECTYPE] = "LNAV" nvd[CLKBIAS] = data.af0 # clock bias nvd[CLKDRIFT] = data.af1 # clock drift nvd[CLKRATE] = data.af2 # clock drift rate @@ -805,134 +813,372 @@ def _convert_rawnav(self, data: RawNav): nvb[6][2] = "" # - Spare nvb[6][3] = "" # - Spare - def _convert_ionocorr(self, data: Any, gnss: str): + def _convert_rawnav_gps_cnav(self, data: RawNav): + """ + Format RawNav GPS CNAV broadcast orbit blocks. + + :param RawNav data: RawNav object containing data \ + collated from UBX RXM-SFRBX messages or other \ + raw NAV subframe sources. + """ + + self.__app.set_current_epoch(data.epoch, NAV) + iodc = data.toc # TODO confirm + self._navdata[(data.svcode, iodc)] = {} + nvd = self._navdata[(data.svcode, iodc)] + + nvd[EPOCH] = data.epoch + nvd[RECTYPE] = "CNAV" + nvd[CLKBIAS] = data.af0n # clock bias + nvd[CLKDRIFT] = data.af1n # clock drift + nvd[CLKRATE] = data.af2n # clock drift rate + nvd[BOD] = [] + nvb = nvd[BOD] + for _ in range(8): # broadcast orbit data blocks * 8 + nvb.append(["", "", "", ""]) # 4X,4D19.12 + # Multiply by pi to convert semicircles to radians + # BROADCAST ORBIT - 1 + nvb[0][0] = data.adot # - Issue of Data, Ephemeris SFR10 + nvb[0][1] = data.crs # - Crs (meters) SFR11 + nvb[0][2] = data.deltan0 * pi # - Delta n (radians/sec) SFR10 + nvb[0][3] = data.m0 * pi # - M0 (radians) SFR10 + # BROADCAST ORBIT - 2 + nvb[1][0] = data.cuc # - Cuc (radians) SFR11 + nvb[1][1] = data.e # - e Eccentricity SFR10 + nvb[1][2] = data.cus # - Cus (radians) SFR11 + # nvb[1][3] = data.sqrta # - sqrt(a) (sqrt(m)) SFR37 + nvb[1][3] = sqrt(AREF - data.deltaa) # - sqrt(a) (sqrt(m)) SFR10 + # BROADCAST ORBIT - 3 + nvb[2][0] = data.top # - Toe Time of Ephemeris (sec of NAVIC week) SFR10 + nvb[2][1] = data.cic # - Cic (radians) SFR11 + nvb[2][2] = data.omega0 * pi # - OMEGA0 (radians) SFR11 + nvb[2][3] = data.cis # - Cis (radians) SFR11 + # BROADCAST ORBIT - 4 + nvb[3][0] = data.i0 * pi # - i0 (radians) SFR11 + nvb[3][1] = data.crc # - Crc (meters) SFR11 + nvb[3][2] = data.omega * pi # - omega (radians) SFR10 + # nvb[3][3] = data.omegadot * pi # - OMEGA DOT (radians/sec) SFR37 + nvb[3][3] = ( + OMEGADOTREF - data.deltaomegadot + ) * pi # - OMEGA DOT (radians/sec) SFR11 + # BROADCAST ORBIT - 5 + nvb[4][0] = data.idot * pi # - IDOT (radians/sec) SFR11 + nvb[4][1] = data.deltan0 * pi # SFR10 + nvb[4][2] = data.uraned0 # - user range error NED0 SFRCLK + nvb[4][3] = data.uraned1 # - user range error NED1 SFRCLK + # BROADCAST ORBIT - 6 + nvb[5][0] = data.uraed # - user range error ED SFR10 + nvb[5][1] = ( + data.l1health | (data.l2health << 1) | (data.l5health << 2) + ) # - L1,L2,L5 health SFR10 + nvb[5][2] = data.tgd # - TGD (seconds) SFR30 + nvb[5][3] = data.uraned2 # user range error NED2 SFRCLK + # BROADCAST ORBIT - 7 + nvb[6][0] = data.iscl1ca # iono delay SFR30 + nvb[6][1] = data.iscl2c # iono delay SFR30 + nvb[6][2] = data.iscl5i5 # iono delay SFR30 + nvb[6][3] = data.iscl5q5 # iono delay SFR30 + # BROADCAST ORBIT - 8 + nvb[7][0] = data.tow # time of transmission TODO confirm value? + nvb[7][1] = data.wn # continuous week number SFR10 + nvb[7][2] = data.integrity | (data.l2phase << 1) | (data.alert << 2) + + def _convert_ionocorr(self, data: RawNav): """ - Format ionospheric correction header blocks. + Format ionospheric correction blocks. - :param Any data: data containing ionospheric corrections - :param str gnss: gnss code + RINEX 3 places these as IONOSPHERIC CORR header lines + RINEX 4 places these as ION Navigation record types + + :param RawNav data: data containing ionospheric corrections """ - if gnss == GPS: + msgtype = "LNAV" if data.sigid == "1C" else "CNAV" + tot = data.top if msgtype == "CNAV" else data.tot + if self._rinex_version < RINEX4: # RINEX 3.05 # timemark is tow converted to hour of day # and then to A-X character tm = chr(int((data.tow % 86400) / 3600) + 65) - self._ionocorr["GPSA"] = { - "parm1": data.alpha0, - "parm2": data.alpha1, - "parm3": data.alpha2, - "parm4": data.alpha3, - "timemark": tm, - "svid": data.svid, - } - self._ionocorr["GPSB"] = { - "parm1": data.beta0, - "parm2": data.beta1, - "parm3": data.beta2, - "parm4": data.beta3, - "timemark": tm, - "svid": data.svid, - } + self._ionocorr["GPSA"] = format_iono_corr( + corrtype="GPSA", + svid=data.svid, + timemark=tm, + parm1=data.alpha0, + parm2=data.alpha1, + parm3=data.alpha2, + parm4=data.alpha3, + ) + self._ionocorr["GPSB"] = format_iono_corr( + corrtype="GPSB", + svid=data.svid, + timemark=tm, + parm1=data.beta0, + parm2=data.beta1, + parm3=data.beta2, + parm4=data.beta3, + ) + else: # RINEX 4.02 + self._ionocorr[(data.svcode, data.epoch)] = format_ion( + svcode=data.svcode, + msgtype=msgtype, + msgsubtype="", + epoch=data.epoch, + a0=data.alpha0, + a1=data.alpha1, + a2=data.alpha2, + a3=data.alpha3, + b0=data.beta0, + b1=data.beta1, + b2=data.beta2, + b3=data.beta3, + ) - def _convert_timecorr(self, data: Any, gnss: str): + def _convert_eopcorr(self, data: RawNav): """ - Format time correction header blocks. + Format earth orientation correction blocks. - :param Any data: data containing time corrections - :param str gnss: gnss code + RINEX 4 places these as EOP Navigation record types + + TODO check epoch and repetition + + :param RawNav data: data containing eop corrections """ - if gnss == GPS: - self._timecorr["GPUT"] = { - "a0": data.a0, - "a1": data.a1, - "timeref": data.tot, - "weekno": data.wnt, - "svcode": f"{data.gnss}{data.svcode:02d}", - "source": 0, - } + if self._rinex_version < RINEX4: # # RINEX 3.05 + return + + msgtype = "LNAV" if data.sigid == "1C" else "CNAV" + self._eopcorr[(data.svcode, data.epoch)] = format_eop( + svcode=data.svcode, + msgtype=msgtype, + msgsubtype="", + epoch=data.epoch, + tom=data.teop, + xp=data.pmx, + dxpdt=data.pmxdot, + dxpdt2=0, + yp=data.pmy, + dypdt=data.pmydot, + dypdt2=0, + deltaut1=data.deltautgps, + ddeltaut1dt=data.deltautgpsdot, + d2deltaut1dt2=0, + ) - def _convert_rxmsfrbx(self, data: UBXMessage): + def _convert_timecorr(self, data: RawNav): """ - Process UBX RXM-SFRBX raw NAV subframe messages - using RawNav utility method. + Format time correction blocks. + + RINEX 3 places these as TIME SYSTEM CORR header lines + RINEX 4 places these as STO Navigation record types - :param UBXMessage data: UBX RXM-SFRBX message + TODO check epoch and repetition + + :param RawNav data: data containing time corrections + """ + + msgtype = "LNAV" if data.sigid == "1C" else "CNAV" + tot = data.top if msgtype == "CNAV" else data.tot + if self._rinex_version < RINEX4: # RINEX 3.05 + self._timecorr["GPUT"] = format_time_corr( + corrtype="GPUT", + svcode=data.svcode, + source="0", + timeref=tot, + weekno=data.wn, + a0=data.a0, + a1=data.a1, + ) + else: # RINEX 4.02 + self._timecorr[(data.svcode, data.epoch)] = format_sto( + svcode=data.svcode, + msgtype=msgtype, + msgsubtype="", + epoch=data.epoch, + timecode="GPUT", + sbasid="", + utcid="UTC(USNO)", + tot=tot, + a0=data.a0, + a1=data.a1, + a2=getattr(data, "a2", 0), + ) + + def _collate_rxmsfrbx(self, data: UBXMessage): + """ + Collate raw NAV subframes from UBX RXM-SFRBX message + using RawNav utility method. + + CURRENTLY ONLY GPS LNAV AND CNAV IMPLEMENTED, BUT + METHOD READILY EXTENSIBLE. + + :param UBXMessage data: UBX RXM-SFRBX message """ sfrdata = RawNav.process_rxm_sfrbx(data) - # NB: remove gnss and sigid constraints once other GNSS are - # transcribed - if sfrdata.get("gnss", "") != GPS or sfrdata.get("sigid", "") not in ( - "1C", + + # GPS LNAV + if sfrdata.get("gnss", "") == GPS and sfrdata.get("sigid", "") == "1C": + self._collate_rxmsfrbx_gps_lnav(sfrdata) + # GPS CNAV + elif sfrdata.get("gnss", "") == GPS and sfrdata.get("sigid", "") in ( "2L", + "2S", + "5I", + "5Q", ): - return + self._collate_rxmsfrbx_gps_cnav(sfrdata) + # GPS CNV2 + elif sfrdata.get("gnss", "") == GPS and sfrdata.get("sigid", "") in ( + "1S", + "1L", + ): + pass # TODO self._collate_rxmsfrbx_gps_cnv2(sfrdata) + # GALILEO FNAV/INAV + elif sfrdata.get("gnss", "") == GAL: + pass # TODO self._collate_rxmsfrbx_gal_fnav(sfrdata) + # GLONASS FDMA + elif sfrdata.get("gnss", "") == GLO and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_glo_fdma(sfrdata) + # GLONASS L1OC CDMA + elif ( + sfrdata.get("gnss", "") == GLO + and sfrdata.get("sigid", "") in ("??",) + and self._rinex_version >= RINEX4 + ): + pass # TODO self._collate_rxmsfrbx_glo_l1oc(sfrdata) + # GLONASS L3OC CDMA + elif ( + sfrdata.get("gnss", "") == GLO + and sfrdata.get("sigid", "") in ("??",) + and self._rinex_version >= RINEX4 + ): + pass # TODO self._collate_rxmsfrbx_glo_l3oc(sfrdata) + # QZSS LNAV + elif sfrdata.get("gnss", "") == QZS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_qzs_lnav(sfrdata) + # QZSS CNAV + elif sfrdata.get("gnss", "") == QZS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_qzs_cnav(sfrdata) + # QZSS CNV2 + elif sfrdata.get("gnss", "") == QZS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_qzs_cnv2(sfrdata) + # BEIDOU D1/D2 + elif sfrdata.get("gnss", "") == BDS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_bds_d1d2(sfrdata) + # BEIDOU CNV1 + elif sfrdata.get("gnss", "") == BDS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_bds_cnv1(sfrdata) + # BEIDOU CNV2 + elif sfrdata.get("gnss", "") == BDS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_bds_cnv2(sfrdata) + # BEIDOU CNV3 + elif sfrdata.get("gnss", "") == BDS and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_bds_cnv3(sfrdata) + # SBAS + elif sfrdata.get("gnss", "") == SBA: + pass # TODO self._collate_rxmsfrbx_sba(sfrdata) + # NAVIC LNAV + elif sfrdata.get("gnss", "") == IRN and sfrdata.get("sigid", "") in ("??",): + pass # TODO self._collate_rxmsfrbx_irn_lnav(sfrdata) + # NAVIC L1NV + elif ( + sfrdata.get("gnss", "") == IRN + and sfrdata.get("sigid", "") in ("??",) + and self._rinex_version >= RINEX4 + ): + pass # TODO self._collate_rxmsfrbx_irn_l1nv(sfrdata) + + def _collate_rxmsfrbx_gps_lnav( + self, sfrdata: dict[str, str | int | float | NoneType] + ): + """ + Collate raw GPS LNAV subframes from UBX RXM-SFRBX message + using RawNav utility method. + + :param dict[str, str | int | float | NoneType] sfrdata: raw subframe data + """ gnss = sfrdata["gnss"] svid = sfrdata["svid"] sigid = sfrdata["sigid"] subframeid = sfrdata["subframeid"] - svcode = sfrdata.get("svcode", 0) + pageid = sfrdata.get("pageid", 0) subframe = sfrdata["subframe"] try: - self._navframes[(gnss, svid)] = self._navframes.get( - (gnss, svid), RawNav(gnss, svid, sigid) + self._navframes[(gnss, svid, sigid)] = self._navframes.get( + (gnss, svid, sigid), RawNav(gnss, svid, sigid) ) - nav = self._navframes[(gnss, svid)] + nav = self._navframes[(gnss, svid, sigid)] if subframeid == 1: # clock parameters, sv health, etc. - nav.parse(subframe, GPS_LNAV_SUBFRAME_1) + nav.parse(subframe, GPS_LNAV_SUBFRAME_1, GPS_SFRACQ_MAP) elif subframeid == 2: # ephemerides - nav.parse(subframe, GPS_LNAV_SUBFRAME_2) + nav.parse(subframe, GPS_LNAV_SUBFRAME_2, GPS_SFRACQ_MAP) elif subframeid == 3: # ephemerides - nav.parse(subframe, GPS_LNAV_SUBFRAME_3) + nav.parse(subframe, GPS_LNAV_SUBFRAME_3, GPS_SFRACQ_MAP) elif subframeid == 4: - if svcode == 56: # page 18, ionospheric corrections - navit = RawNav(gnss, svid, sigid) - navit.parse(subframe, GPS_LNAV_SUBFRAME_4_P18) - self._convert_ionocorr(navit, gnss) - self._convert_timecorr(navit, gnss) - if nav.sfracq & 0b111 == TARGET_SFR: - self._convert_rawnav(self._navframes.pop((gnss, svid))) + if pageid == 56: # page 18, ionospheric & time corrections + nav.parse(subframe, GPS_LNAV_SUBFRAME_4_P18, GPS_SFRACQ_MAP) + self._convert_ionocorr(nav) + self._convert_timecorr(nav) + # when all the relevant subframes have been acquired, create a NAV record + if nav.sfracq & 0b111 == TARGET_SFR_GPS_LNAV: # 1,2,3 + self._convert_rawnav_gps_lnav(self._navframes.pop((gnss, svid, sigid))) except RINEXProcessingError as err: raise RINEXProcessingError("Error process RXM-SFRBX data") from err - def _get_external_epoch(self, rt: str) -> datetime: + def _collate_rxmsfrbx_gps_cnav( + self, sfrdata: dict[str, str | int | float | NoneType] + ): """ - Get epoch from external source (e.g. other NAV - message in same data stream). + Collate raw GPS CNAV subframes from UBX RXM-SFRBX message + using RawNav utility method. - :param str rt: RINEX observation type (OBS/NAV) - :return: GNSS epoch - :rtype: datetime + :param dict[str, str | int | float | NoneType] sfrdata: raw subframe data """ - return self.__app.get_current_epoch(rt) + gnss = sfrdata["gnss"] + svid = sfrdata["svid"] + sigid = sfrdata["sigid"] + subframeid = sfrdata["subframeid"] + subframe = sfrdata["subframe"] + + try: + self._navframes[(gnss, svid, sigid)] = self._navframes.get( + (gnss, svid, sigid), RawNav(gnss, svid, sigid) + ) + nav = self._navframes[(gnss, svid, sigid)] + if subframeid == 10: # ephemeris 1 + nav.parse(subframe, GPS_CNAV_SUBFRAME_10, GPS_SFRACQ_MAP) + elif subframeid == 11: # ephemeris 2 + nav.parse(subframe, GPS_CNAV_SUBFRAME_11, GPS_SFRACQ_MAP) + elif subframeid == 30: # clock, iono & group delay + nav.parse(subframe, GPS_CNAV_SUBFRAME_30, GPS_SFRACQ_MAP) + self._convert_ionocorr(nav) + elif subframeid == 32 and self._rinex_version >= RINEX4: # clock, eop + nav.parse(subframe, GPS_CNAV_SUBFRAME_32, GPS_SFRACQ_MAP) + self._convert_eopcorr(nav) + elif subframeid == 33: # clock, utc + nav.parse(subframe, GPS_CNAV_SUBFRAME_33, GPS_SFRACQ_MAP) + self._convert_timecorr(nav) + # when all the relevant subframes have been acquired, create a NAV record + if nav.sfracq & 0b111 == TARGET_SFR_GPS_CNAV: # 10,11,30 + self._convert_rawnav_gps_cnav(self._navframes.pop((gnss, svid, sigid))) + except RINEXProcessingError as err: + raise RINEXProcessingError("Error process RXM-SFRBX data") from err - def _get_epoch( - self, wno: int, toc: int, gnss: Literal["G", "E", "C", "J", "I"] - ) -> datetime: + def _get_external_epoch(self, rt: str) -> datetime: """ - Get epoch from message wno, tow. + Get epoch from external source (e.g. other NAV + message in same data stream). - :param int wno: week number in GNSS time system - :param int toc: time of week in seconds in GNSS time system - :param Literal["G","E","C","J","I"] gnss: GNSS time system to use + :param str rt: RINEX observation type (OBS/NAV) :return: GNSS epoch :rtype: datetime """ - epoch = wnotow2utc( - wno=wno, - tow=int(toc * 1000), - ls=0, # use GPS time, not UTC time - gnss=gnss, - autoroll=True, - modwno=False, - ) - self.__app.set_current_epoch(epoch, NAV) - return epoch + return self.__app.get_current_epoch(rt) def _get_rmc_epoch(self, data: NMEAMessage) -> datetime | NoneType: """ diff --git a/src/pygnssutils/rinex_conv_obs.py b/src/pygnssutils/rinex_conv_obs.py index 8f80d40..187fce0 100644 --- a/src/pygnssutils/rinex_conv_obs.py +++ b/src/pygnssutils/rinex_conv_obs.py @@ -25,7 +25,7 @@ from logging import getLogger from typing import Any, Literal -from pynmeagps import NMEAMessage, llh2ecef, wnotow2utc +from pynmeagps import NMEAMessage, llh2ecef from pyrtcm import RTCMMessage from pyubx2 import UBXMessage @@ -69,6 +69,7 @@ format_sys_phaseshift, format_sys_scalefactor, format_timefirstlast, + get_epoch, get_obscode, get_svcode_ubx, ) @@ -219,6 +220,7 @@ def _format_observation_epoch( numobs: int | str, epochflag: int | str = "0", clkoffset: float | str = "", + picosecond: int | str = "", ) -> str: """ Format observation epoch. @@ -227,6 +229,7 @@ def _format_observation_epoch( :param int | str numobs: number of observations in this epoch :param int | str epochflag: epoch flag :param float | str clkoffset: clock offset + :param int | str picosets: picoseconds (RINEX 4 only) :return: formatted string :rtype: str """ @@ -237,12 +240,16 @@ def _format_observation_epoch( if epoch == "": # event return f">{'':>39}{epochflag:>3}{numobs:>3}{'':>21}\n" # A1 ... 2X,I1 I3 # epoch + if isinstance(picosecond, int): + pico = f"{picosecond:05d}" + else: + pico = f"{picosecond:>5}" return ( f">{epoch.year:>5}{epoch.month:>3}{epoch.day:>3}" f"{epoch.hour:>3}{epoch.minute:>3}" f"{FRNX(epoch.second + epoch.microsecond/1000000,11,7)}" - f"{epochflag:>3}{nummeas:>3}{'':>6}{FRNX(clkoffset,15,12)}\n" - ) # A1 1X,I4 4(1X,I2.2) F11.7 2X,I1 I3 6X F15.12 + f"{epochflag:>3}{nummeas:>3}{'':>6}{FRNX(clkoffset,15,12)} {pico}\n" + ) # A1 1X,I4 4(1X,I2.2) F11.7 2X,I1 I3 6X F15.12 (1X,I5.5) def _format_obs_data(self, obsdata: dict[datetime, dict] | str = ""): """ @@ -412,7 +419,6 @@ def get_nmea_pos(self, data: NMEAMessage): self._approxpos = [x, y, z] except (AttributeError, TypeError) as err: raise (err) from err - # print(f"something went wrong {err}") def convert_ubx_rxmrawx(self, data: UBXMessage): """ @@ -427,9 +433,7 @@ def convert_ubx_rxmrawx(self, data: UBXMessage): def geta(att: str, i: int): return getattr(data, f"{att}_{i+1:02d}") - epoch = wnotow2utc( - wno=data.week, tow=int(data.rcvTow * 1000), ls=None, gnss=GPS, autoroll=True - ) + epoch, _ = get_epoch(data.week, data.rcvTow, GPS) if epoch != self.__app.get_current_epoch(OBS): self.__app.set_current_epoch(epoch, OBS) self._obsdata[epoch] = {} diff --git a/src/pygnssutils/rinex_globals.py b/src/pygnssutils/rinex_globals.py index 13864b0..efccd50 100644 --- a/src/pygnssutils/rinex_globals.py +++ b/src/pygnssutils/rinex_globals.py @@ -42,7 +42,8 @@ RINEX_ERROR = 99 RINEX_NORECS = 1 RINEX_OK = 0 -RINEXVERSIONS = ["3.05", "4.02"] # only 3.05 supported in Alpha release +RINEX4 = "4.00" +RINEXVERSIONS = ["3.05", "4.02"] RINEXVER_DEFAULT = RINEXVERSIONS[0] SBA = "S" TIME_BEIDOU = "BDT" @@ -57,22 +58,41 @@ """RINEX File Types.""" # scaling factors +P2_N4 = 0.0625 # 2**-4 P2_N5 = 0.03125 # 2**-5 +P2_N6 = 0.015625 # 2**-6 +P2_N8 = 0.00390625 # 2**-8 +P2_N9 = 0.001953125 # 2**-9 +P2_N14 = 6.103515625e-05 # 2**-14 +P2_N15 = 3.0517578125e-05 # 2**-15 +P2_N16 = 1.52587890625e-05 # 2**-16 P2_N19 = 1.9073486328125e-06 # 2**-19 P2_N20 = 9.5367431640625e-07 # 2**-20 P2_N21 = 4.76837158203125e-07 # 2**-21 P2_N23 = 1.1920928955078125e-07 # 2**-23 P2_N24 = 5.960464477539063e-08 # 2**-24 +P2_N25 = 2.9802322387695312e-08 # 2**-25 P2_N27 = 7.450580596923828e-09 # 2**-27 P2_N29 = 1.862645149230957e-09 # 2**-29 P2_N30 = 9.313225746154785e-10 # 2**-30 P2_N31 = 4.656612873077393e-10 # 2**-31 +P2_N32 = 2.3283064365386963e-10 # 2**-32 P2_N33 = 1.1641532182693481e-10 # 2**-33 +P2_N34 = 5.820766091346741e-11 # 2**-34 +P2_N35 = 2.9103830456733704e-11 # 2**-35 +P2_N37 = 7.275957614183426e-12 # 2**-37 P2_N38 = 3.637978807091713e-12 # 2**-38 P2_N43 = 1.1368683772161603e-13 # 2**-43 +P2_N44 = 5.684341886080802e-14 # 2**-44 +P2_N48 = 3.552713678800501e-15 # 2**-48 P2_N50 = 8.881784197001252e-16 # 2**-50 +P2_N51 = 4.440892098500626e-16 # 2**-51 P2_N55 = 2.7755575615628914e-17 # 2**-55 +P2_N57 = 6.938893903907228e-18 # 2**-57 +P2_N60 = 8.673617379884035e-19 # 2**-60 +P2_N68 = 3.3881317890172014e-21 # 2**-68 P2_P4 = 16 # 2**4 +P2_P9 = 512 # 2**9 P2_P11 = 2048 # 2**11 P2_P12 = 4096 # 2**12 P2_P14 = 16384 # 2**14 @@ -162,40 +182,42 @@ """ UBXRINEXOBSCODE = { - (0, 0): "1C", # L1_C/A", - (0, 3): "2L", # "L2_CL", - (0, 4): "2S", # "L2_CM", - (0, 6): "5I", # "L5_I", - (0, 7): "5Q", # "L5_Q", - (1, 0): "1C", # "L1_C/A", - (2, 0): "1C", # "E1_C", - (2, 1): "1B", # "E1_B", - (2, 3): "5I", # "E5_aI", - (2, 4): "5Q", # "E5_aQ", - (2, 5): "7I", # "E5_bI", - (2, 6): "7Q", # "E5_bQ", - (2, 8): "6B", # "E6_B", - (2, 9): "6C", # "E6_C", - (3, 0): "2I", # "B1I_D1", - (3, 1): "1C", # "B1I_D2", - (3, 2): "1C", # "B2I_D1", - (3, 3): "2C", # "B2I_D2", - (3, 4): "6C", # "B3I_D1", - (3, 5): "1P", # "B1_Cp", - (3, 6): "1D", # "B1_Cd", - (3, 7): "5P", # "B2_ap", - (3, 8): "5D", # "B2_ad", - (3, 10): "6I", # "B3I_D2", - (5, 0): "1C", # "L1_C/A", - (5, 1): "1Z", # "L1_S", - (5, 4): "2S", # "L2_CM", - (5, 5): "2L", # "L2_CL", - (5, 8): "5I", # "L5_I", - (5, 9): "5Q", # "L5_Q", - (5, 12): "1B", # "L1_CB", - (6, 0): "1C", # "L1_OF", - (6, 2): "2C", # "L2_OF", - (7, 0): "5A", # "L5_A", + (0, 0): "1C", # GPS L1 C/A Legacy LNAV + (0, 1): "1S", # GPS L1C D Data code CNV2 (not implemented by u-blox) + (0, 2): "1L", # GPS L1C P Pilot code CNV2 (not implemented by u-blox) + (0, 3): "2L", # GPS L2C L Civil Long-length code CNAV + (0, 4): "2S", # GPS L2C M Civil Moderate code CNAV + (0, 6): "5I", # GPS L5 I In-phase code CNAV + (0, 7): "5Q", # GPS L5 Q Quadrature code CNAV + (1, 0): "1C", # SBA L1 C/A + (2, 0): "1C", # GAL E1_C + (2, 1): "1B", # GAL E1_B + (2, 3): "5I", # GAL E5_aI + (2, 4): "5Q", # GAL E5_aQ + (2, 5): "7I", # GAL E5_bI + (2, 6): "7Q", # GAL E5_bQ + (2, 8): "6B", # GAL E6_B + (2, 9): "6C", # GAL E6_C + (3, 0): "2I", # BDS B1I_D1 + (3, 1): "1C", # BDS B1I_D2 + (3, 2): "1C", # BDS B2I_D1 + (3, 3): "2C", # BDS B2I_D2 + (3, 4): "6C", # BDS B3I_D1 + (3, 5): "1P", # BDS B1_Cp + (3, 6): "1D", # BDS B1_Cd + (3, 7): "5P", # BDS B2_ap + (3, 8): "5D", # BDS B2_ad + (3, 10): "6I", # BDS B3I_D2 + (5, 0): "1C", # QZS L1_C/A + (5, 1): "1Z", # QZS L1_S + (5, 4): "2S", # QZS L2_CM + (5, 5): "2L", # QZS L2_CL + (5, 8): "5I", # QZS L5_I + (5, 9): "5Q", # QZS L5_Q + (5, 12): "1B", # QZS L1_CB + (6, 0): "1C", # GLO L1_OF + (6, 2): "2C", # GLO L2_OF + (7, 0): "5A", # IRN L5_A } """ UBX Signal ID -> RINEX Observation Code Lookup. TODO CHECK THIS MAPPING!!! @@ -324,6 +346,41 @@ - value is (Frequency Band, Signal, Phase Alignment Frequency Band) """ +EPHNAVTYPES = { + (BDS, "B1C"): "CNV1", # B1C + (BDS, "B1I_GEO"): "D2", # B1I + (BDS, "B1I_MEOIGSO"): "D1", # B1I + (BDS, "B2a"): "CNV2", # B2a + (BDS, "B2b"): "CNV3", # B2b + (BDS, "B2I_GEO"): "D2", # B2I + (BDS, "B2I_MEOIGSO"): "D1", # B2I + (BDS, "B3I_GEO"): "D2", # B3I + (BDS, "B3I_MEOIGSO"): "D1", # B3I + (GAL, "E1"): "INAV", # E1 + (GAL, "E5a"): "FNAV", # E5a + (GAL, "E5b"): "INAV", # E5b + (GLO, "L1 C/A"): "FDMA", # L1 C/A + (GLO, "L1OC"): "L1OC", # L1OC + (GLO, "L3OC"): "L3OC", # L3OC + (GPS, "L1 C/A"): "LNAV", # L1 C/A + (GPS, "L1C"): "CNV2", # L1C + (GPS, "L2C"): "CNAV", # L2C + (GPS, "L5"): "CNAV", # L5 + (IRN, "L1"): "L1NV", # L1 + (IRN, "L5/S SPS"): "LNAV", # L5/S SPS + (QZS, "L1 C/A"): "LNAV", # L1 C/A + (QZS, "L1 C/B"): "LNAV", # L1 C/B + (QZS, "L1C"): "CNV2", # L1C + (QZS, "L2C"): "CNAV", # L2C + (QZS, "L5"): "CNAV", # L5 + (SBA, "L1"): "SBAS", # L1 +} +""" +EPH Navigation Message Types. + +- key is (RINEX GNSS Code, Signal Code) +""" + RINEX_METOBS = { "PR": "Pressure (mbar)", "TD": "Dry temperature (deg Celsius)", diff --git a/src/pygnssutils/rinex_helpers.py b/src/pygnssutils/rinex_helpers.py index 35466bb..cb4efdf 100644 --- a/src/pygnssutils/rinex_helpers.py +++ b/src/pygnssutils/rinex_helpers.py @@ -20,7 +20,7 @@ from types import NoneType from typing import Literal -from pynmeagps import leapsecond +from pynmeagps import leapsecond, utc2wnotow, wnotow2utc from pygnssutils.rinex_globals import ( BDS, @@ -80,6 +80,32 @@ def DRNX(num: float | int | str, length: int, sig: int) -> str: return f"{num:>{length}}" +def get_epoch( + wno: int, tow: int, gnss: Literal["G", "E", "C", "J", "I"] +) -> tuple[datetime, int]: + """ + Get epoch and non-modular week number for given modular wno and tow. + + :param int wno: modular week number + :param int two: time of week in seconds + :param Literal['G', 'E', 'C', 'J', 'I'] gnss: gnss code + :return epoch, non-modular wno + :return tuple[datetime, int] + """ + + epoch = wnotow2utc( + wno=wno, + tow=int(tow * 1000), + ls=None, + gnss=gnss, + autoroll=True, + modwno=True, + ) + # convert week number to non-modular + wn, _, _ = utc2wnotow(utc=epoch, gnss=gnss, modwno=False) + return epoch, wn + + def get_fithours(iodc: int, fit: int, gnss: str) -> int | str: """ Get FIT interval in hours for given IODC and fit flag. @@ -195,6 +221,11 @@ def format_filename( interval: int | float, outputpath: Path = Path("."), source: Literal["R", "S", "N", "U"] = "R", + form: Literal["IGS", "GNSS"] = "GNSS", + site: str = "SITE", + marker: int = 0, + receiver: int = 0, + country: str = "USA", ) -> Path: """ Format output file name using RINEX long filename format. @@ -208,10 +239,20 @@ def format_filename( :param int | float interval: observation interval in seconds :param Path | str outputpath: fully-qualified file path (".") :param Literal["R","S","N","U"] source: source of observations (R = Receiver) + :param Literal["IGS","GNSS"] form; filename format to use + :param str site: site (for IGS format only) + :param int marker: marker number (for IGS format only) + :param int receiver: receiver number (for IGS format only) + :param str country: country code (for IGS format only) :return: fully-qualified output file path :rtype: Path """ + if form == "IGS": + station = f"{site:<4}{marker}{receiver}{country:<3}" + else: + station = "pygpsdata" + if startepoch == EPOCHMAX or endepoch == EPOCHMIN: period = TIME_UNDEFINED else: @@ -232,7 +273,7 @@ def format_filename( rtu = rinextype.upper() src = "S" if source == "N" else source # treat NTRIP as Stream return Path.joinpath( - outputpath, f"pygpsdata_{src}_{start}_{period}_{frequency}{gnu}{rtu}.rnx" + outputpath, f"{station}_{src}_{start}_{period}_{frequency}{gnu}{rtu}.rnx" ) @@ -375,13 +416,12 @@ def format_approxpos(approxpos: list[float] | str = "") -> str: :rtype: str """ - if approxpos != "": - x, y, z = approxpos - return ( - f"{FRNX(x,14,4)}{FRNX(y,14,4)}" - f"{FRNX(z,14,4)}{'':<18}APPROX POSITION XYZ\n" - ) # 3F14.4 - return "" + if approxpos == "": + return "" + x, y, z = approxpos + return ( + f"{FRNX(x,14,4)}{FRNX(y,14,4)}" f"{FRNX(z,14,4)}{'':<18}APPROX POSITION XYZ\n" + ) # 3F14.4 def format_antennadeltahen( @@ -417,13 +457,12 @@ def format_antennadeltaxyz(deltaxyz: list[float] | str = "") -> str: :rtype: str """ - if deltaxyz != "": - x, y, z = deltaxyz - return ( - f"{FRNX(x,14,4)}{FRNX(y,14,4)}" - f"{FRNX(z,14,4)}{'':<18}ANTENNA: DELTA X/Y/Z\n" - ) # 3F14.4 - return "" + if deltaxyz == "": + return "" + x, y, z = deltaxyz + return ( + f"{FRNX(x,14,4)}{FRNX(y,14,4)}" f"{FRNX(z,14,4)}{'':<18}ANTENNA: DELTA X/Y/Z\n" + ) # 3F14.4 def format_sys_antennaphasecentre( @@ -475,13 +514,12 @@ def format_antennabsight(bsight: list[float] | str = "") -> str: :rtype: str """ - if bsight != "": - x, y, z = bsight - return ( - f"{FRNX(x,14,4)}{FRNX(y,14,4)}{FRNX(z,14,4)}" - f"{'':<18}ANTENNA: B.SIGHT XYZ\n" - ) # 3F14.4 - return "" + if bsight == "": + return "" + x, y, z = bsight + return ( + f"{FRNX(x,14,4)}{FRNX(y,14,4)}{FRNX(z,14,4)}" f"{'':<18}ANTENNA: B.SIGHT XYZ\n" + ) # 3F14.4 def format_antennazerodirazi(azi: float | str = "") -> str: @@ -493,9 +531,9 @@ def format_antennazerodirazi(azi: float | str = "") -> str: :rtype: str """ - if azi != "": - return f"{FRNX(azi,14,4)}{'':<46}ANTENNA: ZERODIR AZI\n" # F14.4 - return "" + if azi == "": + return "" + return f"{FRNX(azi,14,4)}{'':<46}ANTENNA: ZERODIR AZI\n" # F14.4 def format_antennazerodirxyz(zerodir: list[float] | str = "") -> str: @@ -507,13 +545,12 @@ def format_antennazerodirxyz(zerodir: list[float] | str = "") -> str: :rtype: str """ - if zerodir != "": - x, y, z = zerodir - return ( - f"{FRNX(x,14,4)}{FRNX(y,14,4)}" - f"{FRNX(z,14,4)}{'':<18}ANTENNA: ZERODIR XYZ\n" - ) # 3F14.4 - return "" + if zerodir == "": + return "" + x, y, z = zerodir + return ( + f"{FRNX(x,14,4)}{FRNX(y,14,4)}" f"{FRNX(z,14,4)}{'':<18}ANTENNA: ZERODIR XYZ\n" + ) # 3F14.4 def format_centermass(centermass: list[float] | str = "") -> str: @@ -526,13 +563,12 @@ def format_centermass(centermass: list[float] | str = "") -> str: :rtype: str """ - if centermass != "": - x, y, z = centermass - return ( - f"{FRNX(x,14,4)}{FRNX(y,14,4)}" - f"{FRNX(z,14,4)}{'':<18}CENTER OF MASS: XYZ\n" - ) # 3F14.4 - return "" + if centermass == "": + return "" + x, y, z = centermass + return ( + f"{FRNX(x,14,4)}{FRNX(y,14,4)}" f"{FRNX(z,14,4)}{'':<18}CENTER OF MASS: XYZ\n" + ) # 3F14.4 def format_obstypes(obstypes: dict | str = "") -> str: @@ -597,9 +633,9 @@ def format_interval(interval: float | str = "") -> str: :rtype: str """ - if interval != "": - return f"{FRNX(interval, 10, 3)}{'':<50}INTERVAL\n" # F10.3 - return "" + if interval == "": + return "" + return f"{FRNX(interval, 10, 3)}{'':<50}INTERVAL\n" # F10.3 def format_obstime(utc: datetime, source: str = "GPS") -> str: @@ -643,9 +679,9 @@ def format_clockoffset(offset: int | str = "") -> str: :rtype: str """ - if offset not in (0, ""): - return f"{offset:<6}{'':<54}RCV CLOCK OFFS APPL\n" # I6 - return "" + if offset in (0, ""): + return "" + return f"{offset:<6}{'':<54}RCV CLOCK OFFS APPL\n" # I6 def format_sys_dcbsapplied(dcbsapplied: dict | str = "") -> str: @@ -914,85 +950,193 @@ def format_numsats(numsats: int | str = "") -> str: return "" -def format_iono_corr(ionocorr: dict | str = "") -> str: +def format_iono_corr( + svid: int, + corrtype: str, + timemark: str, + parm1: float, + parm2: float, + parm3: float, + parm4: float, +) -> str: """ - Format Ionospheric Corrections. - - Format of ionocorr dict:: + Format Ionospheric Corrections (RINEX 3). - ionocorr = { - corrtype (str) : { - "parm1" : parm1 (float), - "parm2" : parm2 (float), - "parm3" : parm3 (float), - "parm4" : parm4 (float), - "timemark" : timemark (str), - "svid" : svid (str), - }, - } - - :param dict | str ionocorr: ionospheric correction dictionary + :param str svid: SV id + :param str corrtype: correction type + :param str timemark: time mark + :param float parm1: ionospheric correction parameter 1 + :param float parm2: ionospheric correction parameter 2 + :param float parm3: ionospheric correction parameter 3 + :param float parm4: ionospheric correction parameter 4 :return: formatted string :rtype: str """ - if ionocorr == "": - ionocorr = {} + # A4,1X 4D12.4 1X,A1 1X,I2 + out = ( + f"{corrtype:<4} {DRNX(parm1,12,4)}{DRNX(parm2,12,4)}{DRNX(parm3,12,4)}" + f"{DRNX(parm4,12,4)} {timemark} {svid:02d} IONOSPHERIC CORR\n" + ) + return out - out = "" - for corrtype, values in ionocorr.items(): - parm1 = values["parm1"] - parm2 = values["parm2"] - parm3 = values["parm3"] - parm4 = values["parm4"] - timemark = values["timemark"] - svid = values["svid"] - # A4,1X 4D12.4 1X,A1 1X,I2 - out += ( - f"{corrtype:<4} {DRNX(parm1,12,4)}{DRNX(parm2,12,4)}{DRNX(parm3,12,4)}" - f"{DRNX(parm4,12,4)} {timemark} {svid:02d} IONOSPHERIC CORR\n" - ) + +def format_ion( + svcode: str, + msgtype: str, + msgsubtype: str, + epoch: datetime, + a0: float, + a1: float, + a2: float, + a3: float, + b0: float, + b1: float, + b2: float, + b3: float, +) -> str: + """ + Format Ionospheric Correction record (RINEX 4). + + param str svcode: SV code + param str msgtype: message type + param str msgsubtype: sub message type + param datetime epoch: epoch + param float a0: ionospheric correction coefficient a0 + param float a1: ionospheric correction coefficient a1 + param float a2: ionospheric correction coefficient a2 + param float a3: ionospheric correction coefficient a3 + param float b0: ionospheric correction coefficient b0 + param float b1: ionospheric correction coefficient b1 + param float b2: ionospheric correction coefficient b2 + param float b3: ionospheric correction coefficient b3 + """ + + # A1 1X,A3 1X,A1 A2 1X,A4 1X,A4 + # 4X,I4, 5(1X,I2.2), 3E19.12 + # 4X,4E19.12 + # 4X,E19.12 + out = ( + f"> ION {svcode:<3} {msgtype:<4} {msgsubtype:<4}\n" + f" {epoch.year:04d} {epoch.month:02d} {epoch.day:02d} {epoch.hour:02d} " + f"{epoch.minute:02d} {epoch.second:02d}{DRNX(a0, 19,12)}{DRNX(a1, 19,12)}" + f"{DRNX(a2, 19,12)}\n" + f" {DRNX(a3, 19,12)}{DRNX(b0, 19,12)}{DRNX(b1, 19,12)}{DRNX(b2, 19,12)}\n" + f" {DRNX(b3, 19,12)}\n" + ) return out -def format_time_corr(timecorr: dict | str = "") -> str: +def format_eop( + svcode: str, + msgtype: str, + msgsubtype: str, + epoch: datetime, + tom: float, + xp: float, + dxpdt: float, + dxpdt2: float, + yp: float, + dypdt: float, + dypdt2: float, + deltaut1: float, + ddeltaut1dt: float, + d2deltaut1dt2: float, +) -> str: """ - Format Time System Corrections. + Format Earth Orientation record (RINEX 4). - Format of timecorr dict:: + param str svcode: SV code + param str msgtype: message type + param str msgsubtype: sub message type + param datetime epoch: epoch - timecorr = { - corrtype (str) : { - "a0" : a0 (float), - "a1" : a1 (float), - "timeref" : timeref (int), - "weekno" : weekno (int), - "svcode" : svcode (str), - "source" : source (str), - }, - } + """ + + # A1 1X,A3 1X,A1 A2 1X,A4 1X,A4 + # 4X,I4, 5(1X,I2.2), 3E19.12 + # 4X,A19, 3E19.12 + # 4X,4E19.12 + out = ( + f"> EOP {svcode:<3} {msgtype:<4} {msgsubtype:<4}\n" + f" {epoch.year:04d} {epoch.month:02d} {epoch.day:02d} {epoch.hour:02d} " + f"{epoch.minute:02d} {epoch.second:02d}{DRNX(xp, 19,12)}{DRNX(dxpdt, 19,12)}" + f"{DRNX(dxpdt2, 19,12)}\n" + f" {'':<19}{DRNX(yp, 19,12)}{DRNX(dypdt, 19,12)}{DRNX(dypdt2, 19,12)}\n" + f" {DRNX(tom,19,12)}{DRNX(deltaut1, 19,12)}{DRNX(ddeltaut1dt, 19,12)}{DRNX(d2deltaut1dt2, 19,12)}\n" + ) + return out - :param dict | str timecorr: time corrections dictionary + +def format_time_corr( + svcode: str, + corrtype: str, + timeref: int, + weekno: int, + source: str, + a0: float, + a1: float, +) -> str: + """ + Format Time System Offset (RINEX 3). + + :param str svcode: SV code + :param str corrtype: correction type + :param int timeref: time reference + :param int weekno: week number + :param str source: time source + :param float a0: a0 clock offset + :param float a1: a1 clock offset :return: formatted string :rtype: str """ - if timecorr == "": - timecorr = {} + # A4,1X D17.10 D16.9 1XI6 1XI4 1X,A5,1X I2,1X + out = ( + f"{corrtype:<4} {DRNX(a0,17, 10)}{DRNX(a1,16, 9)} " + f"{timeref:>6} {weekno:>4} {svcode:>5} {source:>2} TIME SYSTEM CORR\n" + ) + return out - out = "" - for corrtype, values in timecorr.items(): - a0 = values["a0"] - a1 = values["a1"] - timeref = values["timeref"] - weekno = values["weekno"] - svcode = values["svcode"] - source = values["source"] - # A4,1X D17.10 D16.9 1XI6 1XI4 1X,A5,1X I2,1X - out += ( - f"{corrtype:<4} {DRNX(a0,17, 10)}{DRNX(a1,16, 9)} " - f"{timeref:<6} {weekno:<4} {svcode:>5} {source:>2} TIME SYSTEM CORR\n" - ) + +def format_sto( + svcode: str, + msgtype: str, + msgsubtype: str, + epoch: datetime, + timecode: str, + sbasid: str, + utcid: str, + tot: float, + a0: float, + a1: float, + a2: float, +) -> str: + """ + Format System Time Offset (RINEX 4). + + param str svcode: SV code + param str msgtype: message type + param str msgsubtype: sub message type + param datetime epoch: epoch + param str timecode: timecode + param str sbasid: sbas id (blank if not sbas) + param str utcid: UTC id + param float tot: time of transmission + param float a0: clock offset coefficient a0 + param float a1: clock offset coefficient a1 + param float a2: clock offset coefficient a2 + """ + + # A1 1X,A3 1X,A1 A2 1X,A4 1X,A4 + # 4X,I4, 5(1X,I2.2), 1X,A18 (left justified) 1X,A18 (left justified) 1X,A18 (left justified) + # 4X,4E19.12 + out = ( + f"> STO {svcode:<3} {msgtype:<4} {msgsubtype:<4}\n" + f" {epoch.year:04d} {epoch.month:02d} {epoch.day:02d} {epoch.hour:02d} " + f"{epoch.minute:02d} {epoch.second:02d} {timecode:<18} {sbasid:<18} {utcid:<18}\n" + f" {DRNX(tot, 19,12)}{DRNX(a0, 19,12)}{DRNX(a1, 19,12)}{DRNX(a2, 19,12)}\n" + ) return out @@ -1079,6 +1223,70 @@ def format_met_sensorpos(senspos: list[float] | str, obstype: str = "PR") -> str ) +def format_doi(doi: str = "") -> str: + """ + Format Digital Object Identifier (DOI). + + e.g. "https://doi.org/10.1000/182" + + :param str doi: digital object identifier + :return: formatted string + :rtype: str + """ + + if doi == "": + return "" + return f"{doi:<{DATAWIDTH}}DOI\n" + + +def format_licenseofuse(lou: str = "") -> str: + """ + Format License of Use. + + e.g. "CC BY 04 ; https://creativecommons.org/licenses/by/4.0/" + + :param str lou: license of use descriptor + :return: formatted string + :rtype: str + """ + + if lou == "": + return "" + return f"{license:<{DATAWIDTH}}LICENSE OF USE\n" + + +def format_stationinfo(station: str = "") -> str: + """ + Format Station Information. + + :param str station: station info + :return: formatted string + :rtype: str + """ + + if station == "": + return "" + return f"{station:<{DATAWIDTH}}STATION INFORMATION\n" + + +def format_nav_typesvmssg( + rectyp: str, sv: str, msgtyp: str, msgsubtyp: str = "" +) -> str: + """ + Format navigation record type, SV and message type/subtype. + + :param str rectyp: record type e.g. EPH + :param str sv: SV e.g. G04 + :param str msgtyp: message type e.g. LNAV + :param str msgsubtyp: message subtype + :return: formatted string + :rtype: str + """ + + # A1 1X,A3 1X,A3 1X,A4 1X,A4 + return f"> {rectyp:<3} {sv:<3} {msgtyp:<4} {msgsubtyp:<4}\n" + + def format_headerend() -> str: """ Format header end. diff --git a/src/pygnssutils/rinex_subframes_gps.py b/src/pygnssutils/rinex_subframes_gps.py index 2a23d28..ff8f628 100644 --- a/src/pygnssutils/rinex_subframes_gps.py +++ b/src/pygnssutils/rinex_subframes_gps.py @@ -1,8 +1,6 @@ """ rinex_subframes_gps.py -!!! WORK IN PROGRESS !!! - GPS NAV Subframe definitions. https://archive.gps.gov/technical/icwg/IS-GPS-200N.pdf @@ -19,19 +17,21 @@ NOTE: if extracting raw subframe data from UBX RXM-SFRBX messages, note that the native 30-bit GPS words are padded to 32-bit dwrds `BUT` the padding -treatment depends on the signal type. For L1 (LNAV) signals, the final +treatment depends on the signal type. For L1 C/A (LNAV) signals, the final 2 bits of each 32-bit dwrd are padding and can simply be stripped off, -but for other signal types things are more complicated, and the relevant -treatment does not appear to be specified in any public domain UBX protocol -specification. +but for other signal types (L2C / L5 (CNAV), etc.) things are more complicated, +and the relevant treatment does not appear to be documented in any public domain +UBX protocol specification. There have been attempts to 'reverse engineer' the treatment. See, for example: https://portal.u-blox.com/s/question/0D52p00008HKD1kCAH/why-are-the-sfrbx-messages-words-32-bits-but-in-isgps200h-the-words-are-specified-as-being-30-bits-long https://github.com/semuconsulting/pyubx2/blob/master/src/pyubx2/ubxtypes_get.py -GPS L1 : pads 2 bits at the end of each word (30 bit message + 2 bit padding) x 10 words = 320 bits total -GPS L5: pads 20 bits at end (300 bit message + 20 bit padding using the next frame) = 320 bits total +GPS L1 : pads 2 bits at the end of each word (30 bit message + 2 bit \ + padding) x 10 words = 320 bits total +GPS L5: pads 20 bits at end (300 bit message + 20 bit padding using \ + the next frame) = 320 bits total Created on 6 Oct 2025 @@ -44,29 +44,50 @@ from pygnssutils.rawnav import PREAMBLE, SFR, TOC, TOW, VALPREAMBLE, WN, S, U from pygnssutils.rinex_globals import ( + P2_N4, P2_N5, + P2_N6, + P2_N8, + P2_N9, + P2_N14, + P2_N15, + P2_N16, P2_N19, P2_N20, P2_N21, P2_N23, P2_N24, + P2_N25, P2_N27, P2_N29, P2_N30, P2_N31, + P2_N32, P2_N33, + P2_N34, + P2_N35, + P2_N37, P2_N38, P2_N43, + P2_N44, + P2_N48, P2_N50, + P2_N51, P2_N55, + P2_N57, + P2_N60, + P2_N68, P2_P4, + P2_P9, P2_P11, P2_P12, P2_P14, P2_P16, ) -# TODO complete scaling factors for all subframe definitions +# ********************************************************************** +# LNAV - L1 C/A +# ********************************************************************** # attribute_name: (bit offset, bit length, bit encoding, scaling) GPS_LNAV_TLM = { @@ -217,7 +238,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "_word3_10": (68, 232, U, 0), } @@ -228,7 +249,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "e": (68, 16, U, P2_N21), "_parity3": (84, 6, U, 0), # word4 @@ -266,7 +287,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "toa": (68, 8, U, 0), "wna": (76, 8, U, 0), "_parity3": (84, 6, U, 0), @@ -321,7 +342,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "_reserved2": (68, 16, U, 0), "_parity3": (84, 6, U, 0), # word4 @@ -357,7 +378,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "_reserved2": (68, 16, U, 0), "_parity3": (84, 6, U, 0), # word4 @@ -392,7 +413,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "alpha0": (68, 8, S, P2_N30), "alpha1": (76, 8, S, P2_N27), "_parity3": (84, 6, U, 0), @@ -436,7 +457,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "sv65asc": (68, 4, U, 0), "sv66asc": (72, 4, U, 0), "sv67asc": (76, 4, U, 0), @@ -503,7 +524,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "avail": (68, 2, U, 0), "erd1": (70, 6, U, 0), "erd2": (76, 6, U, 0), @@ -568,7 +589,7 @@ **GPS_LNAV_HOW, # word3 "dataid": (60, 2, U, 0), - "svcode": (62, 6, U, 0), + "pageid": (62, 6, U, 0), "_reserved2": (68, 16, U, 0), "_parity3": (84, 6, U, 0), # word4 @@ -638,17 +659,332 @@ GPS_LNAV_SUBFRAME_5_P23 = GPS_LNAV_SUBFRAME_5_P01 GPS_LNAV_SUBFRAME_5_P24 = GPS_LNAV_SUBFRAME_5_P01 +# ********************************************************************** +# CNAV - L2C, L5 +# ********************************************************************** -""" -Utility to apply and validate offsets and word partitions -""" -# offset = 0 -# w = 0 -# for key, (_, len, typ, sca) in GPS_LNAV_SUBFRAME_45_GENERIC.items(): -# if not offset % 30: -# w += 1 -# print(f"# word{w}") -# print(f'"{key}": ({offset},{len},{typ},{sca}),') -# offset += len -# if w != 10: -# print("!!! check attribute lengths !!!") +GPS_CNAV_TLM = { + VALPREAMBLE: 0b10001011, # optional, used to validate preamble value + PREAMBLE: (0, 8, U, 0), + "prn": (8, 6, U, 0), + SFR: (14, 6, U, 0), + TOW: (20, 17, U, 0), + "alert": (37, 1, U, 0), +} + +GPS_CNAV_CLOCK = { + "top": (38, 11, U, 300), + "uraned0": (49, 5, S, 0), + "uraned1": (54, 3, U, 0), + "uraned2": (57, 3, U, 0), + TOC: (61, 11, U, 300), + "af0n": (71, 26, S, P2_N35), + "af1n_msb": (97, 3, S, P2_N48), + "af1n_lsb": (100, 17, S, P2_N48), + "af2n": (117, 10, S, P2_N60), +} + +GPS_CNAV_PARITY = { + "_parity": (276, 24, U, 0), +} + +GPS_CNAV_RAP = { + "prn": (0, 6, U, 0), + "deltaa": (6, 8, S, P2_P9), # Relative to Aref = 26,559,710 meters, meters + "omega0": (14, 7, S, P2_N6), # semi-circles + "phi0": (21, 7, S, P2_N6), # M0 + omega, semi-circles + "l1health": (28, 1, U, 0), + "l2health": (29, 1, U, 0), + "l5health": (30, 1, U, 0), +} +# 31 bit reduced almanac packet + +GPS_CNAV_CDC = { + "prn": (0, 8, U, 0), + "deltaaf0": (8, 13, U, P2_N35), + "deltaaf1": (21, 8, U, P2_N51), + "udra": (29, 5, S, 0), +} +# 34 bit clock differential correction + +GPS_CNAV_EDC = { + "prn": (0, 8, U, 0), + "deltaalpha": (8, 14, S, P2_N34), + "deltabeta": (22, 14, S, P2_N34), + "deltalambda": (36, 15, S, P2_N32), + "deltai": (51, 12, S, P2_N31), + "deltaomega": (63, 12, S, P2_N32), + "deltaa": (75, 12, S, P2_N9), + "udradot": (87, 5, S, 0), +} +# 92 bit ephemeris differential correction + +GPS_CNAV_SUBFRAME_10 = { + **GPS_CNAV_TLM, + WN: (38, 13, U, 0), + "l1health": (51, 1, U, 0), + "l2health": (52, 1, U, 0), + "l5health": (53, 1, U, 0), + "top": (54, 11, U, 300), + "uraed": (65, 5, U, 0), + "toe": (70, 11, U, 300), + "deltaa_msb": (81, 19, S, P2_N9), + "deltaa_lsb": (100, 7, S, P2_N9), + "adot": (107, 25, S, P2_N21), + "deltan0": (132, 17, S, P2_N44), + "deltan0dot": (149, 23, S, P2_N57), + "m0_msb": (172, 28, S, P2_N32), + "m0_lsb": (200, 5, S, P2_N32), + "e": (205, 33, U, P2_N34), + "omega": (238, 33, S, P2_N32), + "integrity": (271, 1, U, 0), + "l2phase": (272, 1, U, 0), + "_reserved1": (273, 3, U, 0), + **GPS_CNAV_PARITY, +} # Ephemeris 1 + +GPS_CNAV_SUBFRAME_11 = { + **GPS_CNAV_TLM, + "toe": (38, 11, U, 300), + "omega0": (49, 33, S, P2_N32), + "i0_msb": (82, 18, S, P2_N32), + "i0_lsb": (100, 15, S, P2_N32), + "deltaomegadot": (115, 17, S, P2_N44), + "idot": (132, 15, S, P2_N44), + "cis": (147, 16, S, P2_N30), + "cic": (163, 16, S, P2_N30), + "crs_msb": (179, 21, S, P2_N8), + "crs_lsb": (200, 3, S, P2_N8), + "crc": (203, 24, S, P2_N8), + "cus": (227, 21, S, P2_N30), + "cuc": (248, 21, S, P2_N30), + "_reserved1": (269, 7, U, 0), + **GPS_CNAV_PARITY, +} # Ephemeris 2 + +GPS_CNAV_SUBFRAME_12 = { + **GPS_CNAV_TLM, + "wna": (38, 13, U, 0), + "toa": (51, 8, U, P2_P12), + "rap1": (59, 31, U, 0), # GPS_CNAV_RAP + "rap2_msb": (90, 10, U, 0), + "rap2_lsb": (100, 21, U, 0), + "rap3": (121, 31, U, 0), + "rap4": (152, 31, U, 0), + "rap5_msb": (183, 17, U, 0), + "rap5_lsb": (200, 14, U, 0), + "rap6": (214, 31, U, 0), + "rap7": (245, 31, U, 0), + **GPS_CNAV_PARITY, +} # Reduced Almanac + +GPS_CNAV_SUBFRAME_13 = { + **GPS_CNAV_TLM, + "topd": (38, 11, U, 300), + "tod": (49, 11, U, 300), + "dctype1": (60, 1, U, 0), + "cdc1": (61, 34, U, 0), # GPS_CNAV_CDC + "dctype2": (95, 1, U, 0), + "cdc2_msb": (96, 4, U, 0), + "cdc2_lsb": (100, 30, U, 0), + "dctype3": (130, 1, U, 0), + "cdc3": (131, 34, U, 0), + "dctype4": (165, 1, U, 0), + "cdc4": (166, 34, U, 0), + "dctype5": (200, 1, U, 0), + "cdc5": (201, 34, U, 0), + "dctype6": (235, 1, U, 0), + "cdc6": (236, 34, U, 0), + "_reserved1": (270, 6, U, 0), + **GPS_CNAV_PARITY, +} # Clock Differential Correction + +GPS_CNAV_SUBFRAME_14 = { + **GPS_CNAV_TLM, + "topd": (38, 11, U, 300), + "tod": (49, 11, U, 300), + "dctype1": (60, 1, U, 0), + "edc1_msb": (61, 39, U, 0), # GPS_CNAV_EDC + "edc1_lsb": (100, 53, U, 0), + "dctype2": (153, 1, U, 0), + "edc2_msb": (154, 46, U, 0), + "edc2_lsb": (200, 46, U, 0), + "_reserved1": (226, 30, U, 0), + **GPS_CNAV_PARITY, +} # Ephemeris Differential Correction + +GPS_CNAV_SUBFRAME_15 = { + **GPS_CNAV_TLM, + "text": (38, 232, U, 0), + "textpage": (270, 4, U, 0), + "_reserved1": (274, 2, U, 0), + **GPS_CNAV_PARITY, +} # Text + +GPS_CNAV_SUBFRAME_30 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "tgd": (127, 13, S, P2_N35), + "iscl1ca": (140, 13, S, P2_N35), + "iscl2c": (153, 13, S, P2_N35), + "iscl5i5": (166, 13, S, P2_N35), + "iscl5q5": (179, 13, S, P2_N35), + "alpha0": (192, 8, U, P2_N30), # where is scaling defined for CNAV? + "alpha1": (200, 8, U, P2_N27), # have assumed same as LNAV + "alpha2": (208, 8, U, P2_N24), + "alpha3": (216, 8, U, P2_N24), + "beta0": (224, 8, U, P2_P11), + "beta1": (232, 8, U, P2_P14), + "beta2": (240, 8, U, P2_P16), + "beta3": (248, 8, U, P2_P16), + "wno": (256, 8, U, 0), + "_reserved1": (264, 12, U, 0), + **GPS_CNAV_PARITY, +} # Clock, IONO & Group Delay + +GPS_CNAV_SUBFRAME_31 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "wna": (127, 13, U, 0), + "toa": (140, 8, U, P2_P12), + "rap1": (148, 31, U, 0), # GPS_CNAV_RAP + "rap2_msb": (179, 21, U, 0), + "rap2_lsb": (200, 10, U, 0), + "rap3": (210, 31, U, 0), + "rap4": (241, 31, U, 0), + "_reserved1": (272, 4, U, 0), + **GPS_CNAV_PARITY, +} # Clock & Reduced Almanac + +GPS_CNAV_SUBFRAME_32 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "teop": (127, 16, U, P2_P4), + "pmx": (143, 21, S, P2_N20), + "pmxdot": (164, 15, S, P2_N21), + "pmy": (179, 21, S, P2_N20), + "pmydot": (200, 15, S, P2_N21), + "deltautgps": (215, 31, S, P2_N23), + "deltautgpsdot": (246, 19, S, P2_N25), + "_reserved1": (265, 11, U, 0), + **GPS_CNAV_PARITY, +} # Clock & EOP + +GPS_CNAV_SUBFRAME_33 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "a0": (127, 16, S, P2_N35), + "a1": (43, 13, S, P2_N51), + "a2": (156, 7, S, P2_N68), + "deltatls": (163, 8, S, 1), + "tot": (171, 16, U, P2_P4), + "wnot": (187, 13, U, 1), + "wnlsf": (200, 13, U, 1), + "dn": (213, 4, U, 1), + "deltatlsf": (217, 8, S, 1), + "_reserved1": (225, 52, U, 0), + **GPS_CNAV_PARITY, +} # Clock & UTC + +GPS_CNAV_SUBFRAME_34 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "topd": (127, 11, U, 300), + "tod": (138, 11, U, 300), + "dctype": (149, 1, U, 0), + "cdc": (150, 34, U, 0), # GPS_CNAV_CDC + "edc_msb": (184, 16, U, 0), # GPS_CNAV_EDC + "edc_lsb": (200, 76, U, 0), + "_reserved1": (225, 52, U, 0), + **GPS_CNAV_PARITY, +} # Clock & Differential Correction + +GPS_CNAV_SUBFRAME_35 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "tggto": (127, 16, U, P2_P4), + "wnggto": (143, 13, U, 1), + "gnssid": (156, 3, U, 0), + "a0ggto": (159, 16, U, P2_N35), + "a1ggto": (175, 13, U, P2_N51), + "a2ggto": (188, 7, U, P2_N68), + "_reserved1": (195, 5, U, 0), + "_reserved2": (200, 76, U, 0), + **GPS_CNAV_PARITY, +} # Clock & GGTO (GPS/GNSS Time Offset) + +GPS_CNAV_SUBFRAME_36 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "text_msb": (127, 73, U, 0), + "text_msb": (200, 71, U, 0), + "textpage": (271, 4, U, 0), + "_reserved1": (275, 1, U, 0), + **GPS_CNAV_PARITY, +} # Clock & Text + +GPS_CNAV_SUBFRAME_37 = { + **GPS_CNAV_TLM, + **GPS_CNAV_CLOCK, + "wna": (127, 13, U, 0), + "toa": (140, 8, U, P2_P12), + "prna": (148, 6, U, 0), + "l1health": (154, 1, U, 0), + "l2health": (155, 1, U, 0), + "l5health": (156, 1, U, 0), + "e": (157, 11, U, P2_N16), + "omegai": (168, 11, S, P2_N14), + "omegadot": (179, 11, S, P2_N33), + "sqrta_msb": (190, 10, U, P2_N4), + "sqrta_lsb": (200, 7, U, P2_N4), + "omega0": (107, 16, S, P2_N15), + "omega": (223, 16, S, P2_N15), + "m0": (239, 16, S, P2_N15), + "af0": (255, 11, S, P2_N20), + "af1": (266, 10, S, P2_N37), + **GPS_CNAV_PARITY, +} # Clock & Midi Almanac + +GPS_CNAV_SUBFRAME_40 = { + **GPS_CNAV_TLM, + "gnssid": (38, 4, U, 0), + "wnism": (42, 13, U, 1), + "towism": (55, 6, U, 4), + "tcorrel": (61, 4, U, 0), + "bnom": (65, 4, U, 0), + "lambdanom": (69, 4, U, 0), + "rsat": (73, 4, U, 0), + "pconst": (77, 4, U, 0), + "mfd": (81, 4, U, 0), + "servicelevel": (85, 3, U, 0), + "mask": (88, 63, U, 0), + "filler": (151, 94, U, 0), + "_ism_crc": (244, 32, U, 0), + **GPS_CNAV_PARITY, +} # Integrity Support Message + +# mapping for subframe acquisition mask sfracq +GPS_SFRACQ_MAP = { + # LNAV + 1: 1, + 2: 2, + 3: 4, + 4: 8, + 5: 16, + # CNAV + 10: 1, + 11: 2, + 30: 4, + 12: 8, + 13: 16, + 14: 32, + 15: 64, + 31: 128, + 32: 256, + 33: 512, + 34: 1024, + 35: 2048, + 36: 4096, + 37: 8192, + 40: 16384, +} diff --git a/tests/test_rawnav.py b/tests/test_rawnav.py index f3d2282..980f777 100644 --- a/tests/test_rawnav.py +++ b/tests/test_rawnav.py @@ -17,6 +17,7 @@ GPS_LNAV_SUBFRAME_1, GPS_LNAV_SUBFRAME_3, GPS_LNAV_SUBFRAME_2, + GPS_SFRACQ_MAP, ) SUBFRAME1 = { @@ -68,7 +69,7 @@ def testconstructor(self): epoch=datetime(2026, 4, 18, 9, 42, 13, tzinfo=timezone.utc), ) data = int("10001011" + "1111" + "10100101" * 36, 2) - raw.parse(data, sfr, sequence=False) + raw.parse(data, sfr, sfrmap=GPS_SFRACQ_MAP, sequence=False) print(f'"{raw}",') #self.assertEqual(str(raw), EXPECTED_RESULTS[i]) #self.assertEqual((raw.gnss, raw.svid, raw.sigid, raw.epoch),("G",32,"1C",datetime(2026, 4, 18, 9, 42, 13, tzinfo=timezone.utc))) diff --git a/tests/test_rinex.py b/tests/test_rinex.py index 1010f1a..fafcf35 100644 --- a/tests/test_rinex.py +++ b/tests/test_rinex.py @@ -21,6 +21,7 @@ from pygnssutils.rinex_helpers import ( DRNX, FRNX, + get_epoch, adjust_time_units, format_antennabsight, format_antennadeltahen, @@ -39,6 +40,7 @@ format_glonassphasebias, format_headerend, format_interval, + format_eop, format_iono_corr, format_leapseconds, format_marker, @@ -58,7 +60,10 @@ format_time_corr, format_timefirstlast, format_version, + format_nav_typesvmssg, listify, + format_sto, + format_ion, ) SENSORTYPES = { @@ -178,6 +183,31 @@ def testformat_filename(self): self.assertTrue( str(res), "/Users/steve/Downloads/pygpsdata_R_999912310000_00U_00U_GN.rnx" ) + firstobs = datetime(2026, 3, 14, 12, 4, 6) + lastobs = firstobs + timedelta(minutes=60) + interval=15 + res = format_filename( + rinextype="O", + gnssfilter=[GPS], + startepoch=firstobs, + endepoch=lastobs, + interval=interval, + outputpath=Path("/Users/steve/Downloads"), + form="IGS", + site="SITE", + marker=1, + receiver=1, + country="GBR", + ) + # print(res) + self.assertTrue( + str(res), "/Users/steve/Downloads/SITE11GBR_R_202603141204_60M_15S_GO.rnx" + ) + + def test_getepoch(self): + res = get_epoch(366,411634,"G") + # print(res) + self.assertEqual(res, (datetime(2026, 4, 16, 18, 20, 16, tzinfo=timezone.utc), 2414)) def testformat_antennabsight(self): res = format_antennabsight() @@ -294,44 +324,114 @@ def testformat_iono_corr(self): "GPSA 1.2481e-07 5.0391e-06 2.3771e-07 1.2346e-13 A 02 IONOSPHERIC CORR\n" "GPSB 1.2481e-07 5.0391e-06 2.3771e-07 1.2346e-13 B 14 IONOSPHERIC CORR\n" ) - ionocorr = { - "GPSA": { - "parm1": 0.00000012481234567890, - "parm2": 0.0000050391234567890, - "parm3": 0.00000023771234567890, - "parm4": 0.00000000000012345678909, - "timemark": "A", - "svid": 2, - }, - "GPSB": { - "parm1": 0.00000012481234567890, - "parm2": 0.0000050391234567890, - "parm3": 0.00000023771234567890, - "parm4": 0.00000000000012345678909, - "timemark": "B", - "svid": 14, - }, - } - res = format_iono_corr(ionocorr) + res = format_iono_corr( + svid=2, + timemark="A", + corrtype="GPSA", + parm1=0.00000012481234567890, + parm2=0.0000050391234567890, + parm3=0.00000023771234567890, + parm4=0.00000000000012345678909, + ) + res += format_iono_corr( + svid=14, + timemark="B", + corrtype="GPSB", + parm1=0.00000012481234567890, + parm2=0.0000050391234567890, + parm3=0.00000023771234567890, + parm4=0.00000000000012345678909, + ) # print(res) self.assertEqual(res, EXPECTED_RESULT) + def testformat_ion(self): + EXPECTED_RESULT = ( + "> ION G24 LNAV XXXX\n" + " 2026 05 13 08 34 02 1.234567000000e-12 1.234567000000e-12-1.234567000000e-12\n" + " 1.234567000000e-12 1.234567000000e-12 1.234567000000e-12-1.234567000000e-12\n" + " 1.234567000000e-12\n" + ) + res = format_ion( + svcode="G24", + msgtype="LNAV", + msgsubtype="XXXX", + epoch=datetime(2026,5,13,8,34,2,tzinfo=timezone.utc), + a0=1.234567e-12, + a1=1.234567e-12, + a2=-1.234567e-12, + a3=1.234567e-12, + b0=1.234567e-12, + b1=1.234567e-12, + b2=-1.234567e-12, + b3=1.234567e-12, + ) + #print(res) + self.assertEqual(res, EXPECTED_RESULT) + + def testformat_eop(self): + EXPECTED_RESULT = ( + "> EOP G24 LNAV XXXX\n" + " 2026 05 13 08 34 02 2.082471847534e-01-6.551742553711e-04 0.000000000000e+00\n" + " 3.444433212280e-01-9.121894836426e-04 0.000000000000e+00\n" + " 1.729860000000e+05-1.754972934723e-01 5.635917186737e-04 0.000000000000e+00\n" + ) + res = format_eop( + svcode="G24", + msgtype="LNAV", + msgsubtype="XXXX", + epoch=datetime(2026,5,13,8,34,2,tzinfo=timezone.utc), + tom=1.729860000000e+05, + xp=2.082471847534e-01, + dxpdt=-6.551742553711e-04, + dxpdt2=0.000000000000e+00, + yp=3.444433212280e-01, + dypdt=-9.121894836426e-04, + dypdt2=0.000000000000e+00, + deltaut1=-1.754972934723e-01, + ddeltaut1dt=5.635917186737e-04, + d2deltaut1dt2=0.000000000000e+00, + ) + print(res) + self.assertEqual(res, EXPECTED_RESULT) + def testformat_time_corr(self): EXPECTED_RESULT = "GAUT 3.7252902980e-09 5.329070520e-15 345600 1849 E10 5 TIME SYSTEM CORR\n" - timecorr = { - "GAUT": { - "a0": 0.000000003725290298, - "a1": 0.00000000000000532907052, - "timeref": 345600, - "weekno": 1849, - "svcode": "E10", - "source": "5", - } - } - res = format_time_corr(timecorr) + res = format_time_corr( + svcode="E10", + corrtype="GAUT", + timeref=345600, + weekno=1849, + source="5", + a0=0.000000003725290298, + a1=0.00000000000000532907052 + ) # print(res) self.assertEqual(res, EXPECTED_RESULT) + + def testformat_sto(self): + EXPECTED_RESULT = ( + "> STO G24 LNAV XXXX\n" + " 2026 05 13 08 34 02 GPUT SSSS UTC(USNO) \n" + " 4.567890000000e+05 1.234567000000e-23 1.234567000000e-23-1.234567800000e-12\n" + ) + res = format_sto( + svcode="G24", + msgtype="LNAV", + msgsubtype="XXXX", + epoch=datetime(2026,5,13,8,34,2,tzinfo=timezone.utc), + timecode="GPUT", + sbasid="SSSS", + utcid="UTC(USNO)", + tot=456789, + a0=1.234567e-23, + a1=1.234567e-23, + a2=-1.2345678e-12 + ) + # print(res) + self.assertEqual(res,EXPECTED_RESULT) + def testformat_met_obstypes(self): EXPECTED_RESULT = ( " 10 PR TD HR ZW ZD ZT WD WS PI# / TYPES OF OBSERV\n" @@ -364,6 +464,12 @@ def testformat_met_sensorpos(self): # print(res) self.assertEqual(res, EXPECTED_RESULT) + def testformat_nav_typesvmssg(self): + EXPECTED_RESULT = "> EPH G04 LNAV \n" + res = format_nav_typesvmssg("EPH", "G04", "LNAV") + # print(res) + self.assertEqual(res, EXPECTED_RESULT) + def testadjust_time_units(self): self.assertEqual(adjust_time_units(34), (34, "S")) self.assertEqual(adjust_time_units(345), (6.0, "M"))