Minimum Variance Analysis on the Magnetic Field¶

Here, I try to reproduce the paper, Romaneli et al., 2016, JGR, "Proton cyclotron waves occurrence rate upstream fromMars observed by MAVEN: Associated variability of the Martian upper atmosphere".

Reference: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016JA023270 https://www.researchgate.net/publication/234549390_Minimum_Variance_Analysis_Cluster_Themes

In [31]:
%config InlineBackend.figure_format='retina'
import sys
from datetime import datetime
from rich.console import Console
console = Console(highlight=False)#,force_jupyter=False,width=999)


from cdasws.cdasws import CdasWs, TimeInterval
from cdasws.datarepresentation import DataRepresentation as dr
from cdasws.doi import get_doi_landing_page_url, get_doi_badge_url, get_doi_hdp_url

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

from IPython.core.display import HTML

today = datetime.now().strftime('%Y-%B-%d')

cdas = CdasWs()

def display_link(
    url: str,
    text: str) -> None:
    display(HTML(f'<a href="{url}" target="_blank">{text}</a>'))

def display_image_link(
    url: str,
    image_url: str) -> None:
    display(HTML(f'<a href="{url}" target="_blank"><img src={image_url}/></a>'))
    

# Injects CSS to target the notebook's standard output stream container
display(HTML('''
<style>
div.output_area pre {
    font-size: 16px !important;
}
</style>
'''))
In [2]:
cdas = CdasWs()
In [3]:
# Fetch all observatory groups
obs_groups = cdas.get_observatory_groups()

# Print the names
print("Observatory -> IDs\n")
for indx, group in enumerate(obs_groups):
    #print("'%s'" %(group['Name']),',', end="")
    print("\033[1;32m%s\033[0m -> %s" %(group['Name'], group['ObservatoryId']))
    #print(f'{group["Name"]}')#"'%s'" %(group['Name']),',', end="")
Observatory -> IDs

