diff --git a/.gitignore b/.gitignore index ffda4974..6d48f78f 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,5 @@ __pycache__ **/tmp_test*/* **.profile **.pstat +src/prx/sp3/sp3_files/* +src/prx/precise_corrections/sp3/sp3_files/2023/001/IGS0OPSFIN_20230010000_01D_15M_ORB.SP3 diff --git a/src/prx/precise_corrections/antex/antex_file_discovery.py b/src/prx/precise_corrections/antex/antex_file_discovery.py new file mode 100644 index 00000000..efd1e8b1 --- /dev/null +++ b/src/prx/precise_corrections/antex/antex_file_discovery.py @@ -0,0 +1,162 @@ +from ftplib import FTP +import ftplib +from pathlib import Path +import re + +import georinex +import pandas as pd +import urllib +import fnmatch + +from prx import converters +from prx import util +from prx.util import timestamp_to_mid_day + +""" +Logic overview: +- The code first extracts the GPS week from the provided RINEX observation file. +- It then identifies the most recent ANTEX file available both locally and remotely. +- If both local and remote ANTEX files are the same, the local version is reused, as the local database is considered up to date. +- Otherwise, the latest remote file is downloaded and used. +- Note: If the most recent ANTEX file has a GPS week that is earlier than the observation file's GPS week +(e.g., the remote database has not yet been updated), a warning is displayed to indicate that the most recent ANTEX file +was still used despite being older than the observation file. +""" + +log = util.get_logger(__name__) + +atx_filename = "igs20_????.atx" + + +def date_to_gps_week(date: pd.Timestamp): + """ + Convert a Timestamp to GPS week + + GPS week starts on 01/06/1980 (sunday) + """ + gps_epoch = pd.Timestamp("1980-01-06T00:00:00Z") + delta = date.tz_localize("UTC") - gps_epoch + gps_week = delta.days // 7 + return gps_week + + +def extract_gps_week(filename: str) -> int: + """ + Extract GPS week from filenames like 'igsxx_WWWW.atx'. + """ + match = re.search(r"igs...(\d{4})\.atx", filename) + if match: + return int(match.group(1)) + return -1 + + +def atx_file_database_folder(): + """ + Returns the path to the folder where ATX database files are stored. + """ + db_folder = ( + util.prx_repository_root() / "src/prx/precise_corrections/antex/atx_files" + ) + db_folder.mkdir(exist_ok=True) + return db_folder + + +def find_latest_local_antex_file(db_folder=atx_file_database_folder()): + candidates = list(db_folder.glob(atx_filename)) + if not candidates: + return None + return max(candidates, key=lambda c: extract_gps_week(c)) + + +def list_ftp_directory(server, folder): + ftp = FTP(server) + ftp.login() + ftp.cwd(folder) + dir_list = [] + ftp.dir(dir_list.append) + return [c.split()[-1].strip() for c in dir_list] + + +def fetch_latest_remote_antex_file(): + """ + List the ANTEX files available online and returns the latest + """ + server = "gssc.esa.int" + remote_folder = "/igs/station/general" + candidates = list_ftp_directory(server, remote_folder) + candidates = [c for c in candidates if fnmatch.fnmatch(c, atx_filename)] + if not candidates: + return None + return max(candidates, key=lambda c: extract_gps_week(c)) + + +def check_online_availability(file: str, folder: Path) -> Path | None: + """ + Need to keep the same inputs as try_downloading_sp3_ftp, in order to be able to use `unittest.mock.patch` in tests + """ + server = "gssc.esa.int" + remote_folder = "/igs/station/general" + ftp = FTP(server) + ftp.login() + ftp.cwd(remote_folder) + try: + ftp.size(file) + return folder.joinpath(file) + except ftplib.error_perm: + log.warning(f"{file} not available on {server}") + return None + + +def try_downloading_atx_ftp(file: str, folder: Path): + """ + Download the wanted remote file + """ + server = "gssc.esa.int" + remote_folder = "/igs/station/general" + ftp_file = f"ftp://{server}/{remote_folder}/{file}" + local_file = folder / file + urllib.request.urlretrieve(ftp_file, local_file) + if not local_file.exists(): + log.warning(f"Could not download {ftp_file}") + return None + log.info(f"Downloaded ANTEX file {ftp_file}") + return local_file + + +def get_atx_file(date: pd.Timestamp, db_folder=atx_file_database_folder()): + gps_week = date_to_gps_week(date) + latest_atx_local = find_latest_local_antex_file(db_folder) + latest_atx_remote = fetch_latest_remote_antex_file() + if latest_atx_remote == latest_atx_local: + return latest_atx_local + elif latest_atx_remote and ( + not latest_atx_local + or extract_gps_week(latest_atx_remote) > extract_gps_week(latest_atx_local) + ): + # Download the latest file + atx_file = try_downloading_atx_ftp(latest_atx_remote, db_folder) + if atx_file is not None: + if gps_week > extract_gps_week(atx_file.name): + log.warning( + f"No ANTEX file found for the target GPS week {gps_week} — using the most recent available instead." + ) + return atx_file + + raise FileNotFoundError("No file ANTEX found locally or online.") + + +def discover_or_download_atx_file(observation_file_path=Path()): + """ + Returns the path to a valid antes file (local or downloaded) corresponding to the observation file. + """ + log.info(f"Finding auxiliary files for {observation_file_path} ...") + rinex_3_obs_file = converters.anything_to_rinex_3(observation_file_path) + header = georinex.rinexheader(rinex_3_obs_file) + + t_start = timestamp_to_mid_day( + util.rinex_header_time_string_2_timestamp_ns(header["TIME OF FIRST OBS"]) + - pd.Timedelta(200, unit="milliseconds") + ) + + atx_file = get_atx_file(t_start) + return atx_file diff --git a/src/prx/precise_corrections/antex/test/datasets/test_antex_file_discovery.py b/src/prx/precise_corrections/antex/test/datasets/test_antex_file_discovery.py new file mode 100644 index 00000000..9535a0da --- /dev/null +++ b/src/prx/precise_corrections/antex/test/datasets/test_antex_file_discovery.py @@ -0,0 +1,139 @@ +import os +from pathlib import Path +import shutil +import georinex +import pandas as pd +import pytest +from unittest.mock import patch + +import prx +from prx import util +from prx.precise_corrections.antex import antex_file_discovery as atx + + +@pytest.fixture +def set_up_test(): + test_directory = Path(f"./tmp_test_directory_{__name__}").resolve() + if test_directory.exists(): + # Make sure the expected files has not been generated before and is still on disk due to e.g. a previous + # test run having crashed: + shutil.rmtree(test_directory) + os.makedirs(test_directory) + test_obs_file = test_directory.joinpath("TLSE00FRA_R_20230010100_10S_01S_MO.crx.gz") + # local_database = atx.atx_file_folder(pd.Timestamp('2023-01-01'), Path('src/prx/precise_corrections/atx/test/datasets')) + + shutil.copy( + util.prx_repository_root() + / f"src/prx/test/datasets/TLSE_2023001/{test_obs_file.name}", + test_obs_file, + ) + assert test_obs_file.exists() + + yield {"test_obs_file": test_obs_file} + shutil.rmtree(test_directory) + + +def test_extract_gps_week(): + filenames = [ + "igs20_2134.atx", + "igs_10.atx", + "igs14_abc.atx", + ] + expected_returns = [2134, -1, -1] + gps_week = [atx.extract_gps_week(f) for f in filenames] + assert gps_week == expected_returns + + +def test_download_if_not_local(set_up_test): + """ + Tests that the ANTEX file is downloaded when no local file is available, + but a valid remote file exists. + + Scenario : + - No local ANTEX file available + - A remote ANTEX file exists + -> The function should trigger the download of the remote file. + """ + obs_file = set_up_test["test_obs_file"] + header = georinex.rinexheader(obs_file) + t_start = util.rinex_header_time_string_2_timestamp_ns( + header["TIME OF FIRST OBS"] + ) - pd.Timedelta(200, unit="milliseconds") + + db_folder = set_up_test["test_obs_file"].parent + downloaded_atx = atx.get_atx_file(t_start, db_folder) + + # Ensure the file was downloaded and exists + assert downloaded_atx is not None + assert downloaded_atx.exists() + + +def test_download_if_remote_is_newer(set_up_test): + """ + Tests that the ANTEX file is downloaded when the remote file has a newer GPS week + than the local one. + + Scenario : + - local file = igs20_2370.atx + - remote file = igs20_2375.atx + -> The function should download igs20_2375.atx + """ + latest_local = "igs20_2370.atx" + latest_remote = "igs20_2375.atx" + date = pd.Timestamp("2025-01-01 12:00:00") + db_folder = set_up_test["test_obs_file"].parent + + with ( + patch( # replace download function by online availability check + "prx.precise_corrections.antex.antex_file_discovery.try_downloading_atx_ftp", + new=prx.precise_corrections.antex.antex_file_discovery.check_online_availability, + ), + patch( # simulate that the latest local file found is 'igs20_2370.atx' + "prx.precise_corrections.antex.antex_file_discovery.find_latest_local_antex_file", + return_value=latest_local, + ), + patch( # simulate that the latest remote file found is 'igs20_2375.atx' + "prx.precise_corrections.antex.antex_file_discovery.fetch_latest_remote_antex_file", + return_value=latest_remote, + ), + ): + result = atx.get_atx_file(date, db_folder) + assert result is not None + assert result.name == latest_remote + + +def test_skip_download_if_same_week(set_up_test): + """ + Tests that the ANTEX file is not downloaded when the remote and local files + are from the same GPS week. + + Scenario: + - local file = igs20_2375.atx + - remote file = igs20_2375.atx + → The function should skip download and return the local file + """ + latest_local = "igs20_2375.atx" + latest_remote = "igs20_2375.atx" + date = pd.Timestamp("2025-01-01 12:00:00") + db_folder = set_up_test["test_obs_file"].parent + + with ( + patch( # mock download function + "prx.precise_corrections.antex.antex_file_discovery.try_downloading_atx_ftp" + ) as mock_download, + patch( # simulate that the latest local file found is 'igs20_2375.atx' + "prx.precise_corrections.antex.antex_file_discovery.find_latest_local_antex_file", + return_value=latest_local, + ), + patch( # simulate that the latest remote file found is 'igs20_2375.atx' + "prx.precise_corrections.antex.antex_file_discovery.fetch_latest_remote_antex_file", + return_value=latest_remote, + ), + ): + result = atx.get_atx_file(date, db_folder) + + assert result is not None + # Ensure that the download function was not called + mock_download.assert_not_called() + # Ensure the returned file is the local one + assert result == latest_local diff --git a/src/prx/sp3/README.md b/src/prx/precise_corrections/sp3/README.md similarity index 100% rename from src/prx/sp3/README.md rename to src/prx/precise_corrections/sp3/README.md diff --git a/src/prx/sp3/__init__.py b/src/prx/precise_corrections/sp3/__init__.py similarity index 100% rename from src/prx/sp3/__init__.py rename to src/prx/precise_corrections/sp3/__init__.py diff --git a/src/prx/sp3/evaluate.py b/src/prx/precise_corrections/sp3/evaluate.py similarity index 100% rename from src/prx/sp3/evaluate.py rename to src/prx/precise_corrections/sp3/evaluate.py diff --git a/src/prx/precise_corrections/sp3/sp3_file_discovery.py b/src/prx/precise_corrections/sp3/sp3_file_discovery.py new file mode 100644 index 00000000..eb42bd83 --- /dev/null +++ b/src/prx/precise_corrections/sp3/sp3_file_discovery.py @@ -0,0 +1,191 @@ +# This module aims at create a local database of SP3 ORB and NAV files. +# Upon request of a particular date, # the availability of the files in the local database will be checked, +# and if missing, they will be downloaded from the IGS FTP servers. +# A priority list is defined to provide a preference order of IGS product, in terms of types (final or rapid products) +# and IGS analysis center. + + +import logging +import ftplib +import os +from pathlib import Path +from typing import List, Tuple + +import georinex +import pandas as pd +import urllib + +from prx import converters, util +from prx.util import timestamp_to_mid_day, timestamp_to_gps_week_and_dow + +log = logging.getLogger(__name__) + +### Since gps week 2238 +priority = [ + ("COD", "FIN"), + ("GRG", "FIN"), + ("GFZ", "FIN"), + ("ESA", "FIN"), + ("WUM", "FIN"), + ("JAX", "FIN"), + ("JPL", "FIN"), + ("MIT", "FIN"), + ("COD", "RAP"), + ("GRG", "RAP"), + ("GFZ", "RAP"), + ("ESA", "RAP"), + ("WUM", "RAP"), +] +# WWWW/AAA0PPPTYP_YYYYDDDHHMM_LEN_SMP_CNT.FMT.gz +# PPP : MGX +# CNT : ORB +# FMT : SP3 + +# This priority list is applied starting from GPS week 2238, to select SP3 orbit and CLK files based on the preferred +# analysis centers and product types. Final ("FIN") products are prioritized due to their higher accuracy and reliability. +# When final products are unavailable, rapid ("RAP") products serve as a fallback. +# Among the analysis centers, COD, GRG, and GFZ are placed at the top based on their long-standing reputation for delivering +# precise and complete orbit solutions within the IGS community. ESA, WUM, JAX, JPL, and MIT follow, as they also provide +# high-quality products, but may differ slightly in availability, latency, or consistency. +# Before GPS week 2238, the same type of SP3 files can be found, but they are stored in /{gps_week}/mgex directories, +# requiring a different file discovery logic. + + +def get_index_of_priority_from_filename(filename: str) -> int: + for i, p in enumerate(priority): + if p[0] in filename and p[1] in filename: + return i + + +def build_sp3_filename(date: pd.Timestamp, aaa_typ: (str, str)) -> (str, str): + # aaa_typ : tuple of str (aaa, typ) + # aaa: IGS analysis center + # typ: IGS product type (RAP or FIN) + yyyy = date.year + ddd = f"{date.day_of_year:03d}" + aaa = aaa_typ[0] + typ = aaa_typ[1] + sp3_filename = f"{aaa}0MGX{typ}_{yyyy}{ddd}0000_01D_05M_ORB.SP3.gz" + clk_filename = f"{aaa}0MGX{typ}_{yyyy}{ddd}0000_01D_30S_CLK.CLK.gz" + return sp3_filename, clk_filename + + +def sp3_file_database_folder() -> Path: + """ + Returns the path to the folder where SP3 database files are stored. + """ + db_folder = util.prx_repository_root() / "src/prx/precise_corrections/sp3/sp3_files" + db_folder.mkdir(exist_ok=True) + return db_folder + + +def sp3_file_folder( + date: pd.Timestamp, parent_folder: Path = sp3_file_database_folder() +) -> Path: + """ + Returns the path to the folder where SP3 files for a specific day are stored. + """ + folder = parent_folder / f"{date.year}/{date.day_of_year:03d}" + folder.mkdir(parents=True, exist_ok=True) + return folder + + +def get_local_sp3( + date: pd.Timestamp, file: str, db_folder=sp3_file_database_folder() +) -> Path | None: + candidates = list(sp3_file_folder(date, db_folder).glob(file)) + if len(candidates) == 0: + return None + log.info(f"Found the sp3 local file : {candidates[0]}") + local_file = converters.compressed_to_uncompressed(candidates[0]) + return local_file + + +def check_online_availability(gps_week: int, folder: Path, file: str) -> Path | None: + """ + Need to keep the same inputs as try_downloading_sp3_ftp, in order to be able to use `unittest.mock.patch` in tests + """ + server = "gssc.esa.int" + if gps_week > 2237: + remote_folder = f"/gnss/products/{gps_week}" + else: + remote_folder = f"/gnss/products/{gps_week}/mgex" + ftp = ftplib.FTP(server) + ftp.login() + ftp.cwd(remote_folder) + try: + ftp.size(file) + return folder.joinpath(Path(file).stem) + except ftplib.error_perm: + log.warning(f"{file} not available on {server}") + return None + + +def try_downloading_sp3_ftp(gps_week: int, folder: Path, file: str) -> Path | None: + server = "gssc.esa.int" + if gps_week > 2237: + remote_folder = f"/gnss/products/{gps_week}" + else: + remote_folder = f"/gnss/products/{gps_week}/mgex" + ftp_file = f"ftp://{server}/{remote_folder}/{file}" + local_compressed_file = folder / file + urllib.request.urlretrieve(ftp_file, local_compressed_file) + if not local_compressed_file.exists(): + log.warning(f"Could not download {ftp_file}") + return None + local_file = converters.compressed_to_uncompressed(local_compressed_file) + os.remove(local_compressed_file) + log.info(f"Downloaded sp3 file {ftp_file}") + return local_file + + +def get_sp3_files( + mid_day_start: pd.Timestamp, + mid_day_end: pd.Timestamp, + db_folder=sp3_file_database_folder(), +) -> List[Tuple[Path]]: + sp3_files = [] + date = mid_day_start + gps_week, _ = timestamp_to_gps_week_and_dow(date) + while date <= mid_day_end: + for p in priority: + sp3_filename, clk_filename = build_sp3_filename(date, p) + file_orb = get_local_sp3(date, sp3_filename, db_folder) + file_clk = get_local_sp3(date, clk_filename, db_folder) + if file_orb is None: + file_orb = try_downloading_sp3_ftp( + gps_week, sp3_file_folder(date, db_folder), sp3_filename + ) + if file_clk is None: + file_clk = try_downloading_sp3_ftp( + gps_week, sp3_file_folder(date, db_folder), clk_filename + ) + if file_orb is not None and file_clk is not None: + sp3_files.append((file_orb, file_clk)) + break + # If we reach the end of the priority list without success + if file_orb is None and file_clk is None and p == priority[-1]: + sp3_files.append((None, None)) + date += pd.Timedelta(1, unit="days") + return sp3_files + + +def discover_or_download_sp3_file(observation_file_path=Path) -> List[Tuple[Path]]: + """ + Returns the path to a valid SP3 file (local or downloaded) corresponding to the observation file. + Tries to respect a priority hierarchy: IGS FIN > COD FIN > GRG FIN > ... > IGS ULR. + """ + log.info(f"Finding sp3 files for {observation_file_path} ...") + rinex_3_obs_file = converters.anything_to_rinex_3(observation_file_path) + header = georinex.rinexheader(rinex_3_obs_file) + + t_start = timestamp_to_mid_day( + util.rinex_header_time_string_2_timestamp_ns(header["TIME OF FIRST OBS"]) + - pd.Timedelta(200, unit="milliseconds") + ) + t_end = timestamp_to_mid_day( + util.rinex_header_time_string_2_timestamp_ns(header["TIME OF LAST OBS"]) + ) + + sp3_files = get_sp3_files(t_start, t_end) + return sp3_files diff --git a/src/prx/sp3/test/__init__.py b/src/prx/precise_corrections/sp3/test/__init__.py similarity index 100% rename from src/prx/sp3/test/__init__.py rename to src/prx/precise_corrections/sp3/test/__init__.py diff --git a/src/prx/precise_corrections/sp3/test/datasets/2023/001/COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz b/src/prx/precise_corrections/sp3/test/datasets/2023/001/COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz new file mode 100644 index 00000000..c417cba3 Binary files /dev/null and b/src/prx/precise_corrections/sp3/test/datasets/2023/001/COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz differ diff --git a/src/prx/precise_corrections/sp3/test/datasets/2023/001/GFZ0MGXRAP_20230010000_01D_05M_ORB.SP3.gz b/src/prx/precise_corrections/sp3/test/datasets/2023/001/GFZ0MGXRAP_20230010000_01D_05M_ORB.SP3.gz new file mode 100644 index 00000000..a5524f6b Binary files /dev/null and b/src/prx/precise_corrections/sp3/test/datasets/2023/001/GFZ0MGXRAP_20230010000_01D_05M_ORB.SP3.gz differ diff --git a/src/prx/precise_corrections/sp3/test/datasets/2023/001/JAX0MGXFIN_20230010000_01D_30S_CLK.CLK.gz b/src/prx/precise_corrections/sp3/test/datasets/2023/001/JAX0MGXFIN_20230010000_01D_30S_CLK.CLK.gz new file mode 100644 index 00000000..8f087ef2 Binary files /dev/null and b/src/prx/precise_corrections/sp3/test/datasets/2023/001/JAX0MGXFIN_20230010000_01D_30S_CLK.CLK.gz differ diff --git a/src/prx/precise_corrections/sp3/test/datasets/2023/001/WUM0MGXFIN_20230010000_01D_05M_ORB.SP3.gz b/src/prx/precise_corrections/sp3/test/datasets/2023/001/WUM0MGXFIN_20230010000_01D_05M_ORB.SP3.gz new file mode 100644 index 00000000..4efc65b3 Binary files /dev/null and b/src/prx/precise_corrections/sp3/test/datasets/2023/001/WUM0MGXFIN_20230010000_01D_05M_ORB.SP3.gz differ diff --git a/src/prx/precise_corrections/sp3/test/datasets/2023/001/WUM0MGXFIN_20230010000_01D_30S_CLK.CLK.gz b/src/prx/precise_corrections/sp3/test/datasets/2023/001/WUM0MGXFIN_20230010000_01D_30S_CLK.CLK.gz new file mode 100644 index 00000000..c218e4f9 Binary files /dev/null and b/src/prx/precise_corrections/sp3/test/datasets/2023/001/WUM0MGXFIN_20230010000_01D_30S_CLK.CLK.gz differ diff --git a/src/prx/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB.SP3 b/src/prx/precise_corrections/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB.SP3 similarity index 100% rename from src/prx/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB.SP3 rename to src/prx/precise_corrections/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB.SP3 diff --git a/src/prx/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB_one_sample_removed.SP3 b/src/prx/precise_corrections/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB_one_sample_removed.SP3 similarity index 100% rename from src/prx/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB_one_sample_removed.SP3 rename to src/prx/precise_corrections/sp3/test/datasets/WUM0MGXULT_20220010000_01H_05M_ORB_one_sample_removed.SP3 diff --git a/src/prx/sp3/test/test_evaluate.py b/src/prx/precise_corrections/sp3/test/test_evaluate.py similarity index 98% rename from src/prx/sp3/test/test_evaluate.py rename to src/prx/precise_corrections/sp3/test/test_evaluate.py index 67c70764..c8192458 100644 --- a/src/prx/sp3/test/test_evaluate.py +++ b/src/prx/precise_corrections/sp3/test/test_evaluate.py @@ -1,7 +1,7 @@ import pandas as pd import numpy as np from pathlib import Path -from prx.sp3.evaluate import compute +from prx.precise_corrections.sp3.evaluate import compute import shutil import pytest import os diff --git a/src/prx/precise_corrections/sp3/test/test_sp3_file_discovey.py b/src/prx/precise_corrections/sp3/test/test_sp3_file_discovey.py new file mode 100644 index 00000000..de732790 --- /dev/null +++ b/src/prx/precise_corrections/sp3/test/test_sp3_file_discovey.py @@ -0,0 +1,242 @@ +import os +from pathlib import Path +import shutil +import georinex +import pytest + +import pandas as pd +from unittest.mock import patch +import prx +import prx.util +from prx import util +from prx.precise_corrections.sp3 import sp3_file_discovery as sp3 + + +@pytest.fixture +def set_up_test(): + test_directory = Path(f"./tmp_test_directory_{__name__}").resolve() + if test_directory.exists(): + # Make sure the expected files has not been generated before and is still on disk due to e.g. a previous + # test run having crashed: + shutil.rmtree(test_directory) + os.makedirs(test_directory) + test_obs_file = test_directory.joinpath("TLSE00FRA_R_20230010100_10S_01S_MO.crx.gz") + + shutil.copy( + util.prx_repository_root() + / f"src/prx/test/datasets/TLSE_2023001/{test_obs_file.name}", + test_obs_file, + ) + assert test_obs_file.exists() + test_date = pd.Timestamp("2023-01-01") + sp3_subfolder = sp3.sp3_file_folder(test_date, test_directory) + + # copy all sp3 files + for file in ( + util.prx_repository_root() + / "src" + / "prx" + / "precise_corrections" + / "sp3" + / "test" + / "datasets" + / "2023" + / "001" + ).glob("*"): + test_sp3_file = sp3_subfolder / file.name + shutil.copy( + util.prx_repository_root() + / f"src/prx/precise_corrections/sp3/test/datasets/2023/001/{test_sp3_file.name}", + test_sp3_file, + ) + assert test_sp3_file.exists() + + yield { + "test_obs_file": test_obs_file, + } + shutil.rmtree(test_directory) + + +def test_get_index_of_priority_from_filename(): + index_priority = [] + for sp3_filename in [ + "COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz", + "GFZ0MGXRAP_20230010000_01D_05M_ORB.SP3.gz", + "JAX0MGXFIN_20230010000_01D_30S_CLK.CLK.gz", + "WUM0MGXFIN_20230010000_01D_05M_ORB.SP3.gz", + "WUM0MGXFIN_20230010000_01D_30S_CLK.CLK.gz", + ]: + index_priority.append(sp3.get_index_of_priority_from_filename(sp3_filename)) + assert index_priority == [0, 10, 5, 4, 4] + + +def test_download_single_sp3(tmp_path): + """ + This test aims at validating that if a file is reachable with function `sp3.check_online_availability`, + then it can be downloaded from IGS servers. + + Then, in other tests, only the online availability will be tested, to avoid data exchange with IGS servers. + """ + db_folder = tmp_path / "sp3_files" + t_start = pd.Timestamp("2023-01-01 12:00:00") + gps_week, _ = util.timestamp_to_gps_week_and_dow(t_start) + filename_orb, _ = sp3.build_sp3_filename(t_start, sp3.priority[0]) + file_download = sp3.try_downloading_sp3_ftp( + gps_week, sp3.sp3_file_folder(t_start, db_folder), filename_orb + ) + assert file_download.exists() + + file_available = sp3.check_online_availability( + gps_week, sp3.sp3_file_folder(t_start, db_folder), filename_orb + ) + assert file_available is not None + assert file_available == file_download + + +def test_get_sp3_files(set_up_test): + """ + Assert that prx will download the first SP3 file in the priority list if not present. + """ + obs_file = set_up_test["test_obs_file"] + header = georinex.rinexheader(obs_file) + t_start = util.rinex_header_time_string_2_timestamp_ns(header["TIME OF FIRST OBS"]) + t_end = util.rinex_header_time_string_2_timestamp_ns(header["TIME OF LAST OBS"]) + local_db = set_up_test["test_obs_file"].parent + + with ( + patch( # replace download function by online availability check + "prx.precise_corrections.sp3.sp3_file_discovery.try_downloading_sp3_ftp", + new=prx.precise_corrections.sp3.sp3_file_discovery.check_online_availability, + ), + patch( # simulate that local database is empty + "prx.precise_corrections.sp3.sp3_file_discovery.get_local_sp3", + return_value=None, + ), + ): + sp3_files = sp3.get_sp3_files(t_start, t_end, local_db) + + file_orb = sp3_files[0][0].name + file_clk = sp3_files[0][1].name + assert file_orb is not None + assert file_clk is not None + assert sp3.get_index_of_priority_from_filename( + file_orb + ) == sp3.get_index_of_priority_from_filename(file_clk) + + +def test_get_sp3_files_multiple_days(set_up_test): + """ + Assert that prx will download the first SP3 file in the priority list if not present. + """ + obs_file = set_up_test["test_obs_file"] + header = georinex.rinexheader(obs_file) + t_start = util.rinex_header_time_string_2_timestamp_ns(header["TIME OF FIRST OBS"]) + t_end = t_start + pd.Timedelta(days=1) + local_db = set_up_test["test_obs_file"].parent + + with ( + patch( # replace download function by online availability check + "prx.precise_corrections.sp3.sp3_file_discovery.try_downloading_sp3_ftp", + new=prx.precise_corrections.sp3.sp3_file_discovery.check_online_availability, + ), + patch( # simulate that local database is empty + "prx.precise_corrections.sp3.sp3_file_discovery.get_local_sp3", + return_value=None, + ), + ): + sp3_files = sp3.get_sp3_files(t_start, t_end, local_db) + + assert len(sp3_files) == 2 + for ind_day in range(2): + file_orb = sp3_files[ind_day][0].name + file_clk = sp3_files[ind_day][1].name + assert file_orb is not None + assert file_clk is not None + assert sp3.get_index_of_priority_from_filename( + file_orb + ) == sp3.get_index_of_priority_from_filename(file_clk) + + +def test_download_FIN_when_local_RAP_is_available(set_up_test): + """ + The tested scenario is a local database containing the RAP (rapid) SP3 products, but the FIN (final) products are + available online. + The FIN products should be downloaded + """ + obs_file = set_up_test["test_obs_file"] + header = georinex.rinexheader(obs_file) + t_start = util.rinex_header_time_string_2_timestamp_ns(header["TIME OF FIRST OBS"]) + t_end = t_start + local_db = set_up_test["test_obs_file"].parent + + # remove FIN products from local db + for file in local_db.glob("**/*FIN*"): + file.unlink() + + with ( + patch( # replace download function by online availability check + "prx.precise_corrections.sp3.sp3_file_discovery.try_downloading_sp3_ftp", + new=prx.precise_corrections.sp3.sp3_file_discovery.check_online_availability, + ), + ): + sp3_files = sp3.get_sp3_files(t_start, t_end, local_db) + + file_orb = sp3_files[0][0].name + file_clk = sp3_files[0][1].name + assert "FIN" in file_orb + assert "FIN" in file_clk + assert sp3.get_index_of_priority_from_filename( + file_orb + ) == sp3.get_index_of_priority_from_filename(file_clk) + + +def test_match_CLK_and_ORB(set_up_test): + """ + The tested scenario is a local database containing an ORB file, but not the matching CLK file. + The missing CLK file should be downloaded. + """ + obs_file = set_up_test["test_obs_file"] + header = georinex.rinexheader(obs_file) + t_start = util.rinex_header_time_string_2_timestamp_ns(header["TIME OF FIRST OBS"]) + t_end = t_start + local_db = set_up_test["test_obs_file"].parent + + # find an ORB file without matching CLK file + list_orb = list(local_db.glob("**/*ORB*")) + for ind_orb, local_orb in enumerate(list_orb): + idx_priority = sp3.get_index_of_priority_from_filename(local_orb.name) + _, file_clk_expected = sp3.build_sp3_filename( + t_start, sp3.priority[idx_priority] + ) + if not list(local_db.glob("**/file_clk_expected")): + break + + with ( + patch( # replace priority list with single element + "prx.precise_corrections.sp3.sp3_file_discovery.priority", + new=[sp3.priority[idx_priority]], + ), + patch( # replace download function by online availability check + "prx.precise_corrections.sp3.sp3_file_discovery.try_downloading_sp3_ftp", + new=prx.precise_corrections.sp3.sp3_file_discovery.check_online_availability, + ), + ): + sp3_files = sp3.get_sp3_files(t_start, t_end, local_db) + + file_orb = sp3_files[0][0].name + file_clk = sp3_files[0][1].name + assert sp3.get_index_of_priority_from_filename( + file_orb + ) == sp3.get_index_of_priority_from_filename(file_clk) + + +def test_main_sp3_file_discovery_function(set_up_test): + sp3_files = sp3.discover_or_download_sp3_file(set_up_test["test_obs_file"]) + + assert len(sp3_files) != 0 + assert len(sp3_files[0]) == 2 + file_orb = sp3_files[0][0].name + file_clk = sp3_files[0][1].name + assert sp3.get_index_of_priority_from_filename( + file_orb + ) == sp3.get_index_of_priority_from_filename(file_clk) diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index 970402af..0d7395ac 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -8,7 +8,7 @@ parse_rinex_nav_file, ) from prx.rinex_obs.parser import parse_rinex_obs_file -from prx.sp3 import evaluate as sp3_evaluate +from prx.precise_corrections.sp3 import evaluate as sp3_evaluate from prx.rinex_nav import evaluate as rinex_nav_evaluate from prx import constants, converters, util from prx.util import week_and_seconds_2_timedelta diff --git a/src/prx/test/datasets/TLSE00FRA_R_2024001/TLSE00FRA_R_20243510000_01D_30S_MO.crx.gz b/src/prx/test/datasets/TLSE00FRA_R_2024001/TLSE00FRA_R_20243510000_01D_30S_MO.crx.gz new file mode 100644 index 00000000..ab38b955 Binary files /dev/null and b/src/prx/test/datasets/TLSE00FRA_R_2024001/TLSE00FRA_R_20243510000_01D_30S_MO.crx.gz differ diff --git a/src/prx/test/datasets/TLSE_2023001/igr22340.sp3.Z b/src/prx/test/datasets/TLSE_2023001/igr22340.sp3.Z new file mode 100644 index 00000000..6d3a1c5c Binary files /dev/null and b/src/prx/test/datasets/TLSE_2023001/igr22340.sp3.Z differ diff --git a/src/prx/util.py b/src/prx/util.py index a4a4fe75..4ecd66f1 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -15,6 +15,8 @@ from imohash import imohash from astropy.utils import iers from astropy import time as astrotime +from astropy.coordinates import get_sun, ITRS +import astropy.units from prx import constants @@ -419,6 +421,33 @@ def ecef_2_geodetic(pos_ecef): return [latitude_rad, longitude_rad, altitude_m] +def ecef_2_satellite(pos_ecef: np.array, pos_sat_ecef: np.array, epoch: np.array): + """ + Converts an ecef position to the satellite-fixed coordinate frame + + pos_sat_ecef: shape (n,3) + pos_ecef: shape (n,3) + epoch: shape (n,) + """ + k = -pos_sat_ecef / np.linalg.norm(pos_sat_ecef, axis=1).reshape(-1, 1) + pos_sun_ecef = compute_sun_ecef_position(epoch).T + unit_vector_sun_ecef = (pos_sun_ecef - pos_sat_ecef) / np.linalg.norm( + pos_sun_ecef - pos_sat_ecef, axis=1 + ).reshape(-1, 1) + j = np.cross(k, unit_vector_sun_ecef) + j = j / np.linalg.norm(j, axis=1).reshape(-1, 1) + i = np.cross(j, k) + + rot_mat_ecef2sat = np.stack([i, j, k], axis=1) + pos_sat_frame = np.stack( + [ + rot_mat_ecef2sat[i, :, :] @ (pos_ecef[i, :] - pos_sat_ecef[i, :]) + for i in range(pos_sat_ecef.shape[0]) + ] + ) + return pos_sat_frame, rot_mat_ecef2sat + + def obs_dataset_to_obs_dataframe(ds: xarray.Dataset): # Flatten the xarray DataSet into a pandas DataFrame: logger.info("Converting Dataset into flat Dataframe of observations") @@ -482,3 +511,21 @@ def compute_gps_utc_leap_seconds(yyyy: int, doy: int): break assert ~np.isnan(ls), "GPS leap second could not be retrieved" return ls + + +def timestamp_to_gps_week_and_dow(ts: pd.Timestamp) -> tuple[int, int]: + ts_utc = ts.tz_convert("UTC") if ts.tzinfo else ts.tz_localize("UTC") + delta = ts_utc - constants.cGpstUtcEpoch.tz_localize("UTC") + gps_week = delta.days // 7 + dow = delta.days % 7 # day of week + return gps_week, dow + + +def compute_sun_ecef_position(epochs: np.array) -> np.array: + """ + Compute the Sun's ECEF position using the Astropy library + """ + time = astrotime.Time(epochs, scale="utc") + sun_gcrs = get_sun(time) + sun_ecef = sun_gcrs.transform_to(ITRS(obstime=time)) + return sun_ecef.cartesian.xyz.to(astropy.units.meter).value