Wavelet Analysis¶

reference: https://iopscience.iop.org/article/10.1088/0004-637X/751/1/19
https://paos.colorado.edu/research/wavelets/bams_79_01_0061.pdf

In [1]:
%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
import math as mt
import scipy.fft as sf

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]:
# 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 [3]:
# 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 [4]:
### user-specified ###
myobservatory=['Wind']
myinstrument =['Magnetic Fields (space)']
In [5]:
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']
        #dataset_note = dataset['DatasetLink']['Url']
        #skt_url = f"https://cdaweb.gsfc.nasa.gov/pub/software/cdawlib/0SKELTABLES/{dataset_id}_00000000_v01.skt"
        #response = requests.get(skt_url)
        #skeleton_table_text = response.text

        #doi_info = dataset["Doi"]
        #
        #if doi_info:
        #    # Generate the landing page metadata URL
        #    metadata_url = get_doi_landing_page_url(doi_info)
        #else:
        #    # Fallback if no DOI info is provided directly in the dataset object
        #    metadata_url = f"https://nasa.gov{dataset_id}"

        #print('\033[1;32m%s\033[0m -> %s, %s' %(dataset_id, dataset_label,dataset_note))#, display_link(metadata_url, 'metadata')))
        #print("'\033[1m%s\033[0m' -> %s" %(dataset_id, dataset_label))#, display_link(metadata_url, 'metadata')))
        #print("%s" %(dataset_note))#, display_link(metadata_url, 'metadata')))
        #print('')
        #display(Markdown('\033[1;32m%s\033[0m -> %s' %(dataset_id, dataset_label)))#, display_link(metadata_url, 'metadata')))
        #print(f'{index}. {dataset_id}, {dataset_label}')
        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 Wind
Dataset with the link -> Label

OMNI2_H0_MRG1HR -> OMNI Combined, Definitive, Hourly IMF and Plasma Data, and Energetic Proton Fluxes, Time-Shifted to the Nose of the Earth's Bow Shock, plus Solar and Magnetic Indices - J.H. King, N. Papitashvili (ADNET, NASA GSFC)     
OMNI_COHO1HR_MERGED_MAG_PLASMA -> OMNI Combined merged hourly magnetic field, plasma and ephermis data - J.H. King, N. Papatashvilli (AdnetSystems, NASA GSFC)                                                                                
OMNI_HRO2_1MIN -> OMNI Combined, Definitive 1-minute IMF and Definitive Plasma Data Time-Shifted to the Nose of the Earth's Bow Shock, plus Magnetic Indices - J.H. King, N. Papatashvilli (AdnetSystems, NASA GSFC)                          
OMNI_HRO2_5MIN -> OMNI Combined, Definitive 5-minute IMF and Definitive Plasma, and Energetic Proton Fluxes, Time-Shifted to the Nose  of the Earth's Bow Shock, plus Magnetic Indices - J.H. King, N. Papatashvilli (AdnetSystems, NASA GSFC)
OMNI_HRO_1MIN -> OMNI Combined, Definitive, 1-minute IMF and Plasma Data Time-Shifted to the Nose of the Earth's Bow Shock, plus Magnetic Indices - J.H. King, N. Papatashvilli (AdnetSystems, NASA GSFC)                                     
OMNI_HRO_5MIN -> OMNI Combined, Definitive, 5-minute IMF and Plasma, and Energetic Proton Fluxes, Time-Shifted to the Nose  of the Earth's Bow Shock, plus Magnetic Indices - J.H. King, N. Papatashvilli (AdnetSystems, NASA GSFC)           
WI_H0_MFI -> Wind Magnetic Fields Investigation: 3 sec, 1 min, and hourly Definitive Data. - A. Koval (UMBC, NASA/GSFC)                                                                                                                       
WI_H1_SWE -> Wind Solar Wind Experiment, 92-sec Solar Wind Alpha and Proton Anisotropy Analysis - K. Ogilvie (NASA GSFC)                                                                                                                      
WI_H1_SWE_RTN -> Solar wind proton and alpha parameters, including anisotropic temperatures, derived by non-linear fitting of the measurements and with moment techniques. Data in RTN - Bennett Maruca (University of Delaware)              
WI_H2_MFI -> Wind Magnetic Fields Investigation, High-resolution Definitive Data - A. Koval (UMBC, NASA/GSFC)                                                                                                                                 
WI_H3-RTN_MFI -> Wind Magnetic Fields Investigation: 3 sec, 1 min, and hourly Definitive Data (RTN). - A. Koval (UMBC, NASA/GSFC)                                                                                                             
WI_H4-RTN_MFI -> Wind Magnetic Fields Investigation, High-resolution Definitive Data (RTN) - A. Koval (UMBC, NASA/GSFC)                                                                                                                       
WI_K0_GIFWALK -> Links to Wind KP pre-generated survey and other plots - Polar-Wind-Geotail Ground System (NASA GSFC)                                                                                                                         
WI_K0_MFI -> Wind Magnetic Fields Investigation, Key Parameters [PRELIM] - R. Lepping (NASA/GSFC)                                                                                                                                             
WI_K1-RTN_MFI -> Wind Magnetic Fields Investigation, Key Parameters in RTN [PRELIM] - Andriy Koval (UMBC&GSFC)                                                                                                                                

====================================================================================================
In [6]:
### user-specified ###
mydataset=['WI_H0_MFI'] #the order should be consistent with the list above.
In [7]:
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}.'))
        
        #variables = cdas.get_variables(doi)
        #var_names = []
        #for index, variable in enumerate(variables):
        #    name = variable['Name']
        #    var_names.append(name)
        #    print("'%s'" %(variable['Name']),"->", variable['LongDescription'])
        
        print('='*100)
        print(' ')

Selected dataset: WI_H0_MFI (doi:10.48322/av38-wn55)
Available Period: {'Start': '1994-11-13T00:00:30.000Z', 'End': '2026-09-12T23:59:30.000Z'}
HDP Description
CDAWeb Notes
Data Citation: Koval, A., Lepping, R. P., & & Szabo, A. (2021). Wind Magnetic Field Investigation (MFI) Composite Data [Dataset]. Space Physics Data Facility. https://doi.org/10.48322/AV38-WN55. Accessed on 2026-September-25. 
====================================================================================================
In [8]:
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 WI_H0_MFI
Available from: 1994-11-13 00:00:30+00:00 to 2006-11-23 23:59:30+00:00
Available from: 2006-11-28 00:00:30+00:00 to 2009-02-20 23:59:30+00:00
Available from: 2009-02-24 00:00:30+00:00 to 2014-10-27 23:59:30+00:00
Available from: 2014-11-21 00:00:30+00:00 to 2023-07-14 23:59:30+00:00
Available from: 2023-07-26 00:00:30+00:00 to 2026-07-24 23:59:30+00:00
Available from: 2026-07-26 00:00:30+00:00 to 2026-09-12 23:59:30+00:00
 
'BF1' -> Magnetic field magnitude (1 min)
'BF1LOG' -> Magnetic field magnitude (1 min - log scaled)
'BRMSF1' -> RMS magnitude (1 min)
'BGSM' -> Magnetic field vector in GSM cartesian coordinates (1 min)
'BRMSGSM' -> RMS vector in GSM coordinates (1 min)
'BGSE' -> Magnetic field vector in GSE cartesian coordinates (1 min)
'BGSEa' -> Magnetic field vector in GSE angular coordinates (1 min)
'BRMSGSE' -> RMS vector in GSE coordinates (1 min)
'DIST' -> Distance from the center of the earth (Define Re = 6378km) (1 min)
'PGSM' -> Position vector in GSM coordinates (Define Re = 6378km) (1 min)
'PGSE' -> Position vector in GSE coordinates (Define Re = 6378km) (1 min)
'B3F1' -> Magnetic field magnitude (3 sec)
'B3F1LOG' -> Magnetic field magnitude (3 sec - log scaled)
'B3RMSF1' -> RMS magnitude (3 sec)
'B3GSM' -> Magnetic field vector in GSM cartesian coordinates (3 sec)
'B3RMSGSM' -> RMS vector in GSM coordinates (3 sec)
'B3GSE' -> Magnetic field vector in GSE cartesian coordinates (3 sec)
'B3GSEa' -> Magnetic field vector in GSE angular coordinates (3 sec)
'B3RMSGSE' -> RMS vector in GSE coordinates (3 sec)
'B1F1' -> Magnetic field magnitude (1 hour)
'B1F1LOG' -> Magnetic field magnitude (1 hour - log scaled)
'B1RMSF1' -> RMS magnitude (1 hour)
'B1GSM' -> Magnetic field vector in GSM cartesian coordinates (1 hour)
'B1RMSGSM' -> RMS vector in GSM coordinates (1 hour)
'B1GSE' -> Magnetic field vector in GSE cartesian coordinates (1 hour)
'B1GSEa' -> Magnetic field vector in GSE angular coordinates (1 hour)
'B1RMSGSE' -> RMS vector in GSE coordinates (1 hour)
'DIST1' -> Distance from the center of the earth (Define Re = 6378km) (1 hour)
'P1GSM' -> Position vector in GSM coordinates (Define Re = 6378km) (1 hour)
'P1GSE' -> Position vector in GSE coordinates (Define Re = 6378km) (1 hour)
'DISTV' -> Time axis label: Distance from the center of the earth (Define Re = 6378km) (1 min)
'PGSMV' -> Time axis label: Position vector in GSM coordinates (Define Re = 6378km) (1 min)
'PGSEV' -> Time axis label: Position vector in GSE coordinates (Define Re = 6378km) (1 min)
'DIST1V' -> Time axis label: Distance from the center of the earth (Define Re = 6378km) (1 hour)
'P1GSMV' -> Time axis label: Position vector in GSM coordinates (Define Re = 6378km) (1 hour)
'P1GSEV' -> Time axis label: Position vector in GSE coordinates (Define Re = 6378km) (1 hour)
====================================================================================================
In [9]:
# assign parameteprs
#dataset_id = 'AC_H0_SWE'
myvariables  = {}
mystart_time = {}
myend_time   = {}

### user-specified ###
myvariables[mydataset[0]]=['B3GSM']
mystart_time[mydataset[0]] = '1997-05-23T08:00:00Z'
myend_time[mydataset[0]]   = '1997-05-25T00: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}")

### change the name of the first index to "time" (hope the first variable is Epoch-related).
for dataset in mydataset:
    keys=list(data[dataset].indexes.keys())
    data[dataset] = data[dataset].rename({keys[0]: "time"})

WI_H0_MFI
Donwload Success!
In [10]:
print(data['WI_H0_MFI'])

<xarray.Dataset>
Dimensions:    (time: 48000, cartesian: 3)
Coordinates:
  * time       (time) datetime64[ns] 1997-05-23T08:00:01.500000 ... 1997-05-2...
  * cartesian  (cartesian) <U11 'x-component' 'y-component' 'z-component'
    metavar0   (cartesian) <U8 'Bx (GSM)' 'By (GSM)' 'Bz (GSM)'
Data variables:
    B3GSM      (time, cartesian) float32 2.488 -0.4322 -0.1115 ... 1.404 -5.013
Attributes: (12/35)
    Project:                     ['ISTP>International Solar-Terrestrial Physi...
    Discipline:                  ['Space Physics>Heliospheric Science']
    Source_name:                 ['WIND>Wind Interplanetary Plasma Laboratory']
    Software_version:            ['V05.00']
    Data_type:                   ['H0>3 sec, 1 min, and hourly Definitive Data']
    Descriptor:                  ['MFI>Magnetic Fields Investigation']
    ...                          ...
    Mission_group:               ['Wind', '!___Interplanetary Data near 1 AU']
    Web_sites:                   ['http://wind.nasa.gov']
    spase_DatasetResourceID:     ['spase://NASA/NumericalData/Wind/MFI/PT03S']
    DOI:                         ['https://doi.org/10.48322/av38-wn55']
    CDAWEB_PARENTS:              ['wi_h0_mfi_00000000_v01', 'wi_h0_mfi_199705...
    PARENTS:                     ['wi_h0_mfi_00000000_v01', 'wi_h0_mfi_199705...
In [11]:
print(data['WI_H0_MFI'].dims)

Frozen({'time': 48000, 'cartesian': 3})
In [12]:
### fill out NAN values
data['WI_H0_MFI']['B3GSM']=data['WI_H0_MFI']['B3GSM'].interpolate_na(dim='time',method="linear", fill_value="extrapolate")
In [13]:
plt.figure(figsize=(20,5))
plt.plot(data['WI_H0_MFI']['time'],data['WI_H0_MFI']['B3GSM'][:,0])
plt.plot(data['WI_H0_MFI']['time'],data['WI_H0_MFI']['B3GSM'][:,1])
plt.plot(data['WI_H0_MFI']['time'],data['WI_H0_MFI']['B3GSM'][:,2])
plt.ylabel('magnetic fields [nT]')
#plt.xlim('1998-02-04T00:00:00Z','1998-02-10T00:00:00Z')
plt.ylim(-10,10)
plt.show()
plt.show()

Wavelet Analysis to detect flux ropes [Telloni+2012ApJ]¶

The magnetic helicity definied as $H_m=\int{\bf A}\cdot{\bf B}d^3{\bf r}$ can be rewritten in a reduced form, $$H_m^r(k)=\frac{2{\rm Im}[Y^*(k)\cdot Z(k)]}{k}.$$ Here, ${\bf A}$ is the vector potential, and $Y$ and $Z$ are the Fourier coefficients of the components of the magnetic field $B_y$ and $B_z$.

Using the wavelet transform (see the reference above), we can rewrite the reduced magnetic helicity,

$$H_m^r(k,t)=\frac{2{\rm Im}[W_y^*(k,t)\cdot W_z(k,t)]}{k}.$$

Here, $W_y$ and $W_z$ are the wavelet transforms of $B_y$ and $B_z$, respectively.

We also define the normalized magnetic helicity as,

$$\sigma_m^r(k,t)=\frac{kH_m(k,t)}{|W_y(k,t)|^2+|W_z(k,t)|^2}.$$
In [18]:
### continuous wavelet analysis
dt=3/3600  #cadence in units of hours
nt=len(data['WI_H0_MFI']['time'].values)
print(nt)
lt=nt*dt
tt=np.arange(nt)*dt

### choose the wavelet base
wavelet_base='Paul'

### define scale ###
ns=256
scl=np.logspace(-1,np.log2(lt),ns,base=2)
#print(a)

### define frequency range ###
wk=2*np.pi*sf.fftfreq(nt,dt)


Wy=np.zeros((ns,nt),dtype=complex)
Wz=np.zeros((ns,nt),dtype=complex)
for i in range(ns):
    sw=scl[i]*wk
    if wavelet_base=='Mexican hat':
        m=2
        fftphi0=-1j*1j/np.sqrt(3/4*np.sqrt(np.pi))*(sw)**2*np.exp(-0.5*sw*sw) #mexican hat
    elif wavelet_base=='Morlet':
        w0=6.
        sw2=(sw-w0)*(sw-w0)
        fftphi0=np.pi**(-0.25)*np.heaviside(wk,0.)*np.exp(-0.5*sw*sw)
    elif wavelet_base=='Paul':
        m=4
        fftphi0=2**m/np.sqrt(m*mt.factorial((2*m-1)))*np.heaviside(wk,0)*sw**m*np.exp(-abs(sw)) ### I put abs(swk) to avoid the overflow.
   
    fftphi =np.sqrt(2*np.pi*scl[i]/dt)*fftphi0
    fftby  =sf.fft(data['WI_H0_MFI']['B3GSM'][:,1].values)
    fftbz  =sf.fft(data['WI_H0_MFI']['B3GSM'][:,2].values)
    Wy[i,:]=sf.ifft(fftby*np.conj(fftphi))
    Wz[i,:]=sf.ifft(fftbz*np.conj(fftphi))
    
    if(i==10):
        plt.plot(scl[i]*wk/(2*np.pi),fftphi0)
        plt.xlabel(r'$s\omega$')
        plt.ylabel(r'$\hat{\Psi}_0(s\omega)$')
        plt.xlim(-2,2)
        plt.show()
    
### calculate magnetic helicity
Hm=2*np.imag(np.conj(Wy)*Wz).T*scl
### calculate normalized magnetic helicity
sigm=Hm/scl/(abs(Wy.T)**2+abs(Wz.T)**2)

### scaling ###
zz=Hm*np.abs(sigm)*scl**(-8./3.)

### plotting ###
fig,ax=plt.subplots(figsize=(20,10))
X,Y=np.meshgrid(tt,scl)
#im=ax.imshow(zz.T,extent=[0,nt*dt,a[na-1],a[0]],vmin=-2e3,vmax=2e3,aspect='auto',cmap='jet')
im=ax.contourf(X,Y,zz.T,cmap='jet',levels=100)#,vmin=-4e3,vmax=4e3)
#ax.contour( X,Y,zz.T,colors='k',levels=10)
#fig.colorbar(im,ax=ax)
ax.set_yscale('log',base=2)
ax.set_xlim(0,lt)
ax.set_ylim(scl[ns-1],scl[0])
ax.set_xlabel('time [h]')
ax.set_ylabel('scale [h]')
ax.set_title('Magnetic Helicity')
plt.show()

48000