ACE -> ['AC', 'ACE', 'OMNI (1AU IP Data)']
AMPTE -> ['AMPTE-CCE']
ARTEMIS -> ['THB', 'THC']
AWE -> ['AWE']
Alouette -> ['Alouette-2']
Apollo -> ['ALSEP']
Arase (ERG) -> ['Arase(ERG)', 'ERG']
Balloons -> ['PMC Turbo', 'bar_1A', 'bar_1B', 'bar_1C', 'bar_1D', 'bar_1G', 'bar_1H', 'bar_1I', 'bar_1J', 'bar_1K', 'bar_1M', 'bar_1N', 'bar_1O', 'bar_1Q', 'bar_1R', 'bar_1S', 'bar_1T', 'bar_1U', 'bar_1V', 'bar_2A', 'bar_2B', 'bar_2C', 'bar_2D', 'bar_2E', 'bar_2F', 'bar_2I', 'bar_2K', 'bar_2L', 'bar_2M', 'bar_2N', 'bar_2O', 'bar_2P', 'bar_2Q', 'bar_2T', 'bar_2W', 'bar_2X', 'bar_2Y', 'bar_3A', 'bar_3B', 'bar_3C', 'bar_3D', 'bar_3E', 'bar_3F', 'bar_3G', 'bar_4A', 'bar_4B', 'bar_4C', 'bar_4D', 'bar_4E', 'bar_4F', 'bar_4G', 'bar_4H', 'bar_5A', 'bar_6A', 'bar_7A']
BioSentinel -> ['BioSentinel']
CNOFS -> ['CNOFS']
CRRES -> ['CRRES Satellite']
Cassini -> ['CASSINI', 'Cassini']
Cluster -> ['C1', 'C2', 'C3', 'C4', 'CL', 'CT']
DMSP -> ['DMSP', 'dmsp', 'dmsp-f12', 'dmsp-f13', 'dmsp-f14', 'dmsp-f15', 'dmsp-f16', 'dmsp-f17', 'dmsp-f18', 'dmsp-f6', 'dmsp-f7', 'dmsp-f8', 'dmsp-f9']
DSCOVR -> ['DSCOVR']
Demonstration and Science Experiments (DSX) -> ['DSX']
Dynamics Explorer -> ['DE', 'DE-1', 'DE-2', 'DE1']
Equator-S -> ['EQ']
FAST -> ['FAST']
Formosat -> ['FS-5', 'ROCSAT-1(FORMOSAT-1)']
GOES -> ['G0', 'G10', 'G12', 'GOES-08', 'GOES-09', 'GOES-10', 'GOES-11', 'GOES-12', 'GOES-13', 'GOES-14', 'GOES-15', 'GOES-16', 'GOES-17', 'GOES-18', 'GOES-19', 'GOES10', 'GOES11', 'GOES12', 'GOES13', 'GOES14', 'GOES15', 'GOES16', 'GOES17', 'GOES18', 'GOES8', 'GOES9', 'GOES_10', 'GOES_11', 'GOES_12', 'GOES_6', 'GOES_7', 'GOES_8', 'GOES_9', 'g11']
GOLD -> ['GOLD']
GPS -> ['GPS']
Genesis -> ['Genesis']
Geosynchronous Investigations -> ['STPSat-6']
Geotail -> ['GEOTAIL']
Ground-Based Investigations -> ['CANOPUS', 'DARN', 'FMI', 'SESAME', 'THG_L2']
Hawkeye -> ['HK']
Helio Ephemeris -> ['BEPICOLOMBO', 'BORRELLY', 'COMETGS', 'COMETHMP', 'DAWN', 'EARTH', 'GALILEO', 'GIACOBINI', 'GIOTTO', 'HALEBOPP', 'HALLEY', 'HYAKUTAKE', 'JUICE', 'JUNO', 'JUPITER', 'MARS', 'MERCURY', 'MOON', 'MSL', 'NEPTUNE', 'PHOBOS2', 'PLUTO', 'PSYCHE', 'ROSETTA', 'SAKIGAKE', 'SATURN', 'SOLAR-1', 'SUISEI', 'URANUS', 'VENUS', 'WILD2']
Helios -> ['HELIOS', 'HELIOS-1', 'HELIOS-2', 'HELIOS1', 'HELIOS2', 'Helios-1', 'Helios-2']
IBEX -> ['IBEX']
ICON -> ['ICON', 'ICON ']
IMAGE -> ['IM', 'IMAGE']
IMAP -> ['IMAP']
IMP (All) -> ['I8', 'IMP-8', 'IMP7', 'IMP8', 'OMNI (1AU IP Data)']
ISEE -> ['ISEE 1', 'ISEE-1', 'ISEE-3', 'ISEE1', 'ISEE2']
ISIS -> ['ISIS-1', 'ISIS-2']
ISS -> ['AWE', 'ISS']
Interball -> ['INTERBALL-AURORAL', 'INTERBALL-GROUND', 'INTERBALL-TAIL']
LANL -> ['L1', 'LANL-97A', 'LANL1989_046', 'LANL1990_095', 'LANL1991_080', 'LANL1994_084', 'Lanl-01A', 'Lanl-02A', 'Lanl-97A', 'Lanl1989_046', 'Lanl1990_095', 'Lanl1991_080', 'Lanl1994_084']
MAVEN -> ['MAVEN']
MESSENGER -> ['MESSENGER']
MMS -> ['MMS1', 'MMS2', 'MMS3', 'MMS4', 'mms1', 'mms2', 'mms3', 'mms4']
Mariner -> ['MARINER2', 'Mariner2']
Mars Global Surveyor (MGS) -> ['MGS']
Mars Science Laboratory (MSL) -> ['MSL']
Moon -> ['BGM1']
NOAA -> ['MetOp1', 'MetOp2', 'MetOp3', 'NOAA05', 'NOAA07', 'NOAA08', 'NOAA10', 'NOAA12', 'NOAA14', 'NOAA15', 'NOAA16', 'NOAA18', 'NOAA19']
New Horizons -> ['NEW_HORIZONS', 'NH']
OMNI (Combined 1AU IP Data; Magnetic and Solar Indices) -> ['OMNI (1AU IP Data)']
POES/MetOp -> ['MetOp1', 'MetOp2', 'MetOp3', 'NOAA15', 'NOAA16', 'NOAA18', 'NOAA19']
Parker Solar Probe (PSP) -> ['PSP', 'PSP_FLD']
Pioneer -> ['PIONEER10', 'PIONEER11', 'Pioneer10', 'Pioneer11', 'Pioneer6', 'Pioneer7', 'PioneerVenus']
Polar -> ['PO', 'POLAR', 'Polar', 'po']
REACH -> ['REACH', 'reach']
SAMPEX -> ['SAMPEX']
SOHO -> ['SOHO']
SOLAR -> ['SOLAR-1']
ST5 -> ['ST5']
STEREO -> ['Ahead', 'Behind', 'STA', 'STB', 'STEREO', 'STEREOA', 'STEREOB', 'sta', 'stb']
Smallsats/Cubesats -> ['CIRBE', 'CSSWE', 'Munin', 'SNOE', 'ac6a', 'ac6b', 'ela', 'elb']
Solar Orbiter -> ['SOLO', 'solo']
Sounding Rockets -> ['Dynamo-2', 'Endurance', 'RENU2', 'TOOWINDY']
THEMIS -> ['THA', 'THB', 'THC', 'THD', 'THE', 'THEMIS', 'THG', 'THG_L1', 'THG_L2']
TIMED -> ['TIMED']
TRACERS -> ['TRACERS1', 'TRACERS2', 'TS1', 'TS2']
TSS-1R -> ['TSS-1R']
TWINS -> ['TWINS-', 'TWINS1', 'TWINS2']
Ulysses -> ['ULYSSES', 'Ulysses']
Van Allen Probes (RBSP) -> ['RBSP', 'RBSP-A', 'RBSP-A-RBSPICE', 'RBSP-B', 'RBSP-B-RBSPICE', 'RBSPA', 'RBSPB', 'rbsp-a', 'rbsp-b', 'rbspa', 'rbspb']
Voyager -> ['VG1', 'VG2', 'VGR-1', 'VGR-2', 'VOYAGER1', 'VOYAGER2', 'Voyager 1', 'Voyager 2', 'Voyager1', 'Voyager2']
Wind -> ['OMNI (1AU IP Data)', 'WI', 'WIND', 'Wi', 'Wind', 'wi']
In [4]:
# Fetch all instrument types
instrument_types = cdas.get_instrument_types()

