Skip to content
Snippets Groups Projects
Commit f4545dcd authored by simon-ast's avatar simon-ast
Browse files

Organizing

parent 1f944af9
No related branches found
No related tags found
No related merge requests found
Showing
with 0 additions and 276 deletions
File deleted
File deleted
File deleted
File deleted
File deleted
import os
import typing as tp
import numpy as np
# CDF library (is needed to interface with the measurement data files)
# see https://cdf.gsfc.nasa.gov/
os.environ["CDF_LIB"] = "/usr/local/cdf/lib"
def cdf_slice(cdf_file, key: str):
"""
DOC
:param cdf_file:
:param key:
:return:
"""
return cdf_file[key][...]
def data_generation(cdf_file) -> tp.Dict:
"""
DOC
:param cdf_file:
:return:
"""
data = {
"dqf": cdf_slice(cdf_file, key="general_flag"),
"epoch": cdf_slice(cdf_file, key="Epoch"),
"pos": cdf_slice(cdf_file, key="sc_pos_HCI"),
"vr": cdf_slice(cdf_file, key="vp_moment_RTN")[:, 0],
"dvrhi": cdf_slice(cdf_file, key="vp_moment_RTN_deltahigh")[:, 0],
"dvrlo": cdf_slice(cdf_file, key="vp_moment_RTN_deltalow")[:, 0],
"wp": cdf_slice(cdf_file, key="wp_moment"),
"dwphi": cdf_slice(cdf_file, key="wp_moment_deltahigh"),
"dwplo": cdf_slice(cdf_file, key="wp_moment_deltalow")
}
return data
def array_reduction(data_array: np.ndarray,
index_list: np.ndarray) -> np.ndarray:
"""
Reduction of an array by a list of indices.
:param data_array: NDARRAY,
Array that is to be reduced
:param index_list: NDARRAY,
List of indices that should be deleted
:return: NDARRAY,
The reduced array
"""
# Choose axis=0 to delete rows and not flatten the array
reduced_array = np.delete(data_array, index_list, axis=0)
return reduced_array
def full_reduction(data_dict: tp.Dict,
bad_indices: np.array) -> tp.Dict:
"""
DOC
:param data_dict:
:param bad_indices:
:return:
"""
for key in data_dict.keys():
data_dict[key] = array_reduction(data_dict[key], bad_indices)
return data_dict
import numpy as np
import typing as tp
def general_flag(general_flag_array: np.ndarray) -> np.ndarray:
"""
This function reads out index locations where data general flag is
set to 0 (see documentation for PSP data). This requires correct
handling of CDF array keys
:param general_flag_array: NDARRAY,
cdf_data["general_flag"][...] - Only good when set to 0. The
input here must specifically be this array!
:return: NDARRAY,
An array of the "faulty" indices. These indices correspond to
the index m in the m x n GENERAL CDF array, not a measurement
subarray!
"""
# Notable indices are where values are set to non-zero
# From the data user guide: 0 means no condition present
index_list = np.where(general_flag_array != 0)[0]
return index_list
def meas_failed(meas_array: np.ndarray) -> np.ndarray:
"""
Takes a 1D array and returns indices of failed measurements (as
indicated by a value of -1e30, SWEAP documentation)
:param meas_array: NDARRAY,
Measurement array
:return: NDARRAY,
Index array
"""
index_array = np.where(meas_array <= -0.5e30)[0]
return index_array
def full_meas_eval(data_dict: tp.Dict):
"""DOC"""
col_ind = []
for key in data_dict.keys():
# Skip unnecessary keys (maybe make this more dynamic?)
if key in ["pos", "dqf", "epoch"]:
continue
col_ind.append(meas_failed(data_dict[key]))
# Make an array of unique indices
col_ind = np.unique(np.concatenate(col_ind))
return col_ind
import numpy as np
import typing as tp
# The warnings (for me) here don't seem to matter
from astropy.constants import k_B, m_p
def pos_cart_to_sph(position_matrix: np.ndarray) -> tp.Tuple:
"""
Transforms cartesian coordinates to spherical coordinates.
:param position_matrix: NDARRAY,
(n x 3) matrix of (x, y, z) coordinates
:return: TUPLE,
r = heliocentric distance in KM
theta = inclination in RADIANS
phi = azimuth in RADIANS
"""
# Sanity check: pos_mat MUST BE an (n x 3) matrix
assert position_matrix.shape[1] == 3, \
"PLEASE CHECK POSITION MATRIX FOR COORDINATE CONVERSION!"
x = position_matrix[:, 0]
y = position_matrix[:, 1]
z = position_matrix[:, 2]
r = np.sqrt(x * x + y * y + z * z)
theta = np.arccos(z / r)
phi = np.arctan2(y, x)
return r, theta, phi
def wp_to_temp(thermal_speed: np.ndarray) -> np.ndarray:
"""
Calculate temperature from thermal speed of protons.
:param thermal_speed: NDARRAY,
Thermal speed (wp) of protons
:return: NDARRAY,
Temperature from wp = sqrt(2k_BT / m)
"""
wp_si = thermal_speed * 1e3
return wp_si ** 2 * m_p.value / (2 * k_B.value)
import matplotlib as mpl
import matplotlib.pyplot as plt
def rc_setup():
"""Generalized plot attributes"""
mpl.rcParams["xtick.direction"] = "in"
mpl.rcParams["ytick.direction"] = "in"
mpl.rcParams["xtick.minor.visible"] = "True"
mpl.rcParams["ytick.minor.visible"] = "True"
def plot_r_vr(r, vr, v_ulim, v_lolim, save_ind: str = "no"):
"""
Plots radial velocity wrt heliocentric distance
:param r: NDARRAY,
Distance (must be in km)
:param vr: NDARRAY,
Radial velocity (must be in km/s)
:param v_ulim, v_lolim: NDARRAY,
Upper (lower) uncertainty of vr
:param save_ind: STR,
"yes" if plot should be saved (default: no)
"""
fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(r / 6.957e5, vr)
plt.fill_between(r / 6.957e5, vr + v_ulim, vr - v_lolim, color="grey",
alpha=0.25)
ax.set(
xlabel="r [R$_\\odot$]", ylabel="v$_r$ [km s$^{-1}$]",
ylim=(0, 800)
)
ax.grid(True, alpha=0.5)
if save_ind.lower() == "yes":
plt.savefig("plot_vr.png")
def plot_r_temp(r, T, T_ulim, T_lolim, save_ind: str = "no"):
"""
Plots logarithmic temperature wrt heliocentric distance
:param r: NDARRAY,
Heliocentric distance (mst be in km)
:param T: NDARRAY,
logarithmic Temperature (Kelvin)
:param T_ulim, T_lolim: NDARRAY,
Upper (lower) uncertainty of logT
:param save_ind: STR,
"yes" if plot should be saved (default: no)
"""
fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(r / 6.957e5, T)
plt.fill_between(r / 6.957e5, T_ulim, T - T_lolim,
color="grey", alpha=0.25)
ax.set(
xlabel="r [R$_\\odot$]", ylabel="T [K]", yscale="log"
)
ax.grid(True, alpha=0.5)
if save_ind.lower() == "yes":
plt.savefig("plot_logT.png")
File deleted
File deleted
File deleted
import os
import numpy as np
import matplotlib.pyplot as plt
from spacepy import pycdf
from PSPoperations import data_quality as dq
from PSPoperations import data_transformation as dt
# CDF library (is needed to interface with the measurement data files)
# see https://cdf.gsfc.nasa.gov/
os.environ["CDF_LIB"] = "/usr/local/cdf/lib"
# open CDF file and generate faulty index array
cdf_data = pycdf.CDF("psp_swp_spc_l3i_20211117_v02.cdf")
dqf = cdf_data["general_flag"][...]
bad_ind = dq.general_flag(dqf)
# generate reduced epoch array
time = dq.array_reduction(cdf_data["Epoch"][...], bad_ind)
# generate reduced position array and spherical coordinates
pos = dq.array_reduction(cdf_data["sc_pos_HCI"][...], bad_ind)
r, theta, phi = dt.pos_cart_to_sph(pos)
# generate (prelim.) velocity array (is [0] in RTN right?)
vr = dq.array_reduction(cdf_data["vp_moment_RTN"][...][:, 0], bad_ind)
del_vr_lo = dq.array_reduction(cdf_data["vp_moment_RTN_deltalow"][...][:, 0],
bad_ind)
del_vr_hi = dq.array_reduction(cdf_data["vp_moment_RTN_deltahigh"][...][:, 0],
bad_ind)
plt.figure()
plt.plot(r / 1.49e8, vr)
plt.fill_between(r / 1.49e8, vr + del_vr_hi, vr - del_vr_lo, color="grey",
alpha=0.15)
plt.ylim(0, 1000)
plt.figure()
plt.scatter(time, r / 1.49e8)
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment