# Three-sat FAC density estimation with Swarm¶

Adrian Blagau (Institute for Space Sciences, Bucharest)
Joachim Vogt (Jacobs University Bremen)
Version Feb. 2020

This notebook accompanies the article “Multipoint Field-Aligned Current Estimates With Swarm” by A. Blagau, and J. Vogt, 2019. When used for publications, please acknowledge the authors’ work by citing the paper.

Introduction The notebook makes use of the three-s/c gradient estimation method, developed in Vogt, Albert, and Marghitu, 2009 to estimates the field-aligned current (FAC) density from the Swarm measurements when the satellites are in a close configuration. A detailed discussion on applying the three-s/c method in the context of Swarm (notably selection of events and other critical aspects) is provided in Blagau and Vogt, 2019

In the Input parameters section, the user specifies the interval of analysis and the magnetic model(s) used to compute the magnetic field perturbation. Optionally, the sensor data can be shifted in time (tshift array) and/or can be filtered.

Importing useful libraries (numpy, pandas, matplotlib, …)

# Uncomment if necessary:
# %matplotlib inline
# Uncomment for interactivity:
# %matplotlib widget

import numpy as np
import pandas as pd
from tqdm.auto import tqdm
from scipy import signal
import matplotlib.pyplot as plt
import datetime as dtm
import matplotlib.dates as mdt

import warnings
warnings.filterwarnings('ignore')


Preparing to access ESA’s Swarm mission data and models from VirES environment

from viresclient import SwarmRequest


## Defines some convenience functions¶

Computes the parameters of a low-pass Butterworth filter, if one decides to filter the data

fs = 1  # Sampling frequency
fc = 1./20.  # Cut-off frequency of a 20 s low-pass filter
w = fc / (fs / 2) # Normalize the frequency
butter_ord = 5
bf, af = signal.butter(butter_ord, w, 'low')


Defines recivec3s function to compute the reciprocal vectors for a three-s/c configuration

def recivec3s(R):
Rc = np.mean(R, axis=-2)
Rmeso = R - Rc[:, np.newaxis, :]
r12 = Rmeso[:,1,:] - Rmeso[:,0,:]
r13 = Rmeso[:,2,:] - Rmeso[:,0,:]
r23 = Rmeso[:,2,:] - Rmeso[:,1,:]
nuvec = np.cross(r12, r13)
nuvec_norm = np.linalg.norm(nuvec, axis=-1, keepdims=True)
nuvec = np.divide(nuvec, nuvec_norm)
Q3S = np.stack((np.cross(nuvec,r23), np.cross(nuvec,-r13), np.cross(nuvec,r12)),axis = -2)
Q3S = np.divide(Q3S, nuvec_norm[...,None])
Qtens = np.sum(np.matmul(Q3S[:,:,:,None],Q3S[:,:,None, :]), axis = -3)
Rtens = np.sum(np.matmul(Rmeso[:,:,:,None],Rmeso[:,:,None, :]), axis = -3)
return Rc, Rmeso, nuvec, Q3S, Qtens, Rtens


## Input parameters¶

Specifying the time interval, the satellites, the magnetic field model. Optionally, the sensor data can be shifted in time and/or can be filtered

dtime_beg = '2014-06-29T12:33:00'
dtime_end = '2014-06-29T12:38:00'
sats = ['A','B','C']
Bmodel="CHAOS-all='CHAOS-Core'+'CHAOS-Static'+'CHAOS-MMA-Primary'+'CHAOS-MMA-Secondary'"
#tshift = [0, 10, 4]  # optional time-shift introduced between the Swarm sensors
tshift = [0, 0, 0]
use_filter = False   # 'True' for filtering the data


## Data retrieval and preparation¶

Reads from VirES the sat. position (Rsph), magnetic L1b measurement (Bsw), and magnetic field model (Bmod). The quality flag of VFM data, i.e. Flags_B in the L1b files, is read as well. For a description of Flags_B values see Swarm Level 1b Product Definition from the Swarm documentation. The script identifies (and leaves unchanged) poor quality data points, i.e. when Flags_B > 0, as well as data gaps (and fills them with the values of the nearest neighboring point).