# print the names
for group in instrument_types:
    print("'%s'" %(group['Name']))
' Particles (space)'
'Activity Indices'
'Dust and Debris'
'Electric Fields (space)'
'Engineering'
'Ephemeris/Attitude/Ancillary'
'Ground-Based HF-Radars'
'Ground-Based Imagers'
'Ground-Based Magnetometers, Riometers, Sounders'
'Ground-Based VLF/ELF/ULF, Photometers'
'Housekeeping'
'Imagers (space)'
'Imaging and Remote Sensing (ITM/Earth)'
'Imaging and Remote Sensing (Magnetosphere/Earth)'
'Imaging and Remote Sensing (Sun)'
'Magnetic Fields (Balloon)'
'Magnetic Fields (space)'
'Particle Flux (space)'
'Particles (space)'
'Plasma and Solar Wind'
'Radio and Plasma Waves (space)'
'Space Weather Monitor'
In [5]:
### user-specified ###
myobservatory=['MAVEN']
myinstrument =['Magnetic Fields (space)']
In [6]:
all_datasets = []
for obs in myobservatory:
    print('for %s' %(obs))
    print('Dataset with the link -> Label')
    datasets = cdas.get_datasets(observatoryGroup=obs,instrumentType=myinstrument)
    if datasets:
        all_datasets.extend(datasets)
        
    lines=[]
    for index, dataset in enumerate(datasets):
        dataset_id = dataset["Id"]
        dataset_label = dataset["Label"]
        dataset_note = dataset['Notes']
        lines.append("[bold green][link=%s]%s[/link][/] -> %s" %(dataset_note,dataset_id,dataset_label))
    console.print("\n".join(lines),justify='left', soft_wrap=True)

    print('='*100)
    print(' ')
