aodn/imos-toolbox

Check of ADCP bin mapping algorithm

Closed this issue · 7 comments

BecCowley commented

I've been reading a lot of literature to check the bin mapping algorithm used in adcpBinMappingPP.m. These are the lines I have reviewed:

    number_of_beams = 4;
    %TODO: block function.
    CP = cos(pitch);
    % H[TxB] = P[T] x (+-R[T] x B[B])
    nonMappedHeightAboveSensorBeam1 = CP .* (cos(beamAngle + roll) / cos(beamAngle) .* distAlongBeams);
    nonMappedHeightAboveSensorBeam2 = CP .* (cos(beamAngle - roll) / cos(beamAngle) .* distAlongBeams);

    % H[TxB] = R[T] x (-+P[T] x B[B])
    CR = cos(roll);
    nonMappedHeightAboveSensorBeam3 = CR .* (cos(beamAngle - pitch) / cos(beamAngle) .* distAlongBeams);
    nonMappedHeightAboveSensorBeam4 = CR .* (cos(beamAngle + pitch) / cos(beamAngle) .* distAlongBeams);

I have two concerns.

  1. The algorithm only works for up-facing ADCPs.
  2. The multiplication with CP and CR is incorrect.

Here is what I think should be happening:
First, account for up or down-facing sign for pitch:

if ~isUpwardLooking %for down facing adcps
    sg1 = 1; sg3 = 1;
else %for up looking
    sg1 = 1; sg3 = -1;
end

Then, the algorithm should look like this:

    number_of_beams = 4;
    %TODO: block function.
    % H[TxB] = P[T] x (+-R[T] x B[B])
    nonMappedHeightAboveSensorBeam1 = (cos(beamAngle + sg1*roll) / cos(beamAngle) .* distAlongBeams);
    nonMappedHeightAboveSensorBeam2 = (cos(beamAngle + -sg1*roll) / cos(beamAngle) .* distAlongBeams);

    % H[TxB] = R[T] x (-+P[T] x B[B])
    nonMappedHeightAboveSensorBeam3 = (cos(beamAngle + sg3*pitch) / cos(beamAngle) .* distAlongBeams);
    nonMappedHeightAboveSensorBeam4 = (cos(beamAngle + -sg3*pitch) / cos(beamAngle) .* distAlongBeams);

I have confirmed this with the Scripps processing tools, trigonometry review, and all the literature I can find. However, when I compare this with the bin-mapped data output from the RDI Velocity processing software, the current method used in the toolbox agrees better. I'm going around in circles a bit.

@ggalibert and @ocehugo, maybe you can explain the maths behind the algorithm currently in the toolbox? Am I missing something here? Happy to review in a call if easier.

petejan commented

@BecCowley The multiplication by CP and CR is needed otherwise this is along the beam 1-2 plane in the case of CP, not the vertical (earth) distance from the sensor. Its not clear why its / cos (beamAngle) to me.

BecCowley commented

@petejan can we catch up and go through this please?

sspagnol commented

Hi @BecCowley, can you post the links for the Scripps tools. And how does all this relate to the maths in RDI doc below?

Appendix C - Coordinate XFORM Trig Functions.pdf

BecCowley commented

@sspagnol thanks for the document, I've not seen this before and I really needed it! I will review.
I've put the scripps tools on cloudstor here: https://cloudstor.aarnet.edu.au/plus/s/G0AeoN4Bu98EwYk
Bin-mapping and coordinate transformation is controlled in WHbeam_ProcessFunctionVectorized.m

BecCowley commented

After much investigation into trigonometry, experiments with RDI Velocity software, Scripps software and the toolbox, and with @petejan's help, here are my very brief conclusions:

  • The toolbox bin mapping algorithm is fine, except that it is only set up for down-facing instruments. Any up-facing instruments that have used this bin mapping PP tool should possibly be re-processed. I will submit a bug request to incorporate the bin-mapping for up-facing instruments.
  • RDI Velocity software and the toolbox have identical results with conversion to ENU from beam velocities, when no bin-mapping is included
  • When conversion to ENU is done with bin mapping, the bin mapping in the toolbox gives different results to the RDI Velocity software. I put this down to the method used (toolbox uses an interpolation, RDI uses the nearest cell).
  • The differences between the two bin-mapping methods are mostly within +/-0.2m/s in U and V (up-facing) and +/-0.1m/s (down-facing) in the test dataset I used (EAC0500, 2019-2021 data from upward 300khz and downward 75kHz instruments in single ping mode). Most of the larger offsets (up to a maximum of 0.5m/s) occurred in bins beyond the surface or bottom hits. Some occurred in bins where there is ringing or fish interference.

I can provide test datasets and figures if anyone is interested.
I couldn't code the RDI trig functions in the document that @sspagnol provided to replicate anything I got out of the Velocity software or the toolbox. But, I did have to adapt it a bit as it works from raw hertz and I have the beam velocities. Possibly I missed something....

Also, tried processing the test dataset with the Scripps functions, but again, had to hack it a bit to get it to output what I wanted with no QC filtering, and not gridded onto a standard depth grid. Possibly I didn't hack it properly....

I'm satisfied that the toolbox output matches very nearly what the RDI Velocity software is producing, once the bin mapping has been fixed to allow for up-facing instruments.

hugo-sardi commented

Hello, sorry for the delay, I was not aware of the discussion since I turned off my notifications.

maybe you can explain the maths behind the algorithm currently in the toolbox? Am I missing something here? Happy to review in a call if easier.

Can't say much - didn't touch that bit of code despite some small changes. I recall the equation is quite locked to supposedly upward-looking instruments and have been like that for ages and mostly untouched until v2.6.11+.

For the record, I did a small refactor in the routine since v2.6.10 given the request to implement all the RDI beam2ENU conversions and down looking adcps on your PR. There were regressions tests for these changes and the first unit tests were made available (see tests/PreProcessing), but they are far from complete.
I imagine it wouldn't hurt to test results between v2.6.10 and master to double-check if some of your trouble is actually recent bugs introduced - I may have missed something. If this is the case, it is quite likely a matter of interpolation optimization that was done for speed purposes.

Please also note that the binMapping equation code submitted on the PR (which was baked within the rdibeam2earth) was actually not performing any bin mapping (If you want a reference, check the comments on the test/Preprocessing routines). Moreover, it is/was identical to the one already presented in adcpBinMappingPP. Because of that, I went to invest time first in the beam2Earth conversion implementation details and verification, while this was left for later (as advised).

The toolbox bin mapping algorithm is fine, except that it is only set up for down-facing instruments. Any up-facing instruments that have used this bin mapping PP tool should possibly be re-processed. I will submit a bug request to incorporate the bin-mapping for up-facing instruments.

Did you mean downward facing!?

Finally, I would avoid creating such a code pattern:

if ~isUpwardLooking %for down facing adcps
    sg1 = 1; sg3 = 1;
else %for up looking
    sg1 = 1; sg3 = -1;
end

and would use a "typed in" formula for down/up cases. This way the transformation is explicit instead of implicit. Also, please note the comments in the files "% function block" is a hint that that part of the code should be a function instead. This way it is easy to test and verify.

BecCowley commented

@hugo-sardi thanks for your feedback. Yes, you are correct, I mean that any DOWN-FACING instruments would have to have re-processing considered.
The issue with the binmapping equation in the rdibeam2earth code has been noted and I will re-process these files (and all our single-ping ADCP files) to ensure they are fixed.

Re the niceties of coding, I will leave that to the IMOS coding team to tidy up.