Unexpected behavior when converting vectors

[jpreistad]

I am not sure if this turns out to be a bug, but I just wanted to show some plots I just made, that intuitively seems a bit suspicious to me. I have attached a csv datafile of some SuperDARN gridded line-of-sight velocity data that can be read with the below script to reproduce the plot that is also attached here as a pdf.

ocbpy_test.pdf

ocbpy_test.csv.txt

This apparently weird behavior might be to a wrong use of the library, so I would be very grateful for any guidance if that is the case. My plot show to the left the gridded superdarn vectors plotted on the MLT/MLAT grid. Overlaid is the AMPERE circle fit, where the statistical correction to OCB has been applied (Burrell et al 2020). In this case I think the blue line matches pretty well with structures also seen in the convection. However, when representing this in the OCB frame (right panel), this seems to disorganize things. That is sort of the opposite of what was the intention of using this in the first place. Vectors that appear inside the OCB to the left are outside the OCB to the right, and vice versa. In addition, some vectors show a very large rotation between the two frames as well, which seems unrealistic.

I think I am doing something wrong here, but can not find the source of this discrepancy. My intention is to make statistics of the convection, normalized to the OCB, but I can not proceed before I can sort this out. I will be grateful for any advice here!

I am running the master branch of ocbpy, installed using pip install ocbpy on september 10, 2020 on Ubuntu 20 using Anaconda 3 (python 3.7.6).

Sample code to reproduce plot:

import ocbpy
import pandas as pd
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt

def offset (mlt): #Empirical formula from Burrell et al 2020 Anngeo
    o = 4.01 * (1. - 0.55**2) / (1 + 0.55 * np.cos((2.*np.pi-np.radians(mlt*15.)) + 0.92))
    return o

#Global parameters
ocboffset = 3.34 #mean difference from Ampere radius to OCB
boundary = 74.34 # Based on typical size of oval for selected data interval
time = dt.datetime(2011,11,27,3,54)

#Load SuperDARN data
data = pd.read_csv('ocbpy_test.csv')
#OCBpy stuff
ocb = ocbpy.ocboundary.OCBoundary(instrument='ampere', stime=time, 
                                  etime=time + dt.timedelta(seconds=60), boundary_lat = boundary)
ocb.get_next_good_ocb_ind()
o_lat, o_mlt, r_corr = ocb.normal_coord(data.mlat, data.mlt)
dat_ind = np.arange(0, len(data.mlat))
pgr_vect = ocbpy.ocb_scaling.VectorData(dat_ind, ocb.rec_ind,
                                        data.mlat, data.mlt,
                        aacgm_n=data.vlos_north, aacgm_e=data.vlos_east,
                        aacgm_mag=np.sqrt(data.vlos_north**2 + data.vlos_east**2),
                        dat_name='LoS Velocity',
                        dat_units='m s$^{-1}$',
                        scale_func=ocbpy.ocb_scaling.normal_evar)
pgr_vect.set_ocb(ocb)

#Plotting
fig = plt.figure()
ax = fig.add_subplot(121, polar=True) #Left panel just the gridded SD data and the AMPERE circle fit
theta = np.radians(np.array(data.mlt)*15)
r = np.array(90-data.mlat)
dr = np.array(data.vlos_east)
dt = np.array(data.vlos_north)
ax.quiver(theta, r, dr * np.cos(theta) - dt * np.sin(theta), dr * np.sin(theta) + dt * np.cos(theta))
ax.set_theta_zero_location("S") #This make midnight on bottom, dawn to the right
#plot the Ampere circle fit and adding the statistical correction from Burrell.
x0 = ocb.x[0]
y0 = ocb.y[0]
r = ocb.r[0]
x = np.linspace(-25,25,1000)
#obey the equation (x-x0)**2 + (y-y0)**2 = r**2
y1 = (2*y0 + np.sqrt(4*r**2 - 4*(x-x0)**2))/2
y2 = (2*y0 - np.sqrt(4*r**2 - 4*(x-x0)**2))/2
lat1 = 90 - np.sqrt(x**2 + y1**2)
lat2 = 90 - np.sqrt(x**2 + y2**2)
lat = np.concatenate([lat1,lat2])
mlt1 = np.arctan2(y1, x)*12/np.pi + 6
mlt2 = np.arctan2(y2, x)*12/np.pi + 6
mlt = np.concatenate([mlt1,mlt2])
ax.plot(np.radians(mlt*15), (90-(lat+offset(mlt))))
ax.scatter(np.arctan2(y0,x0)+np.pi/2, np.sqrt(x0**2 + y0**2), s=100)