for MAVEN
Dataset with the link -> Label
MVN_INSITU_KP-4SEC -> MAVEN In-situ Key Parameters - B.M. Jakosky (LASP/U. Colorado)                               
MVN_MAG_L2-SUNSTATE-1SEC -> MAVEN Magnetometer Sun-State 1 Second Magnetic Field  - J. Connerney (NASA GSFC)       
====================================================================================================
 
In [9]:
### user-specified ###
mydataset=['MVN_MAG_L2-SUNSTATE-1SEC'] #the order should be consistent with the list above.
In [11]:
mydois=[]
for index, dataset in enumerate(all_datasets):
    dataset_id = dataset["Id"]
    dataset_label = dataset["Label"]

    if dataset_id in mydataset:  # select this dataset for use below
        doi=dataset['Doi']
        #print(doi)
        mydois.append(dataset['Doi'])
        print(f'\nSelected dataset: {dataset_id} (doi:{doi})')
        print(f"Available Period: {dataset['TimeInterval']}")
        display_image_link(get_doi_landing_page_url(doi), get_doi_badge_url(doi))
        display_link(get_doi_hdp_url(doi), 'HDP Description')
        display_link(dataset['Notes'], 'CDAWeb Notes')
        display(HTML(f'<b>Data Citation</b>: {cdas.get_citation(doi)}.  Accessed on {today}.'))
        
        print('='*100)
        print(' ')
Selected dataset: MVN_MAG_L2-SUNSTATE-1SEC (doi:10.48322/b9da-ph25)
Available Period: {'Start': '2014-10-10T00:00:00.000Z', 'End': '2025-12-06T00:00:00.000Z'}
HDP Description
CDAWeb Notes
Data Citation: Connerney, J. E. P. (2023). MAVEN Magnetometer (MAG) Magnetic Field and Orbital Position, Sun-State and Payload Coordinates, Level 2 (L2), 1 s Data [Dataset]. Space Physics Data Facility. https://doi.org/10.48322/B9DA-PH25. Accessed on 2026-September-23.
====================================================================================================
 
In [12]:
for indx,mydoi in enumerate(mydois):
    variables = cdas.get_variables(mydoi)
    print('for %s' %(mydataset[indx]))
    
    search_window = TimeInterval("1900-01-01T00:00:00Z", "2027-01-01T00:00:00Z")
    filtered_inventory = cdas.get_inventory(mydataset[indx], timeInterval=search_window)

    for interval in filtered_inventory:
        print(f"Available from: \033[1;32m{interval.start}\033[0m to \033[1;32m{interval.end}\033[0m")
    print(' ')
    var_names = []
    for index, variable in enumerate(variables):
        name = variable['Name']
        var_names.append(name)
        print("\033[1m'%s'\033[0m" %(variable['Name']),"->", variable['LongDescription'])
    #    print(f'{index}. {name}, {variable["LongDescription"]}')
    print('='*100)
    print(' ')