dti = pd.date_range(start = dtime_beg, end = dtime_end, freq='s', closed='left')
ndti = len(dti)
nsc = len(sats)
datagaps={}
Rsph = np.full((ndti,nsc,3),np.nan)
Bsw = np.full((ndti,nsc,3),np.nan)
Bmod = np.full((ndti,nsc,3),np.nan)
FlagsB = np.full((ndti,nsc),np.nan)
request = SwarmRequest()
for sc in tqdm(range(nsc)):
request.set_collection("SW_OPER_MAG"+sats[sc]+"_LR_1B")
request.set_products(measurements=["B_NEC","Flags_B"],
models=[Bmodel],
sampling_step="PT1S")
data = request.get_between(start_time = dtime_beg,
end_time = dtime_end,
asynchronous=False, show_progress=False)
print('Used MAG L1B file: ', data.sources)
dat = data.as_dataframe()
datagaps[sats[sc]] = dti.difference(dat.index)
dsi = dat.reindex(index=dti, method='nearest')
# store position, magnetic field and magnetic model vectors in corresponding data matrices
Bsw[:,sc,:] = np.stack(dsi['B_NEC'].values, axis=0)
Bmod[:,sc,:] = np.stack(dsi['B_NEC_CHAOS-all'].values, axis=0)
FlagsB[:,sc] = dsi['Flags_B'].values
# collect all data in a single DataFrame for inspection
names=['Var','Sat','Com'])
dfRsph = pd.DataFrame(Rsph.reshape(-1,nsc*3),columns=colRsph,index=dti)

colBswBmod = pd.MultiIndex.from_product([['Bsw','Bmod'],sats,['N','E','C']],
names=['Var','Sat','Com'])
dfBswBmod = pd.DataFrame(np.concatenate((Bsw.reshape(-1,nsc*3),
Bmod.reshape(-1,nsc*3)),axis=1),
columns=colBswBmod,index=dti)

dfFB = pd.DataFrame(FlagsB.reshape(-1,nsc),columns=sats,index=dti)

RsphBswBmod = pd.merge(dfRsph, dfBswBmod, left_index=True, right_index=True)

Used MAG L1B file:  SW_OPER_MAGA_LR_1B_20140629T000000_20140629T235959_0505_MDR_MAG_LR

Used MAG L1B file:  SW_OPER_MAGB_LR_1B_20140629T000000_20140629T235959_0505_MDR_MAG_LR

Used MAG L1B file:  SW_OPER_MAGC_LR_1B_20140629T000000_20140629T235959_0505_MDR_MAG_LR


Shows the first lines of data structure for inspection

RsphBswBmod.head()

Var Rsph Bsw Bmod
Sat A B C A ... C A B C
Com Lat Lon Radius Lat Lon Radius Lat Lon Radius N ... C N E C N E C N E C
2014-06-29 12:33:00 -80.911390 -65.774670 6851755.55 -80.169390 -61.467951 6900623.63 -81.554000 -63.031595 6851779.47 12732.0097 ... -37338.7848 12734.024209 6188.853291 -37157.458968 12941.555713 5199.868902 -35698.749635 12832.159674 5776.273701 -37342.619704
2014-06-29 12:33:01 -80.972259 -65.660006 6851756.68 -80.230700 -61.386827 6900625.20 -81.614416 -62.898266 6851780.59 12728.3808 ... -37365.0720 12730.191872 6174.737562 -37183.780781 12934.753470 5191.139742 -35728.940568 12829.870849 5757.982971 -37369.111138
2014-06-29 12:33:02 -81.033090 -65.543742 6851757.82 -80.291988 -61.304637 6900626.76 -81.674782 -62.762957 6851781.71 12724.2641 ... -37391.1771 12726.512653 6160.298355 -37210.044344 12928.027208 5182.188523 -35759.085743 12827.751426 5739.281490 -37395.546731
2014-06-29 12:33:03 -81.093880 -65.425848 6851758.95 -80.353253 -61.221361 6900628.32 -81.735100 -62.625625 6851782.82 12720.9187 ... -37417.2302 12722.989087 6145.529018 -37236.250076 12921.378343 5173.011078 -35789.184950 12825.803900 5720.160144 -37421.926831
2014-06-29 12:33:04 -81.154631 -65.306289 6851760.07 -80.414497 -61.136977 6900629.87 -81.795367 -62.486228 6851783.94 12717.7181 ... -37443.4390 12719.623649 6130.422607 -37262.397986 12914.808375 5163.603280 -35819.238413 12824.030626 5700.609473 -37448.251130

