{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Minimum Variance Analysis on Swarm\n",
"\n",
"Adrian Blagau (Institute for Space Sciences, Bucharest) <br>\n",
"Joachim Vogt (Jacobs University Bremen) <br>\n",
"Version May 2021"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook accompanies the article [\"Multipoint Field-Aligned Current Estimates With Swarm\"](https://doi.org/10.1029/2018JA026439) by A. Blagau, and J. Vogt, 2019. When used for publications, please acknowledge the authors' work by citing the paper."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Introduction** The notebook implements the algorithm of Minimum Variance Analysis (MVA) to estimate the planarity and orientation of Field Aligned Current (FAC) structures probed by the Swarm satellites. The FAC orientation can be later used to improve the FAC density estimated through the standard implementation of single-sat. method, that assumes only normal (i.e. oriented along the satellite velocity) current sheets. \n",
"\n",
"The rough interval of analysis is provided by the used at the beginning; the notebook allows to inspect the data and, if needed, to refine the initial time interval. MVA is performed on the magnetic field perturbation projected on the plane perpendicular to the local average magnetic field direction. The ratio of the two eigenvalues provides an indication on the current sheet planarity while the eigenvector associated to the direction of minimum magnetic variance is associated with the current sheet orientation. \n",
"\n",
"The algorithm was designed to work with low-resolution (LR) L1b data. In the LR files the missing data are marked by zero values on all magnetic field components; such points are excluded from the analysis, their value replaced by NaN, and the corresponding timetags printed on the screen. \n",
"\n",
"The MVA results, provided as ASCII file, include the interval of analysis, the eigenvalues with the corresponding eigenvectors in GEOC frame, the angle of FAC inclination wrt sat. velocity measured in the tangential plane, and the components of magnetic field perturbation in the proper (i.e. current sheet) frame. The evolution of magnetic field perturbation, FAC density and inclination, as well as the hodograph of magnetic field perturbation in the plane perpendicular to the average magnetic field are\n",
"plotted on screen and saved as eps file."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importing useful libraries/ functions\n",
"\n",
"The standard python libraries (numpy, pandas, matplotlib, ...) are imported. Also, to easily follow the logical flow in the notebook, the functions listed below are read from the accompanying funcs_single_sat_MVA.py python script:\n",
"- *normvec* : computes the unit-vectors of a vector time series\n",
"- *rotvecax* : rotates a vector by an angle around another vector\n",
"- *sign_ang* : returns the signed angle between two vectors\n",
"- *R_B_dB_in_GEOC* : transforms in GEO Cartesian frame the satellite position (initially in spherical coordinates) and magnetic field & perturbation (initially in NEC) \n",
"- *singleJfac* : computes the single-satellite FAC density"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Uncomment if necessary:\n",
"# %matplotlib inline\n",
"# Uncomment for interactivity:\n",
"# (required for selecting analysis interval)\n",
"# %matplotlib widget\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"from scipy import signal\n",
"from scipy.interpolate import interp1d\n",
"import datetime as dtm\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.dates as mdt\n",
"from matplotlib.widgets import SpanSelector\n",
"from viresclient import SwarmRequest\n",
"\n",
"from funcs_fac import *\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Definitions of convenience functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Function to perform the MVA on an array of 3D vectors. To do the analysis in a plane (e.g. perpendicular to the local average magnetic field), the normal to that plane *cdir* (3D vector) should be specified. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def mva(v, cdir = None):\n",
" # returns results of MVA on the array of 3D vector v. If cdir (3D vector) \n",
" # is provided, the analysis is performed in the plane perpendiculat to it\n",
" v_cov = np.cov(v, rowvar=False, bias=True)\n",
" if cdir is not None:\n",
" ccol = cdir.reshape(3,1)\n",
" cunit = ccol / np.linalg.norm(ccol)\n",
" d_mat = (np.identity(3) - np.dot(cunit,cunit.T))\n",
" v_cov = np.matmul(np.matmul(d_mat, v_cov), d_mat)\n",
" return np.linalg.eigh(v_cov)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Function to manually select the analysis interval by dragging the mouse pointer on the screen."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def onselect(xmin, xmax):\n",
" nrl = len(ax2.lines)\n",
" ax2.lines[nrl-1].remove()\n",
" ax2.lines[nrl-2].remove()\n",
" span_sel['beg'], span_sel['end'] = xmin, xmax\n",
" ax2.axvline(xmin, ls='--', c='r')\n",
" ax2.axvline(xmax, ls='--', c='r')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='input'></a>\n",
"## Input parameters\n",
"As input parameters, the user specifies the time-interval and the satellite. To allow for a refined analysis, the algorithm actually downloads Swarm data recorded on a larger, i.e. half-orbit time span, and offers the possibility to modify the analysis interval specified below. To compute the magnetic field perturbation, the algorithm relies on CHAOS magnetic model(s), but the user can select another model available on the VirES platform. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Time range, satellites, and magnetic model in format '%Y-%m-%dT%H:%M:%S'\n",
"\n",
"dtime_beg = '2015-03-17T08:51:54'\n",
"dtime_end = '2015-03-17T08:57:11'\n",
"\n",
"sat = ['C']\n",
"\n",
"use_filter = True # 'True' for filtering the data\n",
"Bmodel = \"CHAOS-all='CHAOS-Core'+'CHAOS-Static'+'CHAOS-MMA-Primary'+'CHAOS-MMA-Secondary'\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data retrieval and preparation\n",
"\n",
"Reads from VirES the sat. position (*Rsph*), magnetic L1b measurement (*Bnec*), and magnetic field model (*Bmod*). Auxiliary parameters *QDLat*, *QDLon*, and *MLT*, used when plotting the results, are retrieved as well. Time instances when the magnetic measurements assume zero values (indicative of data gaps in L1b low-resolution files) are replace by NaN. Computes the magnetic perturbation (in NEC) and, if no data gaps are present, applies a low-pass filter. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"request = SwarmRequest()\n",
"# identify the right half-orbit interval \n",
"orb = request.get_orbit_number(sat[0], dtime_beg, mission='Swarm')\n",
"orb_beg, orb_end = request.get_times_for_orbits(orb, orb, mission='Swarm', spacecraft=sat[0])\n",
"half_torb = (orb_end - orb_beg)/2.\n",
"dtm_beg = dtm.datetime.fromisoformat(dtime_beg)\n",
"dtm_end = dtm.datetime.fromisoformat(dtime_end)\n",
"if dtm_beg - orb_beg < half_torb:\n",
" large_beg, large_end = orb_beg, orb_beg + half_torb\n",
"else:\n",
" large_beg, large_end = orb_beg + half_torb, orb_end\n",
"\n",
"# download the Swarm data\n",
"request.set_collection(\"SW_OPER_MAG\"+sat[0]+\"_LR_1B\")\n",
"request.set_products(measurements=[\"B_NEC\"], \n",
" auxiliaries=['QDLat','QDLon','MLT'],\n",
" models=[Bmodel],\n",
" sampling_step='PT1S')\n",
"dat = request.get_between(start_time = large_beg, \n",
" end_time = large_end,\n",
" asynchronous=False) \n",
"print('Used MAG L1B file: ', dat.sources[1])\n",
"dat_df = dat.as_dataframe()\n",
"\n",
"# set B_NEC data gaps (zero magnitude in L1b LR files) to NaN. \n",
"# impose no filtering if there are missing data points.\n",
"ind_gaps = np.where(\\\n",
" np.linalg.norm(np.stack(dat_df['B_NEC'].values), axis = 1)==0)[0]\n",
"if len(ind_gaps):\n",
" dat_df, timegaps = GapsAsNaN(dat_df, ind_gaps)\n",
" print('NR. OF MISSING DATA POINTS: ', len(ind_gaps))\n",
" print(timegaps.values)\n",
" print('NO FILTERING IS PERFORMED')\n",
" use_filter = False\n",
"\n",
"ti = dat_df.index\n",
"\n",
"# store position, magnetic field and magnetic model vectors as data array\n",
"Rsph = dat_df[['Latitude','Longitude','Radius']].values\n",
"Bnec = np.stack(dat_df['B_NEC'].values, axis=0)\n",
"Bmod = np.stack(dat_df['B_NEC_CHAOS-all'].values, axis=0) \n",
"dBnec = Bnec - Bmod\n",
"\n",
"# if possible, filter the magnetic field perturbation\n",
"if use_filter:\n",
" fc, fs = 1./20, 1 # Cut-off and sampling frequency in s\n",
" bf, af = signal.butter(5, fc / (fs / 2), 'low')\n",
" dBnec_flt = signal.filtfilt(bf, af, dBnec, axis=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Computes the (Cartesian) position vector *R*, magnetic field *B*, and magnetic field perturbation *dB* in the GEOC frame. The un-filtered and, when possible, filtered current density, i.e. *jb* and respectively *jb_flt*, are computed to better assess the selection of MVA interval. Store quantities in DataFrame structures"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"R, B, dB = R_B_dB_in_GEOC(Rsph, Bnec, dBnec)\n",
"tt, Rmid, jb,_,_,_,_,_ = singleJfac(ti, R, B, dB)\n",
"jb_df = pd.DataFrame(jb, columns=['FAC'], index=tt)\n",
"\n",
"if use_filter:\n",
" R, B, dB_flt = R_B_dB_in_GEOC(Rsph, Bnec, dBnec_flt)\n",
" tt, Rmid, jb_flt,_,_,_,_,_ = singleJfac(ti, R, B, dB_flt)\n",
" jb_flt_df = pd.DataFrame(jb_flt, columns=['FAC_flt'], index=tt)\n",
"\n",
"# eV2d is the satellite velocity in the tangent plane (i.e. perpendicular to position vector)\n",
"V3d = R[1:,:] - R[:-1,:] \n",
"eV3d, eRmid = normvec(V3d), normvec(Rmid)\n",
"eV2d = normvec(np.cross(eRmid, np.cross(eV3d, eRmid)))\n",
"\n",
"dBnec_df = pd.DataFrame(dBnec, columns=['dB N', 'dB E', 'dB C'], index=ti)\n",
"dB_df = pd.DataFrame(dB, columns=['dB Xgeo', 'dB Ygeo', 'dB Zgeo'], index=ti)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Selection of analysis interval\n",
"\n",
"The initial interval of analysis is indicated on a plot that presents the magnetic field perturbation and current density (when possible, the filtered current density is preferred). If the user wants to inspect the data he/she can zoom/pan both on x and y axes after clicking the ✥ button of the interactive window (displayed at the bottom). To select another interval for MVA, click again on the ✥ button, press left mouse button and drag to select a region in the top panel. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tmarg = dtm.timedelta(seconds=45)\n",
"fig, (ax1, ax2) = plt.subplots(2, figsize=(7, 3), sharex='all')\n",
"ax1.plot(dB_df)\n",
"ax1.set_xlim(xmin = dtm_beg - tmarg, xmax = dtm_end + tmarg)\n",
"ax1.axvline(pd.to_datetime(dtime_beg), ls='--', c='k')\n",
"ax1.axvline(pd.to_datetime(dtime_end), ls='--', c='k')\n",
"ax1.axhline(0, ls='--', c='k')\n",
"if use_filter:\n",
" ax2.plot(jb_flt_df)\n",
"else:\n",
" ax2.plot(jb_df)\n",
"ax2.axhline(0, ls='--', c='k')\n",
"ax2.axvline(pd.to_datetime(dtime_beg), ls='--', c='k')\n",
"ax2.axvline(pd.to_datetime(dtime_end), ls='--', c='k')\n",
"span = SpanSelector(ax1, onselect, 'horizontal', useblit=True, \n",
" span_stays=True, rectprops=dict(alpha=0.2, facecolor='red'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Minimum Variance Analysis\n",
"\n",
"Updates the analysis interval if needed and applies the MVA in the plane perpendicular to the average magnetic field direction. The current sheet normal is provided by the eigenvector associated with the minimum magnetic variance, while the sense of direction is chosen according to the sat. velocity"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# update MVA interval if span has been selected\n",
"tmva_int = np.array([dtime_beg, dtime_end], dtype='datetime64')\n",
"if any([s != 0 for s in span.extents]):\n",
" tmva_int = np.array(mdt.num2date(span.extents), dtype='datetime64')\n",
"# select quantities for MVA interval and remove NaN points\n",
"indok = np.where((ti >= tmva_int[0]) & (ti <= tmva_int[1]) & ~np.isnan(dB[:,0]))[0]\n",
"dB_int, B_int = dB[indok,:], B[indok,:]\n",
"B_ave = np.average(B_int, axis = 0)\n",
"B_unit = B_ave / np.linalg.norm(B_ave)\n",
"# apply constrained MVA\n",
"eigval,eigvec = mva(dB_int, cdir=B_unit)\n",
"# select the minvar orientation according to sat. velocity\n",
"eV3d_ave = np.average(eV3d[indok[:-1]], axis = 0) \n",
"mindir = eigvec[:,1] \n",
"if np.sum(mindir*eV3d_ave) < 0:\n",
" mindir = -eigvec[:,1]\n",
"maxdir = np.cross(B_unit, mindir)\n",
"\n",
"print('MVA interval: ', tmva_int)\n",
"print('B_unit= %10.2f'%eigval[0], np.round(B_unit, decimals=4))\n",
"print('mindir= %10.2f'%eigval[1], np.round(mindir, decimals=4))\n",
"print('maxdir= %10.2f'%eigval[2], np.round(maxdir, decimals=4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Transforms the magnetic perturbation in the MVA (proper) frame. It also calculates the FAC inclination wrt sat. velocity in the tangential plane and corrects the FAC density estimation accordingly. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# transform magnetic perturbation in MVA frame\n",
"geo2mva = np.stack((B_unit, mindir, maxdir), axis=1)\n",
"dBmva_df= pd.DataFrame(np.matmul(dB_df.values, geo2mva), \n",
" columns=['dB B', 'dB min', 'dB max'], index=ti)\n",
"\n",
"# compute the FAC inclination wrt sat. velocity in the tangential plane\n",
"eN2d, ang = eV2d.copy(), np.zeros(len(tt))\n",
"eN2d[indok[:-1]] = \\\n",
" normvec(np.cross(eRmid[indok[:-1]], np.cross(mindir, eRmid[indok[:-1]])))\n",
"\n",
"cross_v_n = np.cross(eV2d[indok[:-1]], eN2d[indok[:-1]])\n",
"sign_ang = np.sign(np.sum(eRmid[indok[:-1]]*cross_v_n, axis=-1))\n",
"ang[indok[:-1]] = \\\n",
" np.degrees(np.arcsin(sign_ang*np.linalg.norm(cross_v_n, axis=-1)))\n",
"ang[0:indok[0]] = ang[indok[0]]\n",
"ang[indok[-1]:] = ang[indok[-2]]\n",
"\n",
"# DataFrames with FAC inclination and density corrected for inclination\n",
"ang_df = pd.DataFrame(ang, columns=['ang_v_n'], index=tt)\n",
"tt, Bmid, jb_inc, _, _, _, _, _ = singleJfac(ti, R, B, dB, alpha=ang)\n",
"jb_inc_df = pd.DataFrame(jb_inc, columns=['FAC_inc'], index=tt)\n",
"if use_filter:\n",
" tt, Bmid, jb_inc_flt, _, _, _, _, _ = singleJfac(ti, R, B, dB_flt, alpha=ang)\n",
" jb_inc_flt_df = pd.DataFrame(jb_inc_flt, columns=['FAC_inc_flt'], index=tt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting and saving the MVA results \n",
"\n",
"The plot presents (i) the magnetic field perturbation in NEC, (ii) magnetic field perturbation in the proper frame, (iii) a comparison between the standard and corrected for inclination FAC density estimates (filtered quantities are preferred), and (iv) the angle (in degree) of FAC inclination in the tangential plane wrt satellite velocity vector. The hodograph of magnetic field perturbation in the plane perpendicular to the average magnetic field is presented as well.\n",
"The ASCII file contains the MVA results and the magnetic field perturbation in the proper (MVA associated) frame."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%run -i \"plot_and_save_single_sat_MVA.py\""
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}