for MVN_MAG_L2-SUNSTATE-1SEC
Available from: 2014-10-10 00:00:00+00:00 to 2014-10-16 10:44:23+00:00
Available from: 2014-10-17 00:00:00+00:00 to 2014-10-25 14:42:16+00:00
Available from: 2014-10-29 16:26:03+00:00 to 2014-11-19 03:37:50+00:00
Available from: 2014-11-27 05:16:59+00:00 to 2014-12-26 19:05:19+00:00
Available from: 2014-12-27 01:03:50+00:00 to 2015-03-03 18:51:39+00:00
Available from: 2015-03-04 01:24:38+00:00 to 2015-04-03 14:16:47+00:00
Available from: 2015-04-14 17:05:19+00:00 to 2015-06-29 15:45:53+00:00
Available from: 2015-07-01 20:49:03+00:00 to 2015-11-20 07:46:54+00:00
Available from: 2015-11-21 00:00:01+00:00 to 2016-07-13 00:00:01+00:00
Available from: 2016-07-13 23:59:59+00:00 to 2018-02-28 05:08:25+00:00
Available from: 2018-03-01 16:28:29+00:00 to 2019-09-04 12:23:39+00:00
Available from: 2019-09-12 22:49:24+00:00 to 2019-09-19 00:00:01+00:00
Available from: 2019-09-19 23:59:59+00:00 to 2020-06-07 06:35:25+00:00
Available from: 2020-06-08 17:14:16+00:00 to 2020-06-08 17:20:35+00:00
Available from: 2020-06-09 21:05:08+00:00 to 2020-11-01 00:00:01+00:00
Available from: 2020-11-15 00:00:02+00:00 to 2021-01-06 05:06:16+00:00
Available from: 2021-01-07 07:06:01+00:00 to 2021-01-09 02:28:55+00:00
Available from: 2021-01-10 08:59:10+00:00 to 2021-01-11 22:45:14+00:00
Available from: 2021-01-14 00:00:00+00:00 to 2021-01-14 21:25:31+00:00
Available from: 2021-01-15 19:09:46+00:00 to 2021-02-22 21:07:30+00:00
Available from: 2021-02-23 06:48:35+00:00 to 2021-02-23 22:01:14+00:00
Available from: 2021-02-24 16:26:32+00:00 to 2021-11-01 00:00:00+00:00
Available from: 2021-11-15 00:00:01+00:00 to 2021-11-23 21:01:46+00:00
Available from: 2021-11-24 10:00:31+00:00 to 2022-02-22 15:18:25+00:00
Available from: 2022-05-15 00:00:01+00:00 to 2022-09-15 07:01:47+00:00
Available from: 2022-09-16 00:15:46+00:00 to 2023-02-16 00:00:01+00:00
Available from: 2023-02-21 19:08:24+00:00 to 2023-11-01 00:00:00+00:00
Available from: 2023-11-15 00:00:11+00:00 to 2024-02-01 00:00:00+00:00
Available from: 2024-02-15 00:00:02+00:00 to 2024-05-01 00:00:01+00:00
Available from: 2024-05-15 00:00:01+00:00 to 2024-05-31 23:59:59+00:00
Available from: 2024-06-01 23:59:59+00:00 to 2024-07-01 00:00:00+00:00
Available from: 2024-08-15 00:00:00+00:00 to 2025-10-30 00:00:01+00:00
Available from: 2025-10-30 23:59:59+00:00 to 2025-11-15 23:59:59+00:00
Available from: 2025-11-16 23:59:59+00:00 to 2025-12-06 00:00:00+00:00
 
'OB_B' -> OUTBOARD B J2000
'OB_B_range' -> ---> OUTBOARD B J2000 range
'POSN' -> SC position
'OB_BDPL' -> OUTBOARD BD PAYLOAD
'OB_BDPL_range' -> ---> OUTBOARD BD PAYLOAD range
====================================================================================================
 
In [61]:
# assign parameteprs
#dataset_id = 'AC_H0_SWE'
myvariables  = {}
mystart_time = {}
myend_time   = {}

### user-specified ###
myvariables[mydataset[0]]=['OB_B']
mystart_time[mydataset[0]] = '2014-12-24T08:00:00Z'
myend_time[mydataset[0]]   = '2014-12-24T09:00:00Z'
#######################

# request the data representation as Spacepy
data={}
for dataset in mydataset:
    print(dataset)
    status, data[dataset] = cdas.get_data(dataset, myvariables[dataset], mystart_time[dataset], myend_time[dataset], dataRepresentation = dr.XARRAY)
    if status['http']['status_code'] == 200:
        print('Donwload Success!')
    else:
        print(f"Error! HTTP Status Code: {status}")

for dataset in mydataset:
    keys=list(data[dataset].indexes.keys())
    data[dataset] = data[dataset].rename({keys[0]: "time"})
MVN_MAG_L2-SUNSTATE-1SEC
Donwload Success!
In [98]:
start = datetime(2014,12,24,8,55,24)
end   = datetime(2014,12,24,8,56, 4)

fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(16, 8))
ax[0].plot(data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B']['time'][:], data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B'][:,0],'C0')
ax[1].plot(data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B']['time'][:], data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B'][:,1],'C1')
ax[2].plot(data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B']['time'][:], data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B'][:,2],'C2')
ax[0].set_xlim('2014-12-24T08:50:00Z','2014-12-24T09:00:00Z')
ax[0].set_ylim(3,6)
ax[1].set_ylim(0,3)
ax[2].set_ylim(0,3)
ax[0].set_ylabel(r'$B_x$ [nT]')
ax[1].set_ylabel(r'$B_y$ [nT]')
ax[2].set_ylabel(r'$B_z$ [nT]')