5 rows × 27 columns

Prints a report about the quality of VFM data, i.e. the sum of Flags_B values for each satellite (ideally 0) and the number of missing data points (ideally 0). Print timestamps with nonzero Flags_B values or when data is missing (if applicable).

FBnonzero = dfFB.loc[ (dfFB.A > 0) | (dfFB.B > 0) | (dfFB.C > 0)]
FBnA = dfFB['A'].loc[ (dfFB.A > 0)]
FBnB = dfFB['B'].loc[ (dfFB.B > 0)]
FBnC = dfFB['C'].loc[ (dfFB.C > 0)]
print('Number of points with non-zero quality flag / missing data:')
print('Swarm A: ', FBnA.count(), ' / ', datagaps['A'].size)
print('Swarm B: ', FBnB.count(), ' / ', datagaps['B'].size)
print('Swarm C: ', FBnC.count(), ' / ', datagaps['C'].size)
print()
if dfFB.sum(axis=0).sum(0) :
print('Time stamps with nonzero FlagsB:')
print(FBnonzero)
if datagaps['A'].size :
print('Missing data for Swarm A: ', datagaps['A'].values)
if datagaps['B'].size :
print('Missing data for Swarm B: ', datagaps['B'].values)
if datagaps['C'].size :
print('Missing data for Swarm C: ', datagaps['C'].values)

Number of points with non-zero quality flag / missing data:
Swarm A:  0  /  0
Swarm B:  0  /  0
Swarm C:  0  /  0


Computes sats positions (R), magnetic measurements (B), and magnetic field perturbations (dB) in the global geographic (Cartesian) frame. Takes care of the optional time-shifts introduced between the sensors. Filters the magnetic field perturbations if requested so.

ts = np.around(np.array(tshift) - min(tshift))
tsh3s = np.around(np.mean(ts), decimals=3)
ndt = ndti - max(ts)
dt = dti[:ndt].shift(1000.*tsh3s,freq='ms')  # new data timeline
R = np.full((ndt,nsc,3),np.nan)
B = np.full((ndt,nsc,3),np.nan)
dB = np.full((ndt,nsc,3),np.nan)
for sc in range(nsc):
# prepare conversion to global cartesian frame
clt,slt = np.cos(latsc.flat),np.sin(latsc.flat)
cln,sln = np.cos(lonsc.flat),np.sin(lonsc.flat)
north = np.stack((-slt*cln,-slt*sln,clt),axis=-1)
east = np.stack((-sln,cln,np.zeros(cln.shape)),axis=-1)
center = np.stack((-clt*cln,-clt*sln,-slt),axis=-1)
# store cartesian position vectors in position data matrix R
# store magnetic field measurements in the same frame
bnecsc = Bsw[ts[sc]:ts[sc]+ndt,sc,:]
B[:,sc,:] = np.matmul(np.stack((north,east,center),axis=-1),
bnecsc[...,None]).reshape(bnecsc.shape)
# store magnetic field perturbation in the same frame
dbnecsc = Bsw[ts[sc]:ts[sc]+ndt,sc,:] - Bmod[ts[sc]:ts[sc]+ndt,sc,:]
dB[:,sc,:] = np.matmul(np.stack((north,east,center),axis=-1),
dbnecsc[...,None]).reshape(dbnecsc.shape)
if use_filter:
dB[:,sc,:] = signal.filtfilt(bf, af, dB[:,sc,:], axis=0)
# collect all data in a single DataFrame
colRBdB = pd.MultiIndex.from_product([['R','B','dB'],sats,['x','y','z']],
names=['Var','Sat','Com'])
RBdB = pd.DataFrame(np.concatenate((R.reshape(-1,nsc*3),B.reshape(-1,nsc*3),dB.reshape(-1,nsc*3)),axis=1),
columns=colRBdB,index=dt)


