spedas / pyspedas

Python-based Space Physics Environment Data Analysis Software
https://pyspedas.readthedocs.io/
MIT License
145 stars 58 forks source link

More pseudovariable issues #661

Closed jameswilburlewis closed 5 months ago

jameswilburlewis commented 8 months ago

I attach a plot and an associated crib sheet. The crib sheet ends where it is noted in a comment (ENDS HERE). Just above it the plot is produced by a call to tplot. You can view or run in command line.

The issues I see is that the options: ‘ytitle’ and ‘ysubtitle’ (as marked) as well as ‘ylog’ (not shown) do not work for the pseudovariable that contains spectra: 'des_eflux_omni_fast_wscpot'. I may have reported this before. (However these work for scalar pseudovariables as 'dxs_Tparper_fast' above it.)

Another issue that I see (not previously reported) is that when I define a plot title using tplot_options(‘title’,’MMS3_Overview’), that title appears 4 times in the plot, once correctly at the top and 3 other times incorrectly above each pseudovariable.

These are of less importance, though. The most important thing to fix remains the ability to change colors in scalar pseudovariables, a critical means to compare scientific quantities produced by different methods.

mms_ovr_BE_FPImoms

jameswilburlewis commented 8 months ago
# This PySPEDAS crib sheet (.pyspd) shows how to introduce MMS L2 FPI data
#   and combine it with data from the EFI and FGM instruments. Then produce your
#   own moments and use the spacecraft potential to affect the moments production,
#   as well as to produce energy, theta, phi, pitch-angle and gyrophase spectra.
#
# Run this in your python console. Event chosen is published in: Torbert et al., 2018
# “ Electron‐scale dynamics of the diffusion region during symmetric magnetic reconnection in space,”
# Science 362, 1391–1395 (2018). https://doi.org/10.1126/science.aat2998
#
# Import basic pyspedas, os, and numpi stuff
pip install pyspedas --upgrade
import pyspedas 
from pyspedas import * # import all top level routines from pyspedas
from pyspedas.mms.cotrans import * # mms_qcotrans, others... not enough to import all pyspedas: need this too!
import pytplot
from pytplot import * # and from pytplot including tplot, get_data, options, tplot_options etc.
import numpy as np
from numpy import * # and from numpy including pi, mean, tan, and sqrt
from copy import * # imports deepcopy
import os
from os import * # and from os including chdir, getcwd
chdir(r"C:\My Documents\ucla\teaching\EPSS_265_2023_Fall\Lecture17_Particle_Counts_to_DFs")
print("Now in Directory:", getcwd()) # print cwd
#
# first decide on the spacecraft and time
mmsnum='3'
mysc='mms'+mmsnum
myscu=mysc+'_'
timerange_ascii = ['2017/07/11 22:10:00', '2017/07/11 22:50:00']
timerange_double=time_double(timerange_ascii)
time2zoom = ['2017/07/11 22:31:00', '2017/07/11 22:36:00']
time2zoom_more = ['2017/07/11 22:33:50', '2017/07/11 22:34:20']
# 
## load position and attitude info
#mec_tvars = pyspedas.mms.state(trange=timerange_ascii, \
#            probe=mmsnum,datatypes=['pos', 'vel', 'spinras', 'spindec']) # just to see
# 
# load FGM, EDP, L2 data in DMPA, DSL, GSE/GSM (DMPA and DSL are identical)
#
fgm_tvars = pyspedas.mms.fgm(trange=timerange_ascii,time_clip='True', probe=mmsnum) # default is survey
edp_tvars = pyspedas.mms.edp(trange=timerange_ascii,time_clip='True', probe=mmsnum) # default is survey
tnames(fgm_tvars) # check what you loaded
tnames(edp_tvars) # check what you loaded
edp_tvars_scpot = pyspedas.mms.edp(trange=timerange_ascii,time_clip='True', \
              datatype='scpot',probe=mmsnum) # default is survey
