OpenDrift / opendrift

Open source framework for ocean trajectory modelling
https://opendrift.github.io
GNU General Public License v2.0
231 stars 113 forks source link

Possible bug in LCS computation ? #1269

Closed simonweppe closed 2 months ago

simonweppe commented 2 months ago

Hi @knutfrode , @johannesro

Opening up a PR for review/discussion on this.

In the FTLE computation, one needs to calculate the displacement from each initial position (X,Y), that are then fed to ftle().

It seems that opendrift inverts the order of (lon,lat) when run in backwards mode. So doing this below, would return b_x1,b_y1 matrices that are not ordered the same way as X and Y.

b_x1, b_y1 = proj(o.history['lon'].T[-1].reshape(X.shape),o.history['lat'].T[-1].reshape(X.shape))

This means that when doing this b_x1-X it's not substracting the correct initial positions from the final positions. Flipping the final (lon,lat) arrays seems to yield a consistent matrix ordering, and so correct net displacements fed to ftle

b_x1, b_y1 = proj(o.history['lon'].T[-1][::-1].reshape(X.shape),o.history['lat'].T[-1][::-1].reshape(X.shape))

Doing this suppresses the need to flip the array after the ftle() call; there was a note in the code that it was unclear why it was needed.

Now the confusing part...I ran the example example_double_gyre_LCS.py with and without these modifications ..and results were actually quite consistent. That's probably a reason why this was not identified at the time...but I'm a bit surprised that results are that consistent with different displacement inputs...Keen to hear your insights on this.

simonweppe commented 2 months ago

Forgot to mention that it's only needed for the ALCS case that uses backwards simulations. Ordering is consistent for RLCS/forward simulations.

knutfrode commented 2 months ago

Hi Simon,

I have personally barely used the FTLE functionality after Johannes and myself implemented it long time ago. Thus I would trust that your new suggestion is correct, but agree that it is a bit strange that it does not look much different from current implementation.

Thus unless @johannesro or @mateuszmatu (or others) have some other insights to share, I would be fine with merging this (apparent) correction.

Btw: presently there is some hardcoded FTLE-plotting functionality baked into the plot/animation methods. In the near future, I plan to remove this, and rather use Xarray in combination with more generic methods to achieve the same, as in this example https://opendrift.github.io/gallery/example_river_runoff.html

simonweppe commented 2 months ago

Thanks Knut-Frode. Yes will be good to have insights from @johannesro and @mateuszmatu or any that are using it more often.

Ok noted about the improved xarray links with opendrift, that's really useful.

knutfrode commented 2 months ago

I did not get any feedback, but merging this now, as it looks fine to me.