#Right panel the converted measurements and a circle representing boundary_lat
ax2 = fig.add_subplot(122, polar=True)
theta = np.radians(np.array(o_mlt)*15)
r = np.array(90-o_lat)
dr = np.array(pgr_vect.ocb_e)
dt = np.array(pgr_vect.ocb_n)
ax2.quiver(theta, r, dr * np.cos(theta) - dt * np.sin(theta), dr * np.sin(theta) + dt * np.cos(theta))
ax2.set_theta_zero_location("S")
ax2.plot(np.radians(np.linspace(0,24,100)*15), np.ones(100)*(90-boundary))

[aburrell]

Ok, this does look like a problem that could be a bug. I'll carve out some time to take a look at it this week. Thanks for the comprehensive example!

[jpreistad]

I have tried to identify what causes this behavior but not succeeded. It is very difficult to follow all the quadrant logic. However, it appear to be related to how the ocb east/north components obtain their sign. First, I thought it had to do with the 'calc_ocb_vec_sign' function. However, whats going on at lines 803-806 in ocb_scaling.py, (what I interpret as a projection of vmag along the ocb east and north directions), will also affect the sign (ocb_naz can be in [-180,180] I think?). Altogether this unfortunately lead to the wrong sing in the end in some instances. I really hoped that I would be able to get to the bottom of this, but I have failed to do that so far.

[aburrell]

Ok, so I am working on this now. The first step was to reproduce your figure using the file you supplied for the AACGM coordinates. There are some differences, namely the OCB is not oriented correctly in your AACGMV2 plot. However, there are still vectors that look suspicious in the plot on the right, so I think that difference is not particularly important to this issue. I will start by looking into the two vectors near the OCB just after midnight and then examine the vectors near the OCB pole.

[jpreistad]

Yes, I see that I have implemented the AMPERE OCB correction incorrectly, thank you for pointing that out! I was actually surprised in my plot how well the boundary matched the observed velocity shear around 6 MLT, but that was unfortunately not the case.

However, I really hope that you will be able to figure out what is causing this behavior. It seems like you reproduce the problem pointed out.

[aburrell]

Yes, I'm going through step by step. It looks like the OCB and vect quadrants are being assigned correctly, so now I'm looking at the 6 sign combinations. When this is all done, I'll may write up an explanation of this logic in the docs.

[jpreistad]

I should say that I have tried digging into this myself, a few weeks ago. I regret that I did not document my efforts, but I recall that I also concluded that the quadrant assignment was ok. I think I got stuck in understanding what was happening in the calc_ocb_vector_sign function. This may be a good place to check if the logic is correctly implemented :)

Good luck debugging, I really appreciate your efforts!

[aburrell]

I can confirm that the problem is in the calc_ocb_vector_sign function. A few of the points are getting the incorrect north sign. Wish me luck!

[aburrell]

Ok, I found the bug that caused the obvious error, but I want to go through all of the permutations again and check that I haven't made any more mistakes.

[jpreistad]

That is great, good job!

[aburrell]

Wow. Ok. 16 diagrams later, another double check, and I am as sure as one person can be that have caught all of the errors in the sign logic. However, this means that I am pretty sure there's an error in the OCB vector azimuth angle logic that was (in some cases) made invisible by the sign logic error. Check out the data near 06:00 MLT in the OCB dial:

Hopefully it won't take much longer to go through that logic (Lines 867-904 in ocb_scaling) now that I have all the diagrams drawn, but I may not be able to work on this next week.

[aburrell]

And... this kept bugging me so here we are. Check out the PR that addresses the bug: #88

[jpreistad]

Wow, I really appreciate your efforts! Looking forward to see how this impacts my results, I hope to be able to dig into this again this week.