tnames(edp_tvars_scpot) # check what you loaded, use "scpot" (rel. to plasma), instead of psp (rel.to.probes)
#
mytvars2plot = ['fgm_b_dmpa_srvy_l2','edp_scpot_fast_l2','edp_dce_dsl_fast_l2', 'edp_dce_gse_fast_l2']
mytvars2plot = [mysc + '_' + item for item in mytvars2plot]
tplot(mytvars2plot,xsize=9,ysize=12) # also set the plot size accordingly
#
# go straight to GSM coord's here:
cotrans(name_in=myscu+'edp_dce_gse_fast_l2',name_out=myscu+'edp_dce_gsm_fast_l2',coord_in='gse',coord_out='gsm')
set_coords(myscu+'edp_dce_gsm_fast_l2','GSM')
# interpolate FGS data on EFS times in GSM:
tinterpol(myscu+'fgm_b_gsm_srvy_l2_bvec',myscu+'edp_dce_gsm_fast_l2',newname=myscu+'fgm_b_gsm_srvy_l2_int')
#
# compute V=ExB/B^2
# 
tcrossp(myscu+'edp_dce_gsm_fast_l2',myscu+'fgm_b_gsm_srvy_l2_int',newname='exb_temp')
tdotp(myscu+'fgm_b_gsm_srvy_l2_int',myscu+'fgm_b_gsm_srvy_l2_int',newname='bdotb_temp') # form B^2 as N array
bdotb_temp=deepcopy(get_data('bdotb_temp'))
# below will divide denominator B^2 (i.e., multiply numerator ExB) by 1000., to scale V from m/s to km/s
store_data('bdotb_temp',data={'x':bdotb_temp.times,\
                  'y': bdotb_temp.y[:,np.newaxis]*array([1, 1, 1])/1000.}) # Nx3
pytplot.divide('exb_temp','bdotb_temp',myscu+'Vexb_gsm_srvy_l2')
options(myscu+'Vexb_gsm_srvy_l2','colors',['b','g','r'])
options(myscu+'Vexb_gsm_srvy_l2','labels',['Vx ExB GSM','Vy ExB GSM','Vz ExB GSM'])
set_units(myscu+'Vexb_gsm_srvy_l2','km/s')
options(myscu+'Vexb_gsm_srvy_l2','ytitle','Vexb [km/s]')
#
# see the results
mytvars2plot = [myscu + item for item in ['fgm_b_gsm_srvy_l2','edp_dce_gsm_fast_l2','Vexb_gsm_srvy_l2']]
tlimit(time2zoom_more)   ; tplot(mytvars2plot,xsize=12,ysize=12) # also set the plot size accordingly
tlimit(full='True') ; tplot(mytvars2plot,xsize=12,ysize=12)
#
# load FPI data 
#
# load ground-pre-computed FPI L2 moment data and omni spectra
tvars_fpi=pyspedas.mms.fpi(trange=timerange_double, time_clip='True', probe=mmsnum)
tnames(tvars_fpi) # check what you loaded
print("Number of new tvars loaded:", len(tvars_fpi)) # print number of quantities loaded
#
# rotate velocities to GSM
cotrans(name_in=myscu+'dis_bulkv_gse_fast',name_out=myscu+'dis_bulkv_gsm_fast',coord_in='gse',coord_out='gsm')
set_coords(myscu+'dis_bulkv_gsm_fast','GSM')
cotrans(name_in=myscu+'des_bulkv_gse_fast',name_out=myscu+'des_bulkv_gsm_fast',coord_in='gse',coord_out='gsm')
set_coords(myscu+'des_bulkv_gsm_fast','GSM')
#
# plot fgm, edp, moms, spectra+scpot
store_data(myscu + 'dxs_Nie_fast',data=[myscu + 'dis_numberdensity_fast',myscu + 'des_numberdensity_fast'])
options(myscu + 'dxs_Nie_fast','colors',['b','r']) # this does not work
options(myscu + 'dxs_Nie_fast','ytitle','FPI_Ne,i')
store_data(myscu + 'dxs_Tparper_fast',data= [myscu +  item for item in \
              ['dis_temppara_fast','dis_tempperp_fast','des_temppara_fast','des_tempperp_fast']])
options(myscu + 'dxs_Tparper_fast','ylog','True')
options(myscu + 'dxs_Tparper_fast','yrange',[4.e2,2.e4])
options(myscu + 'dxs_Tparper_fast','ytitle','Tei_parper')
options(myscu + 'edp_scpot_fast_l2','ylog','True')
options(myscu + 'edp_scpot_fast_l2','thick',1.)
store_data(myscu + 'des_eflux_omni_fast_wscpot',data= [myscu +  item for item in \
         ['des_energyspectr_omni_fast','edp_scpot_fast_l2']])
options(myscu + 'des_eflux_omni_fast_wscpot','ytitle','DES_eflux') # does not work
options(myscu + 'des_eflux_omni_fast_wscpot','ysubtitle','[eV]') # does not work
mytvars2plot=[myscu + item for item in \
               ['fgm_b_gsm_srvy_l2','edp_dce_gsm_fast_l2','Vexb_gsm_srvy_l2',\
                'dxs_Nie_fast','dxs_Tparper_fast','dis_bulkv_gsm_fast','des_bulkv_gsm_fast',\
                'dis_energyspectr_omni_fast','des_eflux_omni_fast_wscpot']]
tplot_options('title','MMS3 Overview')
tlimit(full='True') ; tplot(mytvars2plot,xsize=12,ysize=12)
tplot(mytvars2plot,xsize=12,ysize=12,save_png='mms_ovr_BE_FPImoms')
#
############## END HERE ##########

# use peer_en_eflux spectra to plot along with spacecraft potential
#
# rescale spacecraft potential
#
myscpot=get_data(mysc+'_pxxm_pot')
myscpot_md=get_data(mysc+'_pxxm_pot',metadata=True)
myscpot_meta=myscpot_md
store_data(mysc+'_pxxm_pot2',data={'x':myscpot.times,'y':(1.+myscpot.y)*1.15},attr_dict=myscpot_meta)
# 
# include pot2 in spectra
options(mysc+'_peir_en_eflux','yrange',[5.,2.5e4])
options(mysc+'_peer_en_eflux','yrange',[5.,3.e4])
options(mysc+'_pxxm_pot2','yrange',[5.,3.e4])
options(mysc+'_pxxm_pot2','ylog',1)
options(mysc+'_peer_en_eflux','ylog',1)
store_data(mysc+'_peer_en_eflux_pot',data=[mysc+'_peer_en_eflux',mysc+'_pxxm_pot2']) # pseudovariable
options(mysc+'_peer_en_eflux_pot','yrange',[5.,3.e4])
options(mysc+'_peer_en_eflux_pot','ylog',1)
options(mysc+'_peir_en_eflux','ytitle','peir_en_eflux')
options(mysc+'_peer_en_eflux_pot','ytitle','peer_en_eflux_pot')

#
store_data(mysc+'_pexm_density',data=[mysc+'_peim_density',mysc+'_peem_density']) # pseudovariable
options(mysc+'_pexm_density','colors',['b','r'])
options(mysc+'_pexm_density','yrange',[0.05,5])
options(mysc+'_pexm_density','ytitle','Ni,e')
#
options(mysc+'_peim_velocity_gsm','ytitle','peim_Vgsm')
#
mytvars2plot =['fgs_gsm','efs_dot0_gsm','Vexb_dot0_gsm', \
       'pexm_density','peim_velocity_gsm','peir_en_eflux','peer_en_eflux_pot']
mytvars2plot = [mysc+'_'+ item for item in mytvars2plot]
tplot(mytvars2plot,xsize=10,ysize=12) # also set the plot size accordingly
nickssl commented 8 months ago

The file mms_overview_plot.py shows how to resolve some of these issues.

To avoid the multiple titles, use: options(panel01, "name", "")

For line colors, use: options(panel01, "Color", ["black", "green"])

The following seem to work correctly:

options(panel01, "name", "") options(panel01, "Color", ["black", "green"]) options(panel01, "legend_names", ["Themis AE", "Kyoto proxy AE"]) options(panel01, "ytitle", "AE Intex") options(panel01, "ysubtitle", "[nt]") options(panel01, "legend_location", "spedas")

nickssl commented 8 months ago

There are two related bugs:

The following produces a legend which is partially hidden, because the width of the plot is not computed correctly: options(panel01, "legend_location", "spedas") See also: https://github.com/spedas/pyspedas/issues/641

The second problem is that store_data crashes when it is used with two tplot variables and one of them does not exist. IDL handles this correctly. Example: store_data(panel01, data=[d1_1, d1_2]) This can probably be fixed while working on https://github.com/spedas/pyspedas/issues/654

jameswilburlewis commented 5 months ago

More concise code to reproduce the multiple plot titles:

import pyspedas
from pytplot import store_data,tplot,tplot_options
pyspedas.themis.state()
store_data('ps1',['thc_spin_initial_delta_phi','thc_spin_idpu_spinper'])
store_data('ps2',['thc_spin_initial_delta_phi','thc_spin_idpu_spinper'])
store_data('ps3',['thc_spin_initial_delta_phi','thc_spin_idpu_spinper'])
tplot_options('title','plot title')
plotvars=['thc_pos', 'ps1', 'thc_vel', 'ps2', 'thc_pos_gse', 'ps3']
tplot(plotvars)
jameswilburlewis commented 5 months ago

See also: https://github.com/spedas/pyspedas/issues/796 for problems combining line plots with different numbers of traces.

jameswilburlewis commented 5 months ago

The multiple plot titles are fixed now. The other problems have been filed as separate issues, so this one can be closed.