### add a rectangle patch ###
import matplotlib.dates as mdates
import matplotlib.patches as patches
startplt = mdates.date2num(start)
endplt   = mdates.date2num(end)
width = endplt - startplt
rect = patches.Rectangle((startplt, 0), width, 3+2*fig.subplotpars.hspace, color='grey',alpha=.25, clip_on=False,transform=ax[2].get_xaxis_transform())
#rect = plt.Rectangle((start,0), width, height=12, clip_on=False, edgecolor="none", facecolor="blue", linewidth=3,alpha=0.25)
ax[2].add_patch(rect)
#############################

plt.show()

Apply the MVA method to the magnetic field data¶

The purpose of this method is to find a vector that makes the following variance minimum (or maximum), given 3D magnetic field data, $$\sigma^2=\frac{1}{N}\sum_k^N({\bf B}_k\cdot{\bf x}-<{\bf B}_k\cdot{\bf x}>)^2,$$ where $k$ is the data point and $<\cdot>$ is a mean.

To find the vector ${\bf x}$ that minises the variance, the easiest way is probably to solve the following equation, $${\bf Mx}=\lambda{\bf x},$$ where ${\bf M}=<B_iB_j>-<B_i><B_j>$, $\lambda$ is the eigenvalue, and ${\bf x}$ is the eigenvector.

In [105]:
bxyz=data['MVN_MAG_L2-SUNSTATE-1SEC']['OB_B'][:,:].sel(time=slice(start, end)).values

Mb=np.zeros((3,3))

for i in range(3):
    for j in range(3):
        Mb[i,j]=np.mean(bxyz[:,i]*bxyz[:,j])-np.mean(bxyz[:,i])*np.mean(bxyz[:,j])

v,w=np.linalg.eig(Mb)

sorted_indices = np.argsort(v)

v_sort=v[sorted_indices]
w_sort=w[:,sorted_indices]

print(w_sort)

### define new IMFs based on the MVA ###
bnew1=np.dot(bxyz[:,:],w_sort[:,2])
bnew2=np.dot(bxyz[:,:],w_sort[:,1])
bnew3=np.dot(bxyz[:,:],w_sort[:,0])


### display values ###
from IPython.display import Math
display(Math(r'\lambda_1=%.4f,\: \lambda_2=%.4f,\: \lambda_3=%.4f' %(v_sort[2],v_sort[1],v_sort[0])))
display(Math(r'\lambda_1/\lambda_2=%.4f,\: \lambda_2/\lambda_3=%.4f' %(v_sort[2]/v_sort[1],v_sort[1]/v_sort[0])))
display(Math(r'{\bf B}_0=[%.4f, %.4f, %.4f]' %(np.mean(bnew1),np.mean(bnew2),np.mean(bnew3))))

### plot figures ###
plt.figure(figsize=(12,4))

plt.subplot(121)
plt.plot(bnew2[0],bnew1[0],'o')
plt.plot(bnew2[2],bnew1[2],'ro')
plt.plot(bnew2,bnew1,'.-')
plt.xlabel(r'$B_2$ [nT]')
plt.ylabel(r'$B_1$ [nT]')
#plt.xlim(-0.25,1.75)
#plt.ylim(-1.7,0.3)

plt.subplot(122)
plt.plot(bnew3,bnew2)
plt.xlim(-5.7,-3.7)
plt.xlabel(r'$B_3$ [nT]')
plt.ylabel(r'$B_2$ [nT]')
#plt.ylim(-0.25,1.75)

plt.show()
[[-0.84572647 -0.50876732 -0.16094272]
 [-0.49636364  0.86076553 -0.11272021]
 [-0.1958823   0.01544436  0.98050579]]
$\displaystyle \lambda_1=0.4099,\: \lambda_2=0.3083,\: \lambda_3=0.0071$
$\displaystyle \lambda_1/\lambda_2=1.3296,\: \lambda_2/\lambda_3=43.1508$
$\displaystyle {\bf B}_0=[0.7508, -0.7515, -4.7323]$
In [ ]: