MAVEN Data Access Example¶

Tamar Ervin

September 21, 2023

• Downloading MAVEN data with PySPEDAS
• PySPEDAS is unable to read the CDF files, so it will break with the CDF reader but still download the files!

In [ ]:

import glob
import pyspedas
import cdflib
import sunpy
import astrospice
import sunpy 
import sys, os
import datetime
import numpy as np
import pandas as pd
import astropy.units as u
import matplotlib.pyplot as plt

import glob import pyspedas import cdflib import sunpy import astrospice import sunpy import sys, os import datetime import numpy as np import pandas as pd import astropy.units as u import matplotlib.pyplot as plt

Download and Plot the Data¶

In [ ]:

### ------- TIME PERIOD OF INTEREST ------- ###
time_range = ['2021-08-10', '2021-08-11']

### ------- TIME PERIOD OF INTEREST ------- ### time_range = ['2021-08-10', '2021-08-11']

We start by using pyspedas to pull the MAVEN MAG data. There are some weird things with these files so we have to use some other functions to actually read in the files.

In [ ]:

### ------- MAVEN: MAG RTN DATA ------- ###
pyspedas.maven.mag(trange=time_range, datatype='ss1s', downloadonly=True)

### ------- MAVEN: MAG RTN DATA ------- ### pyspedas.maven.mag(trange=time_range, datatype='ss1s', downloadonly=True)

These are the functions to read in the MAVEN MAG files.

Courtesy of Rebecca Jolitz.

In [ ]:

import datetime as dt
def read_sts(filename, method='loadtxt'):

    if "mvn_mag_ql" in filename:
        skip_header = 86
    elif "mvn_mag_l2" in filename:
        skip_header = 150

    if method == "loadtxt":
        b = np.loadtxt(filename, skiprows=skip_header, dtype='f8')
    elif method == "genfromtxt":
        b = np.genfromtxt(
            filename,
            delimiter=(6, 4, 3, 3, 3, 4, 14, 11, 10, 10, 4),
            skip_header=skip_header, dtype='f8')

    return b

def doy_to_utc(year, doy, hour, minute, sec, msec):

    # Sometimes these fields are
    # improperly read as integers,
    # which will cause strange additive behaviors.
    # Convert them to floats to avoid this.
    year = year.astype('float')
    doy = doy.astype('float')
    hour = hour.astype('float')
    minute = minute.astype('float')
    sec = sec.astype('float')
    msec = msec.astype('float')

    # Get the total number of seconds that have transpired
    # since the year started.
    total_sec = (doy - 1)*24*60*60 + hour*60*60 +\
        minute*60 + sec + (1e-3*msec)
    time_utc =\
        [dt.datetime(int(y), 1, 1) + dt.timedelta(seconds=s) for (y, s) in
         zip(year, total_sec)]

    return time_utc


def parse_sts(b, columns=("epoch", "Bx", "By", "Bz")):

    # L1 has 11 columns:
    # Year | DOY | Hour | Min | Sec | Msec | Decimal day |
    # (Outboard magnetic field columns)
    # OB_B_X (nT)| OB_B_Y (nT)| OB_B_Z (nT)| OB_B_RANGE (nT)|

    # L2 has 18 columns:
    # Year | DOY | Hour | Min | Sec | Msec | Decimal day |
    # (Outboard magnetic field columns)
    # OB_B_X (nT)| OB_B_Y (nT)| OB_B_Z (nT)| OB_B_RANGE (nT)|
    # (S/c position columns)
    # POSN_X (km) | POSN_Y (km) | POSN_Z (km)
    # (Outboard dynamic corrections)
    # OB (nT)| OB_BD_Y (nT)| OB_BD_Z (nT)| OB_BD_RANGE (nT)|

    time_utc = doy_to_utc(b[:, 0], b[:, 1], b[:, 2], b[:, 3], b[:, 4], b[:, 5])

    bx = b[:, 7]
    by = b[:, 8]
    bz = b[:, 9]

    return time_utc, bx, by, bz

import datetime as dt def read_sts(filename, method='loadtxt'): if "mvn_mag_ql" in filename: skip_header = 86 elif "mvn_mag_l2" in filename: skip_header = 150 if method == "loadtxt": b = np.loadtxt(filename, skiprows=skip_header, dtype='f8') elif method == "genfromtxt": b = np.genfromtxt( filename, delimiter=(6, 4, 3, 3, 3, 4, 14, 11, 10, 10, 4), skip_header=skip_header, dtype='f8') return b def doy_to_utc(year, doy, hour, minute, sec, msec): # Sometimes these fields are # improperly read as integers, # which will cause strange additive behaviors. # Convert them to floats to avoid this. year = year.astype('float') doy = doy.astype('float') hour = hour.astype('float') minute = minute.astype('float') sec = sec.astype('float') msec = msec.astype('float') # Get the total number of seconds that have transpired # since the year started. total_sec = (doy - 1)*24*60*60 + hour*60*60 +\ minute*60 + sec + (1e-3*msec) time_utc =\ [dt.datetime(int(y), 1, 1) + dt.timedelta(seconds=s) for (y, s) in zip(year, total_sec)] return time_utc def parse_sts(b, columns=("epoch", "Bx", "By", "Bz")): # L1 has 11 columns: # Year | DOY | Hour | Min | Sec | Msec | Decimal day | # (Outboard magnetic field columns) # OB_B_X (nT)| OB_B_Y (nT)| OB_B_Z (nT)| OB_B_RANGE (nT)| # L2 has 18 columns: # Year | DOY | Hour | Min | Sec | Msec | Decimal day | # (Outboard magnetic field columns) # OB_B_X (nT)| OB_B_Y (nT)| OB_B_Z (nT)| OB_B_RANGE (nT)| # (S/c position columns) # POSN_X (km) | POSN_Y (km) | POSN_Z (km) # (Outboard dynamic corrections) # OB (nT)| OB_BD_Y (nT)| OB_BD_Z (nT)| OB_BD_RANGE (nT)| time_utc = doy_to_utc(b[:, 0], b[:, 1], b[:, 2], b[:, 3], b[:, 4], b[:, 5]) bx = b[:, 7] by = b[:, 8] bz = b[:, 9] return time_utc, bx, by, bz

In [ ]:

### READ IN TEXT FILES
mag_file_path = glob.glob(os.path.realpath(os.path.join('maven_data/maven/data/sci', 'mag/l2', '2021', '08', '*.sts'))) ## path to the .sts magnetic field files
time_utc, bx, by, bz = parse_sts(read_sts(mag_file_path[1])) ## read in a specific file

### CALCULATE THE MAGNETIC FIELD MAGNITUDE
Btotal = np.sqrt(bx**2 + by**2 + bz**2)

### READ IN TEXT FILES mag_file_path = glob.glob(os.path.realpath(os.path.join('maven_data/maven/data/sci', 'mag/l2', '2021', '08', '*.sts'))) ## path to the .sts magnetic field files time_utc, bx, by, bz = parse_sts(read_sts(mag_file_path[1])) ## read in a specific file ### CALCULATE THE MAGNETIC FIELD MAGNITUDE Btotal = np.sqrt(bx**2 + by**2 + bz**2)

We now plot the x, y, z and total magnetic field!

In [ ]:

### PLOT THE MAG DATA
fig = plt.figure(figsize=[15, 4])
plt.plot(time_utc, bx, color='tab:blue', label=r'$\rm B_x$')
plt.plot(time_utc, by, color='tab:red', label=r'$\rm B_y$')
plt.plot(time_utc, bz, color='tab:green', label=r'$\rm B_z$')
plt.plot(time_utc, Btotal, color='black', label=r'$\rm |B|$')
plt.ylabel(r'$\rm B \; [nT]$')
plt.legend(loc='upper right', fontsize=18)

### PLOT THE MAG DATA fig = plt.figure(figsize=[15, 4]) plt.plot(time_utc, bx, color='tab:blue', label=r'$\rm B_x$') plt.plot(time_utc, by, color='tab:red', label=r'$\rm B_y$') plt.plot(time_utc, bz, color='tab:green', label=r'$\rm B_z$') plt.plot(time_utc, Btotal, color='black', label=r'$\rm |B|$') plt.ylabel(r'$\rm B \; [nT]$') plt.legend(loc='upper right', fontsize=18)

Now we will look at some particle data!

In [ ]:

### ------- SWIA: ION MOMENTS ------- ###
### download ion particle data
pyspedas.maven.swia(trange=time_range, datatype='onboardsvymom', downloadonly=True)

### ------- SWIA: ION MOMENTS ------- ### ### download ion particle data pyspedas.maven.swia(trange=time_range, datatype='onboardsvymom', downloadonly=True)

We will use a special CDF reader to read in the files.

In [ ]:

### READ IN THE ION (SWIA) FILES
files = glob.glob(os.path.join(os.path.realpath(os.path.join('maven_data', 'maven/data/sci/swi/l2/2021/08')), "*.cdf"), recursive=True)
data_cdf = cdflib.CDF(files[0]) ## read in the first ion file

### PRINT OUT FILE INFO
print(data_cdf.cdf_info())

### READ IN VARIABLE INFO
N = data_cdf.varget("density")
v = data_cdf.varget("velocity")
T = data_cdf.varget("temperature")
p = data_cdf.varget("pressure")

### PLOT THE DATA
fig, axs = plt.subplots(4, figsize=[14, 10])
dt = [N, v, T, p]
ylabels = [r'$\rm N_p$', r'$\rm v_p$', r'$\rm T_p$', r'$\rm p_p$']
for i, ax in enumerate(axs):
    ax.plot(dt[i], linewidth=1)
    ax.set_ylabel(ylabels[i], fontsize=16)

### READ IN THE ION (SWIA) FILES files = glob.glob(os.path.join(os.path.realpath(os.path.join('maven_data', 'maven/data/sci/swi/l2/2021/08')), "*.cdf"), recursive=True) data_cdf = cdflib.CDF(files[0]) ## read in the first ion file ### PRINT OUT FILE INFO print(data_cdf.cdf_info()) ### READ IN VARIABLE INFO N = data_cdf.varget("density") v = data_cdf.varget("velocity") T = data_cdf.varget("temperature") p = data_cdf.varget("pressure") ### PLOT THE DATA fig, axs = plt.subplots(4, figsize=[14, 10]) dt = [N, v, T, p] ylabels = [r'$\rm N_p$', r'$\rm v_p$', r'$\rm T_p$', r'$\rm p_p$'] for i, ax in enumerate(axs): ax.plot(dt[i], linewidth=1) ax.set_ylabel(ylabels[i], fontsize=16)

Read in the EUV Data

In [ ]:

files = glob.glob(os.path.join(os.path.realpath(os.path.join('maven_data', 'maven/data/sci/euv/l2/2016/01')), "*"), recursive=True)
data_cdf = cdflib.CDF(files[0]) ## read in the first EUV file

### PRINT OUT FILE INFO
print(data_cdf.cdf_info())

### Get the spectral irradiance (W/m^2) as a timeseries
dt = data_cdf.varget("data")
time = data_cdf.varget("time_unix")

### PLOT THE SPECTRAL IRRANDIANCE
fig = plt.figure(figsize=[15, 4])
plt.plot(time, dt)
plt.ylabel('Spectral Irradiance')

files = glob.glob(os.path.join(os.path.realpath(os.path.join('maven_data', 'maven/data/sci/euv/l2/2016/01')), "*"), recursive=True) data_cdf = cdflib.CDF(files[0]) ## read in the first EUV file ### PRINT OUT FILE INFO print(data_cdf.cdf_info()) ### Get the spectral irradiance (W/m^2) as a timeseries dt = data_cdf.varget("data") time = data_cdf.varget("time_unix") ### PLOT THE SPECTRAL IRRANDIANCE fig = plt.figure(figsize=[15, 4]) plt.plot(time, dt) plt.ylabel('Spectral Irradiance')