diff --git a/README.md b/README.md index 0ad1840..10c12b3 100644 --- a/README.md +++ b/README.md @@ -23,10 +23,14 @@ From the `prx` repository root, run uv run python src/prx/main.py --observation_file_path ``` -You can specify `--prx_level` according to the type of computation you need. `--prx_level 1` is adapted for DGNSS or RTK -processing, `--prx_level 2` is adapted for SPP processing and is the default value. - -There is an optional argument to select the tropospheric delay model, by adding `--tropo saastamoinen` (default) or `--tropo unb3m` +You can specify `--prx_level` according to the type of computation you need. +`--prx_level 1` is adapted for DGNSS or RTK processing. +`--prx_level 2` is adapted for SPP processing and is the default value. +`--prx_level 3` uses precises corrections from IGS products, adapted for PPP processing. +You can specify the analysis center (among COD, GFZ, GRG, WUM) with the option `--analysis_center GRG`. + +There is an optional argument to select the tropospheric delay model, by adding `--tropo saastamoinen` (default) +or `--tropo unb3m` You might have to add `/src/prx` to your `PYTHONPATH` environment variable if you run into import errors. diff --git a/documents/dev_status.md b/documents/dev_status.md index 0e75f74..786a6a2 100644 --- a/documents/dev_status.md +++ b/documents/dev_status.md @@ -52,18 +52,18 @@ Same as level 1 and 2, with the following additional parameters. The parameters uv run python src/prx/main.py --observation_file_path --prx_level 3 ``` -| Parameters | Name in PRX file | Status | -|------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------|--------| -| Satellite antenna phase center position and velocity | `sat_pos_x_m`, `sat_pos_y_m`, `sat_pos_z_m`,
`sat_vel_x_mps`, `sat_vel_y_mps`, `sat_vel_z_mps` | ✅ | -| Satellite center of mass position (computed using `sp3` orbit files) | `sat_pos_com_x_m`, `sat_pos_com_y_m`, `sat_pos_com_z_m` | ✅ | -| Satellite clock offset and drift (including relativistic effect) (computed using `sp3` clock files) | `sat_clock_offset_m`, `sat_clock_drift_mps` | ✅ | -| Tropospheric delay (computed using `tropex` files) | `tropo_delay_m` | ❌ | -| Tropospheric mapping function | to be completed | ❌ | -| Ionospheric delay (computed using `ionex` files) | `iono_delay_m` | ❌ | -| Satellite code & phase bias (computed using `bia` files) | `sat_code_bias_m` | ❌ | -| Satellite & receiver phase center offset (computed using `antex` files) | `sat_pco_x_m`, `sat_pco_y_m`, `sat_pco_z_m` | ✅ | -| Satellite & receiver phase center variation (computed using `antex` files) | to be completed | ❌ | -| Solid Earth Tide | to be completed | ❌ | +| Parameters | Name in PRX file | Status | +|-----------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------|--------| +| Satellite antenna phase center position and velocity | `sat_pos_x_m`, `sat_pos_y_m`, `sat_pos_z_m`,
`sat_vel_x_mps`, `sat_vel_y_mps`, `sat_vel_z_mps` | ✅ | +| Satellite center of mass position (computed using `sp3` orbit files) | `sat_pos_com_x_m`, `sat_pos_com_y_m`, `sat_pos_com_z_m` | ✅ | +| Satellite clock offset and drift (including relativistic effect) (computed using `sp3` clock files) | `sat_clock_offset_m`, `sat_clock_drift_mps` | ✅ | +| Tropospheric delay (computed using `tropex` files) | `tropo_delay_m` | ❌ | +| Tropospheric mapping function | to be completed | ❌ | +| Ionospheric delay (computed using `ionex` files) | `iono_delay_m` | ❌ | +| Satellite code & phase bias (computed using `bia` files) | `sat_code_bias_m` , `sat_carrier_bias_m` | ✅ | +| Satellite & receiver phase center offset (computed using `antex` files) | `sat_pco_x_m`, `sat_pco_y_m`, `sat_pco_z_m` | ✅ | +| Satellite & receiver phase center variation (computed using `antex` files) | to be completed | ❌ | +| Solid Earth Tide | to be completed | ❌ | diff --git a/src/prx/main.py b/src/prx/main.py index 28e3a3d..99fe0f9 100644 --- a/src/prx/main.py +++ b/src/prx/main.py @@ -11,6 +11,7 @@ import prx.util as util from prx import atmospheric_corrections as atmo from prx.constants import carrier_frequencies_hz, cDegPerRad +from prx.precise_corrections.bia import bia_file_discovery, bia_processing from prx.rinex_obs.parser import parse_rinex_obs_file, get_glonass_slot from prx.util import is_rinex_3_obs_file, is_rinex_3_nav_file, configure_logging from prx.rinex_nav import nav_file_discovery @@ -21,6 +22,7 @@ from prx.precise_corrections.sp3 import sp3_file_discovery from prx.precise_corrections.antex import antex_file_discovery + log = logging.getLogger(__name__) @@ -435,11 +437,12 @@ def build_records_levels_12( def build_records_level_3( - rinex_3_obs_file, - sp3_orbit_files, - atx_file, - approximate_receiver_ecef_position_m, - model_tropo, + rinex_3_obs_file: Path, + sp3_orbit_files: list[Path], + atx_file: Path, + bia_files: list[Path], + approximate_receiver_ecef_position_m: list, + model_tropo: str, ): """ Creates a flat_obs dataframe including columns for prx processing level 3. @@ -528,12 +531,13 @@ def build_records_level_3( }, ) - # Compute broadcast position, velocity, clock offset, clock offset rate and TGDs + # Compute broadcast position, velocity, clock offset, clock offset rate and hardware biases sat_states_per_day = [] - for file in sp3_orbit_files: + for file_sp3, file_bia in zip(sp3_orbit_files, bia_files): # get year and doy from sp3 orb filename - year = int(file.name[11:15]) - doy = int(file.name[15:18]) + year = int(file_sp3.name[11:15]) + doy = int(file_sp3.name[15:18]) + # create query for single day day_query = query.loc[ ( query.query_time_isagpst @@ -547,18 +551,20 @@ def build_records_level_3( if day_query.empty: continue + # compute satellite position, velocity and clock log.info(f"Computing satellite states for {year}-{doy:03d}") - sat_states_per_day.append( - sp3_evaluate.compute( - file, - day_query, - atx_file, - ).assign( # TODO: add hw satellite biases - sat_code_bias_m=np.nan, - sat_carrier_bias_m=np.nan, - ) + sat_states_single_day = sp3_evaluate.compute(file_sp3, day_query, atx_file) + # add satellite hardware biases + sat_bias = bia_processing.compute_sat_hw_biases( + sat_states_single_day, bia_processing.parse_bia_file(file_bia) + ).to_pandas() + sat_states_single_day = sat_states_single_day.merge( + sat_bias, on=["sv", "signal", "query_time_isagpst"], how="left" ) + # collect satellite states per day + sat_states_per_day.append(sat_states_single_day) + sat_states = pd.concat(sat_states_per_day) sat_states = sat_states.rename( columns={ @@ -567,12 +573,11 @@ def build_records_level_3( "query_time_isagpst": "time_of_emission_isagpst", }, ) - ( - sat_states["elevation_rad"], - sat_states["azimuth_rad"], - ) = util.compute_satellite_elevation_and_azimuth( - sat_states[["sat_pos_x_m", "sat_pos_y_m", "sat_pos_z_m"]].to_numpy(), - approximate_receiver_ecef_position_m, + (sat_states["elevation_rad"], sat_states["azimuth_rad"]) = ( + util.compute_satellite_elevation_and_azimuth( + sat_states[["sat_pos_x_m", "sat_pos_y_m", "sat_pos_z_m"]].to_numpy(), + approximate_receiver_ecef_position_m, + ) ) # Compute anything else that is satellite-specific @@ -652,6 +657,7 @@ def process( observation_file_path: Path, prx_level=2, model_tropo="saastamoinen", + analysis_center="COD", joblib_backend: str = "loky", ): t0 = pd.Timestamp.now() @@ -688,11 +694,17 @@ def process( aux_files = {} # define auxiliary files aux_files["sp3_orb"], aux_files["sp3_clk"] = ( - sp3_file_discovery.discover_or_download_sp3_file(rinex_3_obs_file) + sp3_file_discovery.discover_or_download_sp3_file( + rinex_3_obs_file, analysis_center + ) ) aux_files["atx"] = antex_file_discovery.discover_or_download_atx_file( rinex_3_obs_file ) + aux_files["bia"] = [ + bia_file_discovery.discover_or_download_bia_file(sp3_file) + for sp3_file in aux_files["sp3_orb"] + ] # define metadata metadata = build_metadata({"obs_file": rinex_3_obs_file, "nav_file": []}) @@ -704,6 +716,7 @@ def process( rinex_3_obs_file, aux_files["sp3_orb"], aux_files["atx"], + aux_files["bia"], metadata["approximate_receiver_ecef_position_m"], model_tropo, ) @@ -738,6 +751,13 @@ def process( choices=[1, 2, 3], default=2, ) + parser.add_argument( + "--analysis_center", + type=str, + help="Analysis center as source for precise correction (cod, gfz, grg, wum)", + choices=["cod", "gfz", "grg", "wum"], + default="cod", + ) parser.add_argument( "--tropo", type=str, @@ -753,6 +773,7 @@ def process( required=False, ) args = parser.parse_args() + configure_logging(args.log_level) if args.observation_file_path is None: log.error("No observation file path provided.") @@ -760,4 +781,9 @@ def process( if not Path(args.observation_file_path).exists(): log.error(f"Observation file {args.observation_file_path} does not exist.") sys.exit(1) - process(Path(args.observation_file_path), args.prx_level, args.tropo) + process( + Path(args.observation_file_path), + args.prx_level, + args.tropo, + args.analysis_center, + ) diff --git a/src/prx/precise_corrections/bia/bia_file_discovery.py b/src/prx/precise_corrections/bia/bia_file_discovery.py new file mode 100644 index 0000000..43f6b31 --- /dev/null +++ b/src/prx/precise_corrections/bia/bia_file_discovery.py @@ -0,0 +1,107 @@ +import ftplib +import logging +import urllib +from pathlib import Path + +import pandas as pd + +from prx import util + +log = logging.getLogger(__name__) + + +def bia_file_database_folder(): + """ + Returns the path to the folder where ATX database files are stored. + """ + db_folder = util.prx_src_directory() / "precise_corrections/bia/bia_files" + db_folder.mkdir(exist_ok=True, parents=True) + return db_folder + + +def build_bia_file_name(year: int, doy: int, analysis_center: str): + # look-up table providing type of product for each analysic center + ac2type = { + "COD": "FIN", + "GFZ": "RAP", + "GRG": "FIN", + "WUM": "RAP", + } + analysis_center = analysis_center.upper() + return f"{analysis_center}0MGX{ac2type[analysis_center]}_{year}{doy:03d}0000_01D_01D_OSB.BIA.gz" + + +def bia_file_folder(year: int, doy: int): + folder = bia_file_database_folder() / f"{year}/{doy:03d}" + folder.mkdir(parents=True, exist_ok=True) + return folder + + +def get_local_bia_file(year: int, doy: int, analysis_center: str) -> Path | None: + local_file = bia_file_folder(year, doy) / build_bia_file_name( + year, doy, analysis_center + ) + if local_file.exists(): + return local_file + else: + return None + + +def check_online_availability(year: int, doy: int, analysis_center: str) -> Path | None: + """ + Need to keep the same inputs as try_downloading_bia_ftp, in order to be able to use `unittest.mock.patch` in tests + """ + server = "gssc.esa.int" + gps_week, _ = util.timestamp_to_gps_week_and_dow( + pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1) + ) + remote_folder = f"gnss/products/{gps_week}" + file = build_bia_file_name(year, doy, analysis_center) + ftp = ftplib.FTP(server) + ftp.login() + ftp.cwd(remote_folder) + try: + ftp.size(file) + return bia_file_folder(year, doy) / build_bia_file_name( + year, doy, analysis_center + ) + except ftplib.error_perm: + log.warning(f"{file} not available on {server}") + return None + + +def try_downloading_bia_ftp(year: int, doy: int, analysis_center: str) -> Path | None: + server = "gssc.esa.int" + file = build_bia_file_name(year, doy, analysis_center) + gps_week, _ = util.timestamp_to_gps_week_and_dow( + pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1) + ) + remote_folder = f"gnss/products/{gps_week}" + ftp_file = f"ftp://{server}/{remote_folder}/{file}" + local_file = bia_file_folder(year, doy) / build_bia_file_name( + year, doy, analysis_center + ) + 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 bia file: {ftp_file}") + return local_file + + +def discover_or_download_bia_file(sp3_file_path: Path) -> Path | None: + log.info(f"Finding bia files for {sp3_file_path} ...") + year = int(sp3_file_path.name[11:15]) + doy = int(sp3_file_path.stem[15:18]) + analysis_center = sp3_file_path.name[0:3] + + local_file = get_local_bia_file(year, doy, analysis_center) + if local_file: + log.info(f"Found local bia file: {local_file}") + return local_file + + downloaded_file = try_downloading_bia_ftp(year, doy, analysis_center) + if downloaded_file: + return downloaded_file + + return None diff --git a/src/prx/precise_corrections/bia/bia_processing.py b/src/prx/precise_corrections/bia/bia_processing.py new file mode 100644 index 0000000..646eabe --- /dev/null +++ b/src/prx/precise_corrections/bia/bia_processing.py @@ -0,0 +1,130 @@ +from pathlib import Path +import polars as pl +import pandas as pd +from prx import constants, converters, util +from datetime import datetime, timedelta + + +def parse_bia_file(filepath_bia_gz: Path) -> pl.DataFrame: + """ + Parse gzipped BIA file and returns a pl.DataFrame with columns: + - sat_id: + - obs_id: rinex obs identifier + - sat_hw_bias_m: bias value in meters + - start: timestamp + - end: timestamp + + Example: + ┌────────┬────────┬───────────────┬─────────────────────┬─────────────────────┐ + │ sat_id ┆ obs_id ┆ sat_hw_bias_m ┆ start ┆ end │ + │ --- ┆ --- ┆ --- ┆ --- ┆ --- │ + │ str ┆ str ┆ f64 ┆ datetime[ns] ┆ datetime[ns] │ + ╞════════╪════════╪═══════════════╪═════════════════════╪═════════════════════╡ + │ G01 ┆ C1C ┆ -0.425645 ┆ 2023-01-01 00:00:00 ┆ 2023-01-02 00:00:00 │ + │ G01 ┆ C1W ┆ -0.0 ┆ 2023-01-01 00:00:00 ┆ 2023-01-02 00:00:00 │ + │ … ┆ … ┆ … ┆ … ┆ … │ + │ E36 ┆ L5Q ┆ -0.018269 ┆ 2023-01-01 00:00:00 ┆ 2023-01-02 00:00:00 │ + │ E36 ┆ L5X ┆ -0.018269 ┆ 2023-01-01 00:00:00 ┆ 2023-01-02 00:00:00 │ + └────────┴────────┴───────────────┴─────────────────────┴─────────────────────┘ + + """ + + @util.disk_cache.cache(ignore=["filepath_bia_gz"]) + def cached_load(filepath_bia_gz: Path, file_hash: str): + filepath_bia = converters.compressed_to_uncompressed(filepath_bia_gz) + with open(filepath_bia, "r", encoding="cp1250") as f: + sat_id_list = [] + obs1_list = [] + val_list = [] + start_list = [] + end_list = [] + # find beginning of block BIAS/SOLUTION + for line in f: + if line.startswith("+BIAS/SOLUTION"): + break + + for line in f: + if line.startswith("-BIAS/SOLUTION"): # escape loop + break + if line.startswith("*BIAS"): # skip header + continue + bias_type = line[0:4].strip() + station = line[15:24].strip() + if (station == "") and ( + bias_type == "OSB" + ): # keep only OSB and satellite biases + unit = line[64:69].strip() + assert unit == "ns", ( + f"Wrong unit in file. Expected 'ns', read '{unit}'" + ) + sat_id = line[11:14].strip() + obs1 = line[25:29].strip() + estimated_value = ( + float(line[70:91]) + / constants.cNanoSecondsPerSecond + * constants.cGpsSpeedOfLight_mps + ) + start = datetime(int(line[35:39]), 1, 1) + timedelta( + days=int(line[40:43]) - 1, seconds=int(line[44:49]) + ) + end = datetime(int(line[50:54]), 1, 1) + timedelta( + days=int(line[55:58]) - 1, seconds=int(line[59:64]) + ) + sat_id_list.append(sat_id) + obs1_list.append(obs1) + val_list.append(estimated_value) + start_list.append(start) + end_list.append(end) + bia_df = pl.DataFrame( + { + "sat_id": sat_id_list, + "obs_id": obs1_list, + "sat_hw_bias_m": val_list, + "start": start_list, + "end": end_list, + }, + ).with_columns( + pl.col("start").cast(pl.Datetime("ns")), + pl.col("end").cast(pl.Datetime("ns")), + ) + return bia_df + + file_content_hash = util.hash_of_file_content(filepath_bia_gz) + return cached_load(filepath_bia_gz, file_content_hash) + + +def compute_sat_hw_biases(query: pd.DataFrame, bia_df: pl.DataFrame) -> pd.DataFrame: + query_sorted = pl.from_pandas(query).sort("query_time_isagpst") + bia_sorted = bia_df.sort("start") + + sat_code_bias = query_sorted.join_asof( + bia_sorted, + left_on="query_time_isagpst", + right_on="start", + by_left=["sv", "signal"], + by_right=["sat_id", "obs_id"], + strategy="backward", + ).get_column("sat_hw_bias_m") + + sat_carrier_bias = ( + query_sorted.with_columns( + (pl.lit("L") + pl.col("signal").str.slice(1)).alias("signal") + ) + .join_asof( + bia_sorted, + left_on="query_time_isagpst", + right_on="start", + by_left=["sv", "signal"], + by_right=["sat_id", "obs_id"], + strategy="backward", + ) + .get_column("sat_hw_bias_m") + ) + + sat_bias = query_sorted.with_columns( + sat_code_bias.alias("sat_code_bias_m"), + sat_carrier_bias.alias("sat_carrier_bias_m"), + ) + return sat_bias[ + ["sv", "signal", "query_time_isagpst", "sat_code_bias_m", "sat_carrier_bias_m"] + ] diff --git a/src/prx/precise_corrections/bia/test/test_bia_discovery.py b/src/prx/precise_corrections/bia/test/test_bia_discovery.py new file mode 100644 index 0000000..d272d35 --- /dev/null +++ b/src/prx/precise_corrections/bia/test/test_bia_discovery.py @@ -0,0 +1,63 @@ +import shutil +from pathlib import Path + +import pandas as pd +import pytest + +from prx import util +from prx.precise_corrections.sp3.sp3_file_discovery import sp3_file_folder + +from prx.precise_corrections.bia import bia_file_discovery as discovery + + +@pytest.fixture +def set_up_test(tmp_path_factory): + test_sp3_file = ( + sp3_file_folder(pd.Timestamp("2023-01-01")) + / "COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz" + ) + shutil.copy( + util.prx_src_directory().joinpath( + "test", + "datasets", + "TLSE_2023001", + "COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz", + ), + test_sp3_file, + ) + assert test_sp3_file.exists() + + yield { + "year": 2023, + "doy": 1, + "analysis_center": "COD", + "sp3": test_sp3_file, + } + + +def test_bia_file_online_availability(set_up_test): + year = set_up_test["year"] + doy = set_up_test["doy"] + ac = set_up_test["analysis_center"] + assert discovery.check_online_availability( + year, doy, ac + ) == discovery.bia_file_folder(year, doy) / discovery.build_bia_file_name( + year, doy, ac + ) + + +def test_download_bia_file(): + # delete local file to trigger download + path_local = discovery.bia_file_folder(2023, 1) / discovery.build_bia_file_name( + 2023, 1, "COD" + ) + path_local.unlink(missing_ok=True) + file = discovery.discover_or_download_bia_file( + Path("COD0MGXFIN_20230010000_01D_05M_ORB.SP3.gz") + ) + assert file.exists() + + +def test_find_local_bia_file(set_up_test): + file = discovery.discover_or_download_bia_file(set_up_test["sp3"]) + assert file.exists() diff --git a/src/prx/precise_corrections/bia/test/test_bia_processing.py b/src/prx/precise_corrections/bia/test/test_bia_processing.py new file mode 100644 index 0000000..9ade6ba --- /dev/null +++ b/src/prx/precise_corrections/bia/test/test_bia_processing.py @@ -0,0 +1,146 @@ +import shutil + +import numpy as np +import pandas as pd +import polars as pl + +import pytest + +from prx import util, constants +import prx.precise_corrections.bia.bia_processing as bia + + +@pytest.fixture(scope="session") +def input_for_test(tmp_path_factory): + temp_directory = tmp_path_factory.mktemp("test_inputs") + src_dir = util.prx_src_directory() + test_dir = src_dir.joinpath("test", "datasets") + test_files = { + "cod": temp_directory / "COD0MGXFIN_20230010000_01D_01D_OSB.BIA.gz", + "gfz": temp_directory / "GFZ0MGXRAP_20230010000_01D_01D_OSB.BIA.gz", + "grg": temp_directory / "GRG0MGXFIN_20240010000_01D_01D_OSB.BIA.gz", + "wum": temp_directory / "WUM0MGXRAP_20240010000_01D_01D_OSB.BIA.gz", + } + shutil.copy( + test_dir.joinpath("TLSE_2023001", test_files["cod"].name), test_files["cod"] + ) + shutil.copy( + test_dir.joinpath("TLSE_2023001", test_files["gfz"].name), test_files["gfz"] + ) + shutil.copy( + test_dir.joinpath("TLSE00FRA_R_2024001", test_files["grg"].name), + test_files["grg"], + ) + shutil.copy( + test_dir.joinpath("TLSE00FRA_R_2024001", test_files["wum"].name), + test_files["wum"], + ) + for test_file_path in test_files.values(): + assert test_file_path.exists() + yield test_files + shutil.rmtree(temp_directory) + + +def test_bia_parsing(input_for_test): + bia_df = bia.parse_bia_file(input_for_test["cod"]) + # manual check in file + assert ( + bia_df.filter((pl.col("sat_id") == "G01") & (pl.col("obs_id") == "C1C")).item( + 0, "sat_hw_bias_m" + ) + == -1.4198 / constants.cNanoSecondsPerSecond * constants.cGpsSpeedOfLight_mps + ) + assert ( + bia_df.filter((pl.col("sat_id") == "G01") & (pl.col("obs_id") == "C1W")).item( + 0, "sat_hw_bias_m" + ) + == -0.0000 / constants.cNanoSecondsPerSecond * constants.cGpsSpeedOfLight_mps + ) + + +def test_iono_free_code_bias(input_for_test): + """ + The iono-free combination of the satellite hardware code biases should be close to 0 + (or close to a constant across all satellites from the same constellation) + """ + threshold = 1e-4 + obs_used_for_if = { + "cod": { + "G": ["C1W", "C2W"], + "E": ["C1C", "C5Q"], + "C": ["C2I", "C6I"], + "J": ["C1C", "C2L"], + }, + "gfz": { + "G": ["C1W", "C2W"], + "E": ["C1C", "C5Q"], + "C": ["C2I", "C6I"], + # "J": ["C1C", "C2L"], + }, + "grg": { + "G": ["L1W", "L2W"], + "E": ["L1X", "L5X"], + "C": ["L2I", "L6I"], + }, + "wum": { + "G": ["C1W", "C2W"], + "R": ["C1P", "C2P"], + "E": ["C1C", "C5Q"], + "C": ["C2I", "C6I"], + # "J": ["C1X", "C2X"], + }, + } + for ac in obs_used_for_if: + print(f"=== Testing analysis center {ac} ===") + bia_df = ( + bia.parse_bia_file(input_for_test[ac]) + .filter(pl.col("start") == pl.col("start").unique().min()) + .pivot(on="obs_id", index="sat_id", values="sat_hw_bias_m") + ) + for const in obs_used_for_if[ac]: + # Iono-free combination of COD for GPS uses C1W and C2W + freq_id1 = obs_used_for_if[ac][const][0][1] + freq_id2 = obs_used_for_if[ac][const][1][1] + f1 = constants.carrier_frequencies_hz()[const]["L" + freq_id1][1] + f2 = constants.carrier_frequencies_hz()[const]["L" + freq_id2][1] + bia_if_gps = bia_df.filter(pl.col("sat_id").str.starts_with(const)).select( + pl.col("sat_id"), + ( + ( + f1**2 * pl.col(obs_used_for_if[ac][const][0]) + - f2**2 * pl.col(obs_used_for_if[ac][const][1]) + ) + / (f1**2 - f2**2) + ).alias("if_code_bias"), + ) + print( + f"Maximum iono-free code bias for {const}: {bia_if_gps['if_code_bias'].max()} m" + ) + assert (bia_if_gps["if_code_bias"] < threshold).all() + + +def test_retrieve_satellite_biases(input_for_test): + # choose a query where 2 different biases exists for the same sat/sig at different times + query = pd.DataFrame( + { + "sv": ["G03", "G03"], + "signal": ["C5Q", "C5Q"], + "query_time_isagpst": [ + pd.Timestamp("2024-01-01 00:00:00"), + pd.Timestamp("2024-01-01 00:15:00"), + ], + } + ) + + sat_bias = bia.compute_sat_hw_biases( + query, bia.parse_bia_file(input_for_test["wum"]) + ) + + assert "sat_code_bias_m" in sat_bias.columns + assert "sat_carrier_bias_m" in sat_bias.columns + # manual check in bia file + assert sat_bias[ + ["sat_code_bias_m", "sat_carrier_bias_m"] + ].to_numpy() == pytest.approx( + np.array([[1.84671646, 0.38192426], [1.84671646, 0.39412953]]) + ) diff --git a/src/prx/precise_corrections/sp3/sp3_file_discovery.py b/src/prx/precise_corrections/sp3/sp3_file_discovery.py index 7999850..8c479ce 100644 --- a/src/prx/precise_corrections/sp3/sp3_file_discovery.py +++ b/src/prx/precise_corrections/sp3/sp3_file_discovery.py @@ -130,8 +130,9 @@ def try_downloading_sp3_ftp(gps_week: int, folder: Path, file: str) -> Path | No 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(): + try: + urllib.request.urlretrieve(ftp_file, local_compressed_file) + except urllib.request.URLError: log.warning(f"Could not download {ftp_file}") return None local_file = converters.compressed_to_uncompressed(local_compressed_file) @@ -143,6 +144,7 @@ def try_downloading_sp3_ftp(gps_week: int, folder: Path, file: str) -> Path | No def get_sp3_files( mid_day_start: pd.Timestamp, mid_day_end: pd.Timestamp, + analysis_center: str, db_folder=sp3_file_database_folder(), ) -> tuple[list[Path | None], list[Path | None]]: sp3_orb_files = [] @@ -151,31 +153,33 @@ def get_sp3_files( 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_orb_files.append(file_orb) - sp3_clk_files.append(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_orb_files.append(None) - sp3_clk_files.append(None) + if analysis_center == p[0]: + 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_orb_files.append(file_orb) + sp3_clk_files.append(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_orb_files.append(None) + sp3_clk_files.append(None) date += pd.Timedelta(1, unit="days") return sp3_orb_files, sp3_clk_files def discover_or_download_sp3_file( observation_file_path=Path, + analysis_center="COD", ) -> tuple[list[Path | None], list[Path | None]]: """ Returns the path to a valid SP3 file (local or downloaded) corresponding to the observation file. @@ -193,5 +197,5 @@ def discover_or_download_sp3_file( util.rinex_header_time_string_2_timestamp_ns(header["TIME OF LAST OBS"]) ) - sp3_orb_files, sp3_clk_files = get_sp3_files(t_start, t_end) + sp3_orb_files, sp3_clk_files = get_sp3_files(t_start, t_end, analysis_center) return sp3_orb_files, sp3_clk_files 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 index 935e9fe..4573edc 100644 --- a/src/prx/precise_corrections/sp3/test/test_sp3_file_discovey.py +++ b/src/prx/precise_corrections/sp3/test/test_sp3_file_discovey.py @@ -110,7 +110,9 @@ def test_get_sp3_files(set_up_test): return_value=None, ), ): - sp3_orb_files, sp3_clk_files = sp3.get_sp3_files(t_start, t_end, local_db) + sp3_orb_files, sp3_clk_files = sp3.get_sp3_files( + t_start, t_end, "COD", local_db + ) file_orb = sp3_orb_files[0].name file_clk = sp3_clk_files[0].name @@ -141,7 +143,7 @@ def test_get_sp3_files_multiple_days(set_up_test): return_value=None, ), ): - sp3_files = sp3.get_sp3_files(t_start, t_end, local_db) + sp3_files = sp3.get_sp3_files(t_start, t_end, "COD", local_db) assert len(sp3_files) == 2 for ind_day in range(2): @@ -176,7 +178,9 @@ def test_download_FIN_when_local_RAP_is_available(set_up_test): new=prx.precise_corrections.sp3.sp3_file_discovery.check_online_availability, ), ): - sp3_orb_files, sp3_clk_files = sp3.get_sp3_files(t_start, t_end, local_db) + sp3_orb_files, sp3_clk_files = sp3.get_sp3_files( + t_start, t_end, "COD", local_db + ) file_orb = sp3_orb_files[0].name file_clk = sp3_clk_files[0].name @@ -205,6 +209,7 @@ def test_match_CLK_and_ORB(set_up_test): _, file_clk_expected = sp3.build_sp3_filename( t_start, sp3.priority[idx_priority] ) + ac_code = sp3.priority[idx_priority][0] if not list(local_db.glob("**/file_clk_expected")): break @@ -218,7 +223,9 @@ def test_match_CLK_and_ORB(set_up_test): new=prx.precise_corrections.sp3.sp3_file_discovery.check_online_availability, ), ): - sp3_orb_files, sp3_clk_files = sp3.get_sp3_files(t_start, t_end, local_db) + sp3_orb_files, sp3_clk_files = sp3.get_sp3_files( + t_start, t_end, ac_code, local_db + ) file_orb = sp3_orb_files[0].name file_clk = sp3_clk_files[0].name diff --git a/src/prx/rinex_nav/test/datasets/COD0MGXFIN_20220010000_01D_01D_OSB.BIA.gz b/src/prx/rinex_nav/test/datasets/COD0MGXFIN_20220010000_01D_01D_OSB.BIA.gz new file mode 100644 index 0000000..9fa4cd3 Binary files /dev/null and b/src/prx/rinex_nav/test/datasets/COD0MGXFIN_20220010000_01D_01D_OSB.BIA.gz differ diff --git a/src/prx/test/datasets/TLSE00FRA_R_2024001/GRG0MGXFIN_20240010000_01D_01D_OSB.BIA.gz b/src/prx/test/datasets/TLSE00FRA_R_2024001/GRG0MGXFIN_20240010000_01D_01D_OSB.BIA.gz new file mode 100644 index 0000000..abf9be4 Binary files /dev/null and b/src/prx/test/datasets/TLSE00FRA_R_2024001/GRG0MGXFIN_20240010000_01D_01D_OSB.BIA.gz differ diff --git a/src/prx/test/datasets/TLSE00FRA_R_2024001/WUM0MGXRAP_20240010000_01D_01D_OSB.BIA.gz b/src/prx/test/datasets/TLSE00FRA_R_2024001/WUM0MGXRAP_20240010000_01D_01D_OSB.BIA.gz new file mode 100644 index 0000000..ff83ce6 Binary files /dev/null and b/src/prx/test/datasets/TLSE00FRA_R_2024001/WUM0MGXRAP_20240010000_01D_01D_OSB.BIA.gz differ diff --git a/src/prx/test/datasets/TLSE_2023001/COD0MGXFIN_20230010000_01D_01D_OSB.BIA.gz b/src/prx/test/datasets/TLSE_2023001/COD0MGXFIN_20230010000_01D_01D_OSB.BIA.gz new file mode 100644 index 0000000..29a23f9 Binary files /dev/null and b/src/prx/test/datasets/TLSE_2023001/COD0MGXFIN_20230010000_01D_01D_OSB.BIA.gz differ diff --git a/src/prx/test/datasets/TLSE_2023001/GFZ0MGXRAP_20230010000_01D_01D_OSB.BIA.gz b/src/prx/test/datasets/TLSE_2023001/GFZ0MGXRAP_20230010000_01D_01D_OSB.BIA.gz new file mode 100644 index 0000000..756d795 Binary files /dev/null and b/src/prx/test/datasets/TLSE_2023001/GFZ0MGXRAP_20230010000_01D_01D_OSB.BIA.gz differ diff --git a/src/prx/test/test_main.py b/src/prx/test/test_main.py index 55e5d30..30d58b5 100644 --- a/src/prx/test/test_main.py +++ b/src/prx/test/test_main.py @@ -12,6 +12,7 @@ from prx.main import write_prx_file from prx.precise_corrections.antex.antex_file_discovery import atx_file_database_folder from prx.precise_corrections.sp3.sp3_file_discovery import sp3_file_database_folder +from prx.precise_corrections.bia.bia_file_discovery import bia_file_database_folder from prx.rinex_nav import nav_file_discovery log = logging.getLogger(__name__) @@ -38,10 +39,7 @@ def input_for_test_tlse(tmp_path_factory): ) rnx_nav = datasets_directory / "TLSE_2023001/BRDC00IGS_R_20230010000_01D_MN.rnx.gz" for file in [compressed_crx, rnx_nav]: - shutil.copy( - file, - test_directory / file.name, - ) + shutil.copy(file, test_directory / file.name) # copy and uncompress precise correction files to local database os.makedirs(sp3_file_database_folder() / "2023/001/", exist_ok=True) @@ -66,6 +64,11 @@ def input_for_test_tlse(tmp_path_factory): atx_local = shutil.copy(atx, atx_file_database_folder() / atx.name) assert atx_local.exists() + os.makedirs(bia_file_database_folder() / "2023/001/", exist_ok=True) + bia = datasets_directory / "TLSE_2023001/COD0MGXFIN_20230010000_01D_01D_OSB.BIA.gz" + bia_local = shutil.copy(bia, bia_file_database_folder() / "2023/001" / bia.name) + assert bia_local.exists() + yield test_directory / compressed_crx.name shutil.rmtree(test_directory) @@ -89,10 +92,7 @@ def input_for_test_tlse_2024(tmp_path_factory): datasets_directory / "TLSE00FRA_R_2024001/BRDC00IGS_R_20240010000_01D_MN.rnx.gz" ) for file in [compressed_compact_rinex_file, ephemerides_file]: - shutil.copy( - file, - test_directory / file.name, - ) + shutil.copy(file, test_directory / file.name) yield test_directory / compressed_compact_rinex_file.name shutil.rmtree(test_directory) @@ -269,11 +269,7 @@ def test_spp_lsq_nist(input_for_test_nist): == df.time_of_reception_in_receiver_time.min() ] for constellations_to_use in [ - ( - "G", - "E", - "C", - ), + ("G", "E", "C"), ("G", "S"), ("G",), ("E",), @@ -305,11 +301,7 @@ def test_spp_lsq_tlse(input_for_test_tlse): & (df.sat_elevation_deg > 10) ] for constellations_to_use in [ - ( - "G", - "E", - "C", - ), + ("G", "E", "C"), ("G", "S"), ("G",), ("E",), @@ -375,7 +367,7 @@ def test_spp_lsq_tlse_single_freq(input_for_test_tlse): assert np.max(np.abs(velocity_offset)) < 3e-2 -def test_spp_lsq_tlse_with_precise_corrections(input_for_test_tlse): +def test_spp_lsq_tlse_iono_free_with_precise_corrections(input_for_test_tlse): """ Use iono-free combinations considered by IGS conventions (CODE Analysis Center): | Constellation | Frequency pair |