{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"# Single-sat FAC estimation with Swarm"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"Adrian Blagau (Institute for Space Sciences, Bucharest) \n",
"Joachim Vogt (Jacobs University Bremen) \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 single-s/c method to compute the field-aligned current (FAC) and ionospheric radial current (IRC) densities on Swarm. The algorithm offers some advantages over the one used to generate the L2 product since (i) both low (LR) and high resolution (HR) L1b magnetic field data can be used, (ii) input data can be filtered by the user, and (iii) the inclination of FAC sheet can be taken into account provided that this information is known (e.g. as a result of applying the Minimum Variance Analysis).\n",
"\n",
"While the algorithm is able to handle HR data, the presence of data gap (not marked in L1b files) could lead to unphysical results both for un-filtered and filtered analysis. In the LR files the missing data are marked by zero values on all magnetic field components; when such points occur within the analysis interval, the un-filtered analysis replaces them by NaN, prints the corresponding timestamps and prevents the analysis on filtered data. \n",
"\n",
"As input parameters (see the corresponding section), the user specifies the interval of analysis, the satellite, and the solicited resolution of magnetic field data. Optionally, the current sheet inclination can be specified in several ways (see below). The algorithm relies on CHAOS magnetic model(s) to compute the magnetic field perturbation, but the user can select another model available on the VirES platform. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Importing useful libraries (numpy, pandas, matplotlib, ...)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Uncomment if necessary:\n",
"# %matplotlib inline\n",
"# Uncomment for interactivity:\n",
"# %matplotlib widget\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"from scipy import signal\n",
"import matplotlib.pyplot as plt\n",
"import datetime as dtm\n",
"import matplotlib.dates as mdt\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Preparing to access ESA’s Swarm mission data and models from VirES environment"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from viresclient import SwarmRequest"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Definition of convenience functions "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def normvec(v):\n",
" # Given an array of vectors v, the function returns \n",
" # the corresponding array of unit vectors\n",
" return np.divide(v,np.linalg.norm(v,axis=-1).reshape(-1,1))\n",
"\n",
"def rotvecax(v, ax, ang):\n",
" # Rotates vector v by angle ang around a normal vector ax \n",
" # Uses Rodrigues' formula when v is normal to ax\n",
" sa, ca = np.sin(np.deg2rad(ang)), np.cos(np.deg2rad(ang))\n",
" return v*ca[...,np.newaxis] + np.cross(ax, v)*sa[...,np.newaxis]\n",
"\n",
"def sign_ang(V, N, R):\n",
" # returns the signed angle between vectors V and N, perpendicular \n",
" # to R; positive sign corresponds to right hand rule along R\n",
" VxN = np.cross(V, N)\n",
" pm = np.sign(np.sum(R*VxN, axis=-1))\n",
" return np.degrees(np.arcsin(pm*np.linalg.norm(VxN, axis=-1))) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function *R_B_dB_in_GEOC* provides the (Cartesian) position vector *R*, the magnetic field *B*, and the magnetic field perturbation *dB* in the GEOC frame."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def R_B_dB_in_GEOC(Rsph, Bnec, dBnec):\n",
" latsc = np.deg2rad(Rsph[:,0])\n",
" lonsc = np.deg2rad(Rsph[:,1]) \n",
" radsc = 0.001*Rsph[:,2]\n",
" # prepares conversion to global cartesian frame\n",
" clt,slt = np.cos(latsc.flat),np.sin(latsc.flat)\n",
" cln,sln = np.cos(lonsc.flat),np.sin(lonsc.flat)\n",
" north = np.stack((-slt*cln,-slt*sln,clt),axis=-1)\n",
" east = np.stack((-sln,cln,np.zeros(cln.shape)),axis=-1)\n",
" center = np.stack((-clt*cln,-clt*sln,-slt),axis=-1)\n",
" # stores cartesian position vectors in position data matrix R\n",
" R = -radsc[...,None]*center\n",
" # stores magnetic data in GEOC (same frame as for R)\n",
" Bgeo = np.matmul(np.stack((north,east,center),axis=-1),\n",
" Bnec[...,None]).reshape(Bnec.shape)\n",
" dBgeo = np.matmul(np.stack((north,east,center),axis=-1),\n",
" dBnec[...,None]).reshape(dBnec.shape)\n",
" return R, Bgeo, dBgeo"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function *singleJfac* computes FAC and IRC densities, together with the corresponding error estimations, based on the single satellite method. Mandatory parameters are the satellite position *R*, magnetic field *B*, magnetic field perturbation *dB* (all arrays of vectors in GEOC), together with the corresponding timestamps *t*. \\\n",
"The current sheet inclination (optional) can be specified with one of the parameters *alpha* (angle of inclination in the tangential plane wrt satellite velocity vector), *N2d* (projection of sheet normal on the tangential plane), or *N3d* (sheet normal in GEOC). For details see Input parameters section. The time-interval (optional) where the inclination is valid can be provided as well; if not, the information on inclination is assumed valid for the whole interval of analysis. \\\n",
"The error in FAC and IRC densities are estimated based on the value of *er_db*, that specifies the error in magnetic field perturbation. The function also returns an array with FAC inclination in the tangential plane."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def singleJfac(t, R, B, dB, alpha=None, N2d=None, \\\n",
" N3d=None, tincl=None, er_db=0.5):\n",
" # Constructs the differences & values at mid-intervals\n",
" dt = t[1:].values - t[:-1].values\n",
" tmid = t[:-1].values + dt*0.5\n",
" Bmid = 0.5*(B[1:,:] + B[:-1,:]) \n",
" Rmid = 0.5*(R[1:,:] + R[:-1,:])\n",
" diff_dB = dB[1:,:] - dB[:-1,:] \n",
" V3d = R[1:,:] - R[:-1,:]\n",
" Vorb = np.sqrt(np.sum(V3d*V3d, axis=-1)) \n",
" # Defines important unit vectors\n",
" eV3d, eBmid, eRmid = normvec(V3d), normvec(Bmid), normvec(Rmid)\n",
" eV2d = normvec(np.cross(eRmid, np.cross(eV3d, eRmid))) \n",
" # Angle between B and R\n",
" cos_b_r = np.sum(eBmid*eRmid, axis=-1)\n",
" bad_ang = np.abs(cos_b_r) < np.cos(np.deg2rad(60))\n",
" \n",
" # incl is the array of FAC incliation wrt Vsat (in tangential plane) \n",
" if N3d is not None:\n",
" eN3d = normvec(N3d)\n",
" eN2d = normvec(eN3d - np.sum(eN3d*eRmid,axis=-1).reshape(-1,1)*eRmid)\n",
" incl = sign_ang(eV2d, eN2d, eRmid)\n",
" elif alpha is not None:\n",
" incl = alpha if isinstance(alpha, np.ndarray) else \\\n",
" np.full(len(tmid), alpha) \n",
" elif N2d is not None:\n",
" eN2d = normvec(np.cross(eRmid, np.cross(N2d, eRmid)))\n",
" incl = sign_ang(eV2d, eN2d, eRmid)\n",
" else:\n",
" incl = np.zeros(len(tmid))\n",
"\n",
" # considers the validity interval of FAC inclination \n",
" if tincl is not None:\n",
" ind_incl = np.where((tmid >= tincl[0]) & (tmid <= tincl[1]))[0]\n",
" incl[0:ind_incl[0]] = incl[ind_incl[0]]\n",
" incl[ind_incl[-1]:] = incl[ind_incl[-1]]\n",
"\n",
" # working in the tangential plane\n",
" eNtang = normvec(rotvecax(eV2d, eRmid, incl))\n",
" eEtang = normvec(np.cross(eNtang, eRmid))\n",
" diff_dB_Etang = np.sum(diff_dB*eEtang, axis=-1)\n",
" Dplane = np.sum(eNtang*eV2d, axis=-1)\n",
" j_rad= - diff_dB_Etang/Dplane/Vorb/(4*np.pi*1e-7)*1.e-6\n",
" j_rad_er= np.abs(er_db/Dplane/Vorb/(4*np.pi*1e-7)*1.e-6) \n",
" \n",
" # FAC density and error\n",
" j_b = j_rad/cos_b_r\n",
" j_b_er = np.abs(j_rad_er/cos_b_r) \n",
" j_b[bad_ang] = np.nan\n",
" j_b_er[bad_ang] = np.nan \n",
" \n",
" return tmid, Rmid, j_b, j_rad, j_b_er, j_rad_er, incl, np.arccos(cos_b_r)*180./np.pi"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function *GapsAsNaN* sets magnetic data gaps to NaN"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def GapsAsNaN(df_ini, ind_gaps):\n",
" df_out = df_ini.copy()\n",
" df_out['B_NEC'][ind_gaps] = [np.full(3,np.NAN)]*len(ind_gaps)\n",
" return df_out, df_ini.index[ind_gaps]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function *rez_param* provides the sampling step (needed in SwarmRequest) and the data sampling frequency (needed for data filtering) according to the magnetic data resolution"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def rez_param(rez):\n",
" sstep = 'PT1S' if rez=='LR' else 'PT0.019S' # sampling step\n",
" fs = 1 if rez=='LR' else 50 # data sampling freq.\n",
" return sstep, fs"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"\n",
"## Input parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Specifying the time interval, satellite, data resolution, and magnetic field model.\\\n",
"Optionally, the current sheet inclination can be specified by providing just one of the following parameters:\n",
"- *alpha* is the angle (in degree) of FAC inclination in the tangential plane wrt satellite velocity vector. The tangential plane is perpendicular to the satellite position vector. *alpha* can be a single value or series of values at mid-intervals, positive or negative according to the right-hand rotation along the position vector. Implicit value: None\n",
"- *N3d* is the FAC sheet normal in GEOC (three components vector). This is the usual output from MVA. Implicit value: None\n",
"- *N2d* is the projection of FAC sheet normal in the tangential plane (i.e. the GEOC components). It can be a vector or series of vectors at mid-intervals. Implicit value: None\n",
"\n",
"The time interval when FAC inclination is to be considered can be specified in *tincl* (optional); outside this interval, the FAC inclination is assumed to take the values at the beginning/end of *tincl*. When *tincl* is not specified, the information on inclination is assumed valid for the whole interval of analysis "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dtime_beg = '2015-03-17T08:51:00'\n",
"dtime_end = '2015-03-17T08:58:00.1'\n",
"\n",
"sat = ['A']\n",
"\n",
"alpha, N3d, N2d, tincl = None, None, None, None\n",
"## Optional: FAC inclination\n",
"# alpha = -20. \n",
"# N3d = [-0.32780841, 0.82644295, -0.45774851]\n",
"# N2d = [-0.326, 0.828, -0.457]\n",
"# tincl = np.array(['2015-03-17T08:51:54','2015-03-17T08:57:11'],dtype='datetime64')\n",
"\n",
"rez = 'LR' # 'LR' or 'HR' for low or high resolution data\n",
"use_filter = True # 'True' for filtering the data\n",
"\n",
"Bmodel=\"CHAOS-all='CHAOS-Core'+'CHAOS-Static'+'CHAOS-MMA-Primary'+'CHAOS-MMA-Secondary'\""
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Data retrieval and preparation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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. Computes the magnetic perturbation (in NEC) and filters it for later use. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"request = SwarmRequest()\n",
"request.set_collection(\"SW_OPER_MAG\"+sat[0]+\"_\"+rez+\"_1B\")\n",
"request.set_products(measurements=[\"B_NEC\",\"Flags_B\"], \n",
" auxiliaries=['QDLat','QDLon','MLT'],\n",
" models=[Bmodel],\n",
" sampling_step=rez_param(rez)[0])\n",
"data = request.get_between(start_time = dtime_beg, \n",
" end_time = dtime_end,\n",
" asynchronous=False) \n",
"print('Used MAG L1B file: ', data.sources[1])\n",
"dat_df = data.as_dataframe()\n",
"\n",
"# sets missing B_NEC data (zero magnitude in L1b LR files) to NaN. \n",
"# imposes 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",
"nti = len(ti)\n",
"# stores position, magnetic field and magnetic model vectors in corresponding data matrices\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",
"FlagsB = dat_df['Flags_B'].values\n",
"\n",
"if use_filter:\n",
" fc, butter_ord = 1/20, 5 # 20 s cutt-off freq., filter order\n",
" bf, af = signal.butter(butter_ord, fc /(rez_param(rez)[1]/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. Compute FAC, IRC densities, and the corresponding estimation errors. 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, jrad, jb_er, jrad_er, incl, ang_BR = \\\n",
" singleJfac(ti, R, B, dB, alpha=alpha, N2d=N2d, N3d=N3d, tincl=tincl)\n",
"j_df = pd.DataFrame(np.stack((Rmid[:,0], Rmid[:,1], Rmid[:,2], \\\n",
" jb, jrad, jb_er, jrad_er, ang_BR, incl)).transpose(),\\\n",
" columns=['Rmid X','Rmid_Y','Rmid Z','FAC','IRC','FAC_er','IRC_er','ang_BR','incl'], index=tt)\n",
"dB_df = pd.DataFrame(dB, columns=['dB Xgeo', 'dB Ygeo', 'dB Zgeo'], index=ti)\n",
"dB_nec_df = pd.DataFrame(dBnec, columns=['dB N', 'dB E', 'dB C'], index=ti)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Computes filtered FAC, IRC densities if data filtering is desired and possible. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if use_filter:\n",
" R, B, dB_flt = R_B_dB_in_GEOC(Rsph, Bnec, dBnec_flt)\n",
" tt, Rmid, jb_flt, jrad_flt, jb_er_flt, jrad_er_flt, incl, ang_BR = \\\n",
" singleJfac(ti, R, B, dB_flt, alpha=alpha, N2d=N2d, N3d=N3d, tincl=tincl, er_db=0.2)\n",
" jflt_df = pd.DataFrame(np.stack((Rmid[:,0], Rmid[:,1], Rmid[:,2],\\\n",
" jb_flt, jrad_flt, jb_er_flt, jrad_er_flt, ang_BR, incl)).transpose(),\\\n",
" columns=['Rmid X','Rmid_Y','Rmid Z','FAC_flt','IRC_flt',\\\n",
" 'FAC_flt_er','IRC_flt_er','ang_BR','incl'], index=tt) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Reads the single-s/c FAC estimate from the L2 product for comparison"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"request.set_collection('SW_OPER_FAC'+sat[0]+'TMS_2F')\n",
"request.set_products(measurements=[\"FAC\",\"IRC\"], sampling_step=\"PT1S\")\n",
"data = request.get_between(start_time = dtime_beg, \n",
" end_time = dtime_end,\n",
" asynchronous=False) \n",
"print('Used FAC file: ', data.sources[0])\n",
"FAC_L2 = data.as_dataframe()\n",
"FAC_L2.rename(columns={'FAC':\"FAC_L2\", 'IRC':\"IRC_L2\"}, inplace = True)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting and saving the results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plots and saves the current density data as ASCII file. Panel 1: magnetic field perturbation in GEOC, panel 2: un-filtered and (when applicable) filtered FAC, panel 3: un-filtered and (when applicable) filtered IRC, panel 4: comparison with FAC L2, panel 5: angle between *B* and *R* vectors, panel 6: considered angle between FAC normal and satellite velocity in the tangential plane. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%run -i \"plot_and_save_single_sat.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
}