seafloor-geodesy / gnatss

Community Seafloor Global Navigation Satellite Systems - Acoustic (GNSS-A) Transponder Surveying Software
https://gnatss.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
9 stars 12 forks source link

Bug: Mismatch between transponders and replies #296

Open johnbdesanto opened 1 week ago

johnbdesanto commented 1 week ago

I have found the following error after running some data through GNATSS:

Traceback (most recent call last):

  File "/home/jdesanto/miniconda3/envs/gipsyx/bin/gnatss", line 8, in <module>
    sys.exit(app())

  File "/home/jdesanto/miniconda3/envs/gipsyx/lib/python3.10/site-packages/gnatss/cli.py", line 100, in run
    run_gnatss(

  File "/home/jdesanto/miniconda3/envs/gipsyx/lib/python3.10/site-packages/gnatss/main.py", line 110, in run_gnatss
    data_dict = run_solver(config, data_dict, return_raw=return_raw)

  File "/home/jdesanto/miniconda3/envs/gipsyx/lib/python3.10/site-packages/gnatss/solver/run.py", line 33, in run_solver
    process_data, is_converged = prepare_and_solve(all_observations, config, twtt_model=twtt_model)

  File "/home/jdesanto/miniconda3/envs/gipsyx/lib/python3.10/site-packages/gnatss/solver/utilities.py", line 414, in prepare_and_solve
    all_results = perform_solve(

  File "/home/jdesanto/miniconda3/envs/gipsyx/lib/python3.10/site-packages/gnatss/solver/solve.py", line 43, in _calc_tr_vectors
    assert (m, n) == (o, p), f"Mismatch shape found: ({m},{n}) != ({o},{p})"

AssertionError: Mismatch shape found: (3,3) != (2,3)

The issue seems to be that there are pings with fewer replies than transponders in the array. I was under the impression that there was a filter to ignore data that had fewer replies than transponders, but some are leaking through. The gps_solution.csv file for this run was generated using a custom script that was not filtering data for this parameter.

johnbdesanto commented 1 week ago

Update on this bug:

1) I have modified the script I was using to generate gps_solution.csv files to explicitly exclude pings with <3 replies, and this fixed the issue when running GNATSS. However, I think this doesn't really solve the underlying issue since it isn't too hard to imagine a future user trying to process their own data that has pings with empty replies from some transponders. (For a concrete example, many past GNSS-A surveys conducted by Japanese scientists would ping single transponders one at a time in order to get around the issue of distinguish between overlapping replies.)

2) One course of action would be to remove the check in _calc_tr_vectors and allow pings without replies from every transponder be allowed in the inversion. The legacy Chadwell software actually allows for this, but retaining this functionality could have knock-on effects downstream by allowing potentially lower quality pings to be allowed into the inversion as well as making the residuals_enu_components plot more tricky to generate since replies from every transponder are explicitly required for that plot.

3) Another course of action would be to implement a redundant check upstream of _calc_tr_vectors. One sensible place to implement this check imo would be in the get_data_inputs function, defined in ops/data.py. Right now, this function accepts a pandas data frame all_observations as input. I would change this by adding an input parameter to this function that is an integer with the number of transponders. I would add an if statement to the final loop in this function to only write to the output data frame if the number of replies == num_transponders. This would have some knock-on effects with the testing function. The modified code I'm envisioning is something like:

    # Merge all inputs
    for data in zip(
        transmit_xyz, reply_xyz_list, gps_covariance_matrices, observed_delay_list, strict=False
    ):
        if reply_xyz_list.len() == num_transponders:
            data_inputs.append(data)
    return data_inputs

4) Another path forward for maximum flexibility (and complexity) is to allow for a mix of both 2) and 3) above, probably with a new flag in the config file to tell gnatss whether or not to care.

@lsetiawan Thoughts?

lsetiawan commented 1 week ago

Hi John, I really like the option 3. However, rather than setting the integer for the number of transponders, you can get the counts of transponder sets from the send time time and remove ones with less than the max count of transponders. I did have this feature previously, but I think it got lost during the refactor with the new data standard! But I think a separate function in the same data.py file would be good to take in the all_observations dataframe and filter them down.

Here's what it was before:

https://github.com/seafloor-geodesy/gnatss/blob/73cbcac48d725474bd4017004840f8b536a715e0/src/gnatss/main.py#L237-L243