IGS satellite biases

README.md

@@ -23,10 +23,14 @@ From the `prx` repository root, run
 uv run python src/prx/main.py --observation_file_path <path_to_rinex_file> 
 ```
 
-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 `<path to prx root>/src/prx` to your `PYTHONPATH` environment variable if you run
 into import errors.

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 <path_to_rinex_file> --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`,<br>`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`,<br>`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                                                                                   | ❌      |
 
 
 

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,11 +773,17 @@ 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.")
         sys.exit(1)
     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,
+    )

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