Reads two-s/s FAC estimate from the L2 product for later comparison

request.set_collection('SW_OPER_FAC_TMS_2F')
request.set_products(measurements="FAC", sampling_step="PT1S")
data = request.get_between(start_time = dtime_beg,
end_time = dtime_end,
asynchronous=False)
print('Used FAC file: ', data.sources)
FAC_L2 = data.as_dataframe()

Used FAC file:  SW_OPER_FAC_TMS_2F_20140629T000000_20140629T235959_0301


## Computes the FAC density and other parameters¶

Uses recivec3s() function to compute the mesocenter of the Swarm constellation (Rc), the s/c positions in the mesocentric frame (Rmeso), the direction normal to spacecraft plane (nuvec), the reciprocal vectors (Q3S), the reciprocal tensor (Qtens), and the position tensor (Rtens).

Rc, Rmeso, nuvec, Q3S, Qtens, Rtens = recivec3s(R)


Computes the direction of (un-subtracted) local magnetic field Bunit and the orientation of spacecraft plane with respect to Bunit (cosBN and angBN). Register (for later exclusion) data points where B vector is too close to the spacecraft plane, i.e. below a threshold value set in angTHR.

Bunit = B.mean(axis=-2)/np.linalg.norm(B.mean(axis=-2),axis=-1)[...,None]
cosBN = np.matmul(Bunit[...,None,:],nuvec[...,:,None]).reshape(dt.shape)
angBN = pd.Series(np.arccos(cosBN)*180./np.pi, index=dt)
angTHR = 25.
texcl = (np.absolute(angBN) < 90 + angTHR) & (np.absolute(angBN) > 90 - angTHR)


Estimates the curl of magnetic field perturbation curlB. Associates its normal / radial component curlBrad with an electric current of density Jfac, flowing along the magnetic field line.

muo = 4*np.pi*1e-7
CurlB = np.sum( np.cross(Q3S,dB,axis=-1), axis=-2 )
Jfac.loc[texcl] = np.nan
plt.figure(figsize = [10,4])
plt.plot(Jfac)
plt.ylabel('$J_{FAC}$\n[$\mu A/m^2$]')
plt.xlabel('Time')

Text(0.5, 0, 'Time') The condition number CN3 evaluates the stability of the solution (the quality of the Rtens inversion). The error errJ in FAC density due to (mutually uncorrelated and isotropic) instrumental noise is evaluated by considering a flat value db nT for δB

eigval = np.sort(np.linalg.eigvals(Rtens), axis=-1)
CN3 = pd.Series(np.log10(np.divide(eigval[:,2],eigval[:,1])),index=dt)

db = 0.5
if use_filter:
db = 0.2
traceq = np.trace(Qtens, axis1=-2, axis2=-1)
errJ = 1e-6*db/muo*pd.Series(np.sqrt(traceq)/np.absolute(cosBN),index=dt)
errJ.loc[texcl] = np.nan


## Plotting and saving the results¶

The script below produces a plot and two ASCII files, for the input data and for the results. A third ASCII file (input data quality report) is produced if non-zero quality flag and/or missing data points are present in the input data. The generated plot presents:

• the NEC components of the magnetic field perturbation, dB, for all three satellites. A common range along y axes is used

• the logarithm of the condition number, CN3

• angle between the normal to the spacecraft plane and the direction of local ambient magnetic field.

• comparison between the FAC density estimated with the three-s/c method and the ESA L2 FAC product obtained from the dual-s/c method.

• the FAC estimation errors due to instrumental noise

The values for tick labels refer to the position of the mesocenter.

On the bottom, the spacecraft configurations are shown at three instances (i.e. beginning, middle, and end of the interval) projected on the NE plane of the local NEC coordinate frame tied to the mesocenter. The s/c velocity vectors are also indicated by arrows.

%run -i "plot_and_save_three_sat.py" 