lkilcher / dolfyn

A library for oceanographic doppler instruments such as Acoustic Doppler Profilers (ADPs, ADCPs) and Acoustic Doppler Velocimeters (ADVs).
BSD 3-Clause "New" or "Revised" License
43 stars 26 forks source link

Velocity magnitude changes when using 'Declination' #2

Closed lkilcher closed 9 years ago

lkilcher commented 9 years ago

The velocity magnitude changes when motion correction is applied even when motion signals are set to zero (dat.Accel[:] = 0; dat.AngRt[:] = 0). This SHOULD NOT be happening.

For example, using a data file from a 'Turbulence Torpedo' deployment, the following code:

import dolfyn.adv.api as avm
import matplotlib.pyplot as plt
import numpy as np

# Create a default motion correction instance
mc = avm.motion.CorrectMotion()

# Load the data file:
dat = avm.read_nortek('data_file.VEC')

# Specify a declination
# !!! This is crucial, the bug is minor or non-existant without this line.
dat.props['declination'] = 14.5

# Specify the postion/orientation of the ADV head.
dat.props['body2head_vec'] = np.array([0.21, 0, 0.603])  # in meters
dat.props['body2head_rotmat'] = np.array([[0, 0, 1],
                                          [0, 1, 0],
                                          [-1, 0, 0], ])
# Grab the relavent data
dat = dat.subset(slice(22150,103900))

# Copy the data for performing motion correction.
dat2 = dat.copy()
# Zero these out so that I can see the effect of 'orientation signals' only.
dat2.Accel[:] = 0
dat2.AngRt[:] = 0
mc(dat2)
avm.rotate.earth2principal(dat2)

#######
# Compare velocity magnitudes.
fig = plt.figure(2, figsize=[5, 5])
fig.clf()
ax = plt.subplot(211)
ax.plot(-dat.u)
ax.plot(-dat2.u)
ax.axhline(0, color='k', linestyle='--')
ax.set_xlim([0, 10000])
ax.set_ylim([0, 4])
ax.set_ylabel('Streamwise velocity [m/s]')
ax = plt.subplot(212)
ax.plot(np.sqrt((dat2._u ** 2).sum(0)) - np.sqrt((dat._u ** 2).sum(0)))
ax.axhline(0, color='k', linestyle='--', linewidth=2)
ax.set_xlim([0, 10000])
ax.set_ylim((-0.4, 0.8))
ax.set_ylabel('$|u_\mathrm{mc}| - |u_\mathrm{raw}|\,\mathrm{[m/s]}$', size='large')

Generates this figure,

bad_magnitude

So, the question is clearly, why is the magnitude of the velocity different (second panel is non-zero) when the motion signals (third panel) are zero?

Note here that this issue is not particularly large when the declination is not specified. Also note that this appears to be related to something going wrong with the orientation matrix:

In [61]: np.linalg.det(dat.orientmat[:,:,0])
Out[61]: 0.99999988

In [62]: np.linalg.det(dat2.orientmat[:,:,0])
Out[62]: -0.18341382

The determinant of the orientation matrix after motion-correction is non-zero?!! This is no good. Also, the orientation matrix is non-zero even when I comment out the dat.props['declination'] = 14.5 line. In that case, I get:

In [64]: np.linalg.det(dat2.orientmat[:,:,0])
Out[64]: -0.35748771

Why does this problem exist? Why does changing declination mess up the velocity magnitude, but the orientation matrix gets messed up - differently - either way.

lkilcher commented 9 years ago

This problem had to do with applying rotations only to a portion of the rotation matrix (.orientmat[:2,:2]). For some reason I messed up the linear-algebra here. Commit 59028fb51ad4275f688d0384530a3e0c2c4c2ed1 fixed this issue.

The bad-plot above now looks like this,

fixed_magnitude

Yay! :)

lkilcher commented 9 years ago

Also, as a side note, the reason declination was messing up velocity is because it was applied to .orientmat before data was rotated. The .orientmat was also getting messed up during the earth2principal call, but that was an 'after rotation' change to .orientmat, so it didn't mess up the data.