PSP_SPAN-I_data_tutorial_pyspedas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting started with plotting PSP data using pyspedas.\n",
    "\n",
    "Greetings! Below is a guide to plotting plasma moments from the PSP SPAN-I instrument using pySPEDAS, the python implementation of the SPEDAS library originally developed in IDL.\n",
    "\n",
    "We note that the SPAN-I instrument has many caveats and the data remains preliminary due to ongoing calibration work done by the instrument team. \n",
    "\n",
    "For optimal scientific interpretation of noteworthy event, the reader is strongly encouraged to contact a member of the instrument team:<br>\n",
    "Roberto Livi rlivi@berkeley.edu <br>\n",
    "Ali Rahmati rahmati@berkeley.edu <br>\n",
    "Michael McManus mdmcmanus@berkeley.edu <br>\n",
    "Davin Larson davin@berkeley.edu <br>\n",
    "\n",
    "This tutorial was brought to you by: <br>\n",
    "Jaye Verniero jaye.l.verniero@nasa.gov <br>\n",
    "Kristoff Paulson kristoff.paulson@cfa.harvard.edu <br>\n",
    "who may serve as points of contact for this tutorial. <br>\n",
    "\n",
    "The paper describing the SPAN-I instrument in more detail can be found here: <br>https://doi.org/10.3847/1538-4357/ac93f5\n",
    "\n",
    "First, we install pyspedas using instructions from this link: https://github.com/spedas/pyspedas\n",
    "\n",
    "Follow all steps under the Installation Section.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#one could alternatively try:\n",
    "%pip install --upgrade pyspedas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyspedas\n",
    "from pytplot import tplot, store_data, get_data ,tlimit,xlim,ylim,tplot_options,options,split_vec\n",
    "from pyspedas import time_string, time_double\n",
    "\n",
    "#import math functions \n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify time range in the form ['yyyy-mm-dd/hh:mm:ss','yyyy-mm-dd/hh:mm:ss']\n",
    "trange=['2020-01-29','2020-01-30']\n",
    "\n",
    "#specify data type to plot\n",
    "datatype='spi_sf00_l3_mom' #protons\n",
    "spi_vars = pyspedas.psp.spi(trange=trange, datatype=datatype, level='l3', time_clip=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print variable names\n",
    "print(spi_vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is a description of each variable:<br>\n",
    "QUALITY FLAG - Data quality indicator <br>\n",
    "DENS - Density (1/cm^3) <br>\n",
    "VEL_INST - Velocity in SPAN-I instrument coordinates (km/s) <br>\n",
    "VEL_SC - Spacecraft Velocity (km/s)<br>\n",
    "VEL_RTN_SUN - Ion Velocity in RTN coordinates (km/s)<br>\n",
    "TEMP - Ion Temperature (eV)<br>\n",
    "T_TENSOR_INST - Temperature Tensor SPAN-I instrument coordinates <br>\n",
    "EFLUX_VS_ENERGY - Proton Differential Energy Flux vs eV <br>\n",
    "EFLUX_VS_THETA - Proton Differential Energy Flux vs Theta (deflectors) <br>\n",
    "EFLUX_VS_PHI - Proton Differential Energy Flux vs Phi (anodes) <br>\n",
    "SUN-DIST - Distance from Sun (km)<br>\n",
    "VENUS-DIST - Distance from Venus (km)<br>\n",
    "SC_VEL_RTN- Spacecraft Velocity in RTN coordinates (km/s)<br>\n",
    "QUAT_SC_TO_RTN - Quaternion rotation of spacecraft in RTN coordinates<br>\n",
    "MAGF_SC - Magnetic field downsampled to SPAN_I time resolution in spacecraft coordinates (nT)<br>\n",
    "MAGF_INST - Magnetic field downsampled to SPAN_I time resolution in SPAN-I instrument coordinates (nT)<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Indexing in Python starts at 0\n",
    "\n",
    "#first variable name in list\n",
    "print(spi_vars[0])\n",
    "\n",
    "#fist 4 variable names in list\n",
    "print(spi_vars[0:3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that pySPEDAS software is still under development. The variables in the cdf file may not have a prefix, while pySPEDAS may assign one. This is currently an issue that will be updated. For flexibilty, we define a prefix for variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prefix = 'psp_spi_'\n",
    "#prefix = ''\n",
    "\n",
    "#plot variable vs time\n",
    "tplot(prefix+'VEL_RTN_SUN')\n",
    "#set legend names \n",
    "options(prefix+'VEL_RTN_SUN', 'legend_names', ['Vr', 'Vt', 'Vn'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#example plotting multple variables\n",
    "tplot([prefix+'DENS',prefix+'VEL_RTN_SUN'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify shorter time range in format 'yyyy-mm-dd hh:mm:ss'\n",
    "start_day=trange[0]\n",
    "trange1=[start_day+' 18:00:00',start_day+ ' 20:00:00']\n",
    "#set limits \n",
    "xlim(trange1[0],trange1[1])\n",
    "\n",
    "#replot\n",
    "tplot([prefix+'DENS',prefix+'VEL_RTN_SUN'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify limits on y-axis\n",
    "ylim(prefix+'DENS',4E2,2E3)\n",
    "ylim(prefix+'VEL_RTN_SUN',-100,400)\n",
    "\n",
    "#change line thickness\n",
    "options([prefix+'DENS',prefix+'VEL_RTN_SUN'],'thick',2)\n",
    "\n",
    "#change plot size in x and y direction\n",
    "tplot([prefix+'DENS',prefix+'VEL_RTN_SUN'], xsize=15, ysize=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Investigation 1\n",
    "The dates of the fourth orbit around the Sun, also known as encounter 4 (E4) were roughly from 2020-01-23 to 2020-02-04. Plot the distance from the sun, velocity in instrument coordinates, and the magnetic field in instrument coordinates for the whole encounter. What trends do you notice? Are there any similarities or differences between the inbound (coming toward the Sun) vs outbound (moving away from the Sun) of the orbit? Play with different time ranges. Add Density and Temperature. Do you notice any patterns between different variables?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make cells below by choosing insert --> cell below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to make new tplot variables\n",
    "\n",
    "This section demonstrates how to use the functions 'get_data' and 'store_data' to make new tplot variables.\n",
    "The 'get_data' function extracts the data from the tplot Variables stored in memory. From there, one can manipilate variables and make new ones. Those new variables are then stored as tplot variables with the 'store_data'command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#assign data structure to tplot variable\n",
    "dens_data=get_data(prefix+'DENS',metadata=True)\n",
    "\n",
    "#get metadata\n",
    "dens_data['CDF']['VATT']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#turn off metadata, (default is set false)\n",
    "dens_data=get_data(prefix+'DENS')\n",
    "dens_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output states this variable is made of two arrays called \"times\" and \"y.\" Here is how to access those arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times=dens_data.times\n",
    "dens_array = dens_data.y\n",
    "\n",
    "#check the shape of each array\n",
    "print(np.shape(times))\n",
    "print(np.shape(dens_array))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Multiply two variables together, such as density and temperature\n",
    "\n",
    "#define new data structure name for temperature\n",
    "temp_data = get_data(prefix+'TEMP')\n",
    "\n",
    "#in tplot variables, the first dimension is always \"times\", so we do not need to call it again\n",
    "temp_array = temp_data.y\n",
    "#note how temperature and density have same number of elements\n",
    "print(np.shape(temp_array))\n",
    "\n",
    "#We can multply two arrays with same number of elements together\n",
    "NT_p = dens_array*temp_array\n",
    "print(np.shape(NT_p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define Boltzmann's constant (eV/K)\n",
    "k = 8.617333262E-5\n",
    "kNT_p = k*NT_p\n",
    "\n",
    "#make new tplot variable called 'nT_p' with kNT_p on y axis, and time on x-axis\n",
    "store_data('nKT_p',data={'x':times,'y':kNT_p})\n",
    "#plot\n",
    "tplot('nKT_p')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The magitude of the velocity $|\\mathbf{v}|$ is often approximated as fair proxy for the solar wind speed. Here is an example of how to compute $|\\mathbf{v}|$ and overplot it with $|\\mathbf{v}|=(v_R,v_T,v_N)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vel_rtn=get_data(prefix+'VEL_RTN_SUN')\n",
    "print(vel_rtn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After printing variable information, we see that there is the \"times\" array, but this time the \"y\" array has 3 columns. This means that the \"y\" array is a matrix that can be accessed using indexing. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#access time indexed by 100 from of the radial (first) component of v\n",
    "vr_t100=vel_rtn.y[100,0]\n",
    "print(vr_t100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define times\n",
    "vtimes=vel_rtn.times\n",
    "\n",
    "print(vtimes[100])\n",
    "\n",
    "#convert from unix time unsing time_string function\n",
    "print(time_string(vtimes[100]))\n",
    "\n",
    "#see times entries 90-100\n",
    "print(time_string(vtimes[89:99]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the ':' symbol accesses everything in first column (time)\n",
    "#access vr, vt, vn\n",
    "vr = vel_rtn.y[:,0]\n",
    "vt = vel_rtn.y[:,1]\n",
    "vn = vel_rtn.y[:,2]\n",
    "\n",
    "#check all arrays are same shape\n",
    "print(vr.shape,vt.shape,vn.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#define magnitude of velocity vector\n",
    "vmag=np.sqrt(vr**2 + vt**2 + vn**2) #note '**' symbol is for exponents\n",
    "print(np.shape(vmag))\n",
    "\n",
    "\n",
    "#define new tplot variable\n",
    "store_data('vmag',data={'x':vtimes,'y':vmag})\n",
    "\n",
    "#plot\n",
    "tplot('vmag')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#now combine new |v| variable with original velocity vector\n",
    "store_data(\"vel_rtn_mag\",data=[prefix+'VEL_RTN_SUN','vmag'])\n",
    "\n",
    "tplot('vel_rtn_mag')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#note how one can also plot individual components of a vector using the 'split_vec' function\n",
    "split_vec(prefix+'VEL_RTN_SUN')\n",
    "tplot(prefix+'VEL_RTN_SUN_x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Investigation 2\n",
    "\n",
    "How much of the velocity magnitude $|\\mathbf{v}|$ is in the radial (R) direction vs the other T and N components? <br>\n",
    "To estimate, <br>\n",
    "1. Make a new tplot variable where the data displayed is instead $\\left(\\frac{v_R}{|\\mathbf{v}|},\\frac{v_T}{|\\mathbf{v}|},\\frac{v_N}{|\\mathbf{v}|}\\right)$ <br>\n",
    "What do you see? Is a trend noticable from the plots? Play with the time ranges again. <br>\n",
    "Think about how you would calculate what fraction of the velocity was in the radial direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to compare the velocity to the magnetic field. Let's use the instrument coordinate coordinate system for now. <br>\n",
    "2. Make a new tplot variable from 'VEL_INST' that displays $\\left(\\frac{v_x}{|\\mathbf{v}|},\\frac{v_y}{|\\mathbf{v}|},\\frac{v_z}{|\\mathbf{v}|}\\right)$ <br>\n",
    "3. Make a new tplot variable from 'MAGF_INST' that displays $\\left(\\frac{b_x}{|\\mathbf{b}|},\\frac{b_y}{|\\mathbf{b}|},\\frac{b_z}{|\\mathbf{b}|}\\right)$ <br>\n",
    "\n",
    "What do you notice? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (pySPEDAS)",
   "language": "python",
   "name": "pyspedas"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}