Advanced Plotting#
Here are some more sophisticated techniques for visualizing CEDAR data.
import numpy as np
import h5py
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import cartopy.crs as ccrs
import madrigalWeb.madrigalWeb
import os
RTI Plots of LP and AC Data#
Download files:
madrigalUrl='http://cedar.openmadrigal.org'
data = madrigalWeb.madrigalWeb.MadrigalData(madrigalUrl)
user_fullname = 'Student Example'
user_email = 'isr.summer.school@gmail.com'
user_affiliation= 'ISR Summer School 2024'
code = 61 # PFISR
year = 2024
month = 1
day = 8
hour1 = 7
minute1 = 1
hour2 = 13
# list of experiments inside a time period of a day
expList = data.getExperiments(code,year,month,day,hour1,minute1,0,year,month,day,hour2,0,0)
for exp in expList:
print(str(exp))
id: 100278619
realUrl: http://cedar.openmadrigal.org/showExperiment/?experiment_list=100278619
url: http://cedar.openmadrigal.org/madtoc/experiments4/2024/pfa/08jan24a
name: Themis36 - Auroral and convection measurements
siteid: 10
sitename: CEDAR
instcode: 61
instname: Poker Flat IS Radar
startyear: 2024
startmonth: 1
startday: 8
starthour: 7
startmin: 1
startsec: 4
endyear: 2024
endmonth: 1
endday: 8
endhour: 18
endmin: 0
endsec: 0
isLocal: True
madrigalUrl: http://cedar.openmadrigal.org/
PI: Asti Bhatt
PIEmail: asti.bhatt@sri.com
uttimestamp: 1709109883
access: 2
Madrigal version: 3.2
fileList = data.getExperimentFiles(expList[0].id)
for file0 in fileList:
print(os.path.basename(file0.name),'\t', file0.kindat, '\t',file0.kindatdesc)
pfa20240108.001_ac_nenotr_01min.001.h5 1000201 Ne From Power - Alternating Code (E-region) - 1 min
pfa20240108.001_ac_fit_01min.001.h5 2000201 Fitted - Alternating Code (E-region) - 1 min
pfa20240108.001_ac_nenotr_03min.001.h5 1000203 Ne From Power - Alternating Code (E-region) - 3 min
pfa20240108.001_ac_fit_03min.001.h5 2000203 Fitted - Alternating Code (E-region) - 3 min
pfa20240108.001_ac_nenotr_05min.001.h5 1000205 Ne From Power - Alternating Code (E-region) - 5 min
pfa20240108.001_ac_fit_05min.001.h5 2000205 Fitted - Alternating Code (E-region) - 5 min
pfa20240108.001_ac_nenotr_10min.001.h5 1000210 Ne From Power - Alternating Code (E-region) - 10 min
pfa20240108.001_ac_fit_10min.001.h5 2000210 Fitted - Alternating Code (E-region) - 10 min
pfa20240108.001_ac_nenotr_15min.001.h5 1000215 Ne From Power - Alternating Code (E-region) - 15 min
pfa20240108.001_ac_fit_15min.001.h5 2000215 Fitted - Alternating Code (E-region) - 15 min
pfa20240108.001_ac_nenotr_20min.001.h5 1000220 Ne From Power - Alternating Code (E-region) - 20 min
pfa20240108.001_ac_fit_20min.001.h5 2000220 Fitted - Alternating Code (E-region) - 20 min
pfa20240108.001_lp_nenotr_01min.001.h5 1000101 Ne From Power - Long Pulse (F-region) - 1 min
pfa20240108.001_lp_fit_01min.001.h5 2000101 Fitted - Long Pulse (F-region) - 1 min
pfa20240108.001_lp_nenotr_03min.001.h5 1000103 Ne From Power - Long Pulse (F-region) - 3 min
pfa20240108.001_lp_fit_03min.001.h5 2000103 Fitted - Long Pulse (F-region) - 3 min
pfa20240108.001_lp_nenotr_05min.001.h5 1000105 Ne From Power - Long Pulse (F-region) - 5 min
pfa20240108.001_lp_fit_05min.001.h5 2000105 Fitted - Long Pulse (F-region) - 5 min
pfa20240108.001_lp_nenotr_10min.001.h5 1000110 Ne From Power - Long Pulse (F-region) - 10 min
pfa20240108.001_lp_fit_10min.001.h5 2000110 Fitted - Long Pulse (F-region) - 10 min
pfa20240108.001_lp_nenotr_15min.001.h5 1000115 Ne From Power - Long Pulse (F-region) - 15 min
pfa20240108.001_lp_fit_15min.001.h5 2000115 Fitted - Long Pulse (F-region) - 15 min
pfa20240108.001_lp_nenotr_20min.001.h5 1000120 Ne From Power - Long Pulse (F-region) - 20 min
pfa20240108.001_lp_fit_20min.001.h5 2000120 Fitted - Long Pulse (F-region) - 20 min
pfa20240108.001_lp_vvels_01min.001.h5 3000101 Resolved Velocity - Long Pulse (F-region) - 1 min
pfa20240108.001_lp_vvels_03min.001.h5 3000103 Resolved Velocity - Long Pulse (F-region) - 3 min
pfa20240108.001_lp_vvels_05min.001.h5 3000105 Resolved Velocity - Long Pulse (F-region) - 5 min
# Download the file that we need to run these examples
os.makedirs('data', exist_ok=True)
filepath_lp= 'data/pfa20240108.001_lp_fit_01min.001.h5 '
filepath_ac= 'data/pfa20240108.001_ac_fit_01min.001.h5 '
if not os.path.exists(filepath_lp):
fileList = data.getExperimentFiles(expList[0].id)
for file0 in fileList:
if file0.kindatdesc == 'Fitted - Long Pulse (F-region) - 1 min':
file2download = file0.name
break
print('Downloading data file...')
file = data.downloadFile(file2download, filepath_lp,
user_fullname, user_email, user_affiliation,'hdf5')
print('...Done!')
else:
print(f"File {filepath_lp} already downloaded")
if not os.path.exists(filepath_ac):
fileList = data.getExperimentFiles(expList[0].id)
for file0 in fileList:
if file0.kindatdesc == 'Fitted - Alternating Code (E-region) - 1 min':
file2download = file0.name
break
print('Downloading data file...')
file = data.downloadFile(file2download, filepath_ac,
user_fullname, user_email, user_affiliation,'hdf5')
print('...Done!')
else:
print(f"File {filepath_ac} already downloaded")
Downloading data file...
...Done!
Downloading data file...
...Done!
It is possible to create RTI plots of LP and AC data together to visualize the full profile.
with h5py.File(filepath_lp, 'r') as h5:
elms = []
bmid = []
for key in h5['Data/Array Layout'].keys():
elms.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['elm'][0]))
bmid.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['beamid'][0]))
elms = np.array(elms)
bmid = np.array(bmid)
idx = np.argmax(elms)
max_beam = bmid[idx]
bidx = f'Array with beamid={max_beam} '
rng_lp = np.array(h5['Data/Array Layout'][bidx]['range'])
ne_lp = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['ne']).T
utime_lp = np.array(h5['Data/Array Layout'][bidx]['timestamps'])
time_lp = utime_lp.astype('datetime64[s]')
ne_lp = ne_lp[:,np.isfinite(rng_lp)]
rng_lp = rng_lp[np.isfinite(rng_lp)]
with h5py.File(filepath_ac, 'r') as h5:
elms = []
bmid = []
for key in h5['Data/Array Layout'].keys():
elms.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['elm'][0]))
bmid.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['beamid'][0]))
elms = np.array(elms)
bmid = np.array(bmid)
idx = np.argmax(elms)
max_beam = bmid[idx]
bidx = f'Array with beamid={max_beam} '
rng_ac = np.array(h5['Data/Array Layout'][bidx]['range'])
ne_ac = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['ne']).T
utime_ac = np.array(h5['Data/Array Layout'][bidx]['timestamps'])
azm = h5['Data/Array Layout'][bidx]['1D Parameters']['azm'][0]
elm = h5['Data/Array Layout'][bidx]['1D Parameters']['elm'][0]
beamid = h5['Data/Array Layout'][bidx]['1D Parameters']['beamid'][0]
time_ac = utime_ac.astype('datetime64[s]')
ne_ac = ne_ac[:,np.isfinite(rng_ac)]
rng_ac = rng_ac[np.isfinite(rng_ac)]
cutoff_rng = 150.*1000.
aidx_ac = np.argmin(np.abs(rng_ac-cutoff_rng))
aidx_lp = np.argmin(np.abs(rng_lp-cutoff_rng))
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(111)
c = ax.pcolormesh(time_ac, rng_ac[:aidx_ac], ne_ac[:,:aidx_ac].T, vmin=0., vmax=4.e11)
c = ax.pcolormesh(time_lp, rng_lp[aidx_lp:], ne_lp[:,aidx_lp:].T, vmin=0., vmax=4.e11)
ax.set_xlabel('Universal Time')
ax.set_ylabel('Range (m)')
ax.set_title('Beam Number: {:.0f} (az:{:.1f}, el:{:.1f})'.format(beamid, azm, elm))
fig.colorbar(c, label=r'Electron Density (m$^{-3}$)')
<matplotlib.colorbar.Colorbar at 0x7fd62fc89d80>
Plot all ISR Parameters#
Create a figure showing all four ISR parameters (Ne, Ti, Te, Vlos).
with h5py.File(filepath_lp, 'r') as h5:
elms = []
bmid = []
for key in h5['Data/Array Layout'].keys():
elms.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['elm'][0]))
bmid.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['beamid'][0]))
elms = np.array(elms)
bmid = np.array(bmid)
idx = np.argmax(elms)
max_beam = bmid[idx]
bidx = f'Array with beamid={max_beam} '
rng_lp = np.array(h5['Data/Array Layout'][bidx]['range'])
ne_lp = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['ne']).T
te_lp = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['te']).T
ti_lp = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['ti']).T
vlos_lp = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['vo']).T
utime_lp = np.array(h5['Data/Array Layout'][bidx]['timestamps'])
time_lp = utime_lp.astype('datetime64[s]')
ne_lp = ne_lp[:,np.isfinite(rng_lp)]
te_lp = te_lp[:,np.isfinite(rng_lp)]
ti_lp = ti_lp[:,np.isfinite(rng_lp)]
vlos_lp = vlos_lp[:,np.isfinite(rng_lp)]
rng_lp = rng_lp[np.isfinite(rng_lp)]
with h5py.File(filepath_ac, 'r') as h5:
elms = []
bmid = []
for key in h5['Data/Array Layout'].keys():
elms.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['elm'][0]))
bmid.append(np.array(h5['Data/Array Layout'][key]['1D Parameters']['beamid'][0]))
elms = np.array(elms)
bmid = np.array(bmid)
idx = np.argmax(elms)
max_beam = bmid[idx]
bidx = f'Array with beamid={max_beam} '
rng_ac = np.array(h5['Data/Array Layout'][bidx]['range'])
ne_ac = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['ne']).T
te_ac = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['te']).T
ti_ac = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['ti']).T
vlos_ac = np.array(h5['Data/Array Layout'][bidx]['2D Parameters']['vo']).T
utime_ac = np.array(h5['Data/Array Layout'][bidx]['timestamps'])
time_ac = utime_ac.astype('datetime64[s]')
ne_ac = ne_ac[:,np.isfinite(rng_ac)]
te_ac = te_ac[:,np.isfinite(rng_ac)]
ti_ac = ti_ac[:,np.isfinite(rng_ac)]
vlos_ac = vlos_ac[:,np.isfinite(rng_ac)]
rng_ac = rng_ac[np.isfinite(rng_ac)]
cutoff_rng = 150.*1000.
aidx_ac = np.argmin(np.abs(rng_ac-cutoff_rng))
aidx_lp = np.argmin(np.abs(rng_lp-cutoff_rng))
fig = plt.figure(figsize=(12,12))
gs = gridspec.GridSpec(4,1)
# Plot Electron Density
ax = fig.add_subplot(gs[0])
c = ax.pcolormesh(time_ac, rng_ac[:aidx_ac], ne_ac[:,:aidx_ac].T, vmin=0., vmax=4.e11, cmap='viridis')
c = ax.pcolormesh(time_lp, rng_lp[aidx_lp:], ne_lp[:,aidx_lp:].T, vmin=0., vmax=4.e11, cmap='viridis')
# ax.set_xlabel('Universal Time')
ax.set_ylabel('Range (m)')
fig.colorbar(c, label=r'Electron Density (m$^{-3}$)')
# Plot Ion Temperature
ax = fig.add_subplot(gs[1])
c = ax.pcolormesh(time_ac, rng_ac[:aidx_ac], ti_ac[:,:aidx_ac].T, vmin=0., vmax=3.e3, cmap='magma')
c = ax.pcolormesh(time_lp, rng_lp[aidx_lp:], ti_lp[:,aidx_lp:].T, vmin=0., vmax=3.e3, cmap='magma')
# ax.set_xlabel('Universal Time')
ax.set_ylabel('Range (m)')
fig.colorbar(c, label=r'Ion Temperature (K)')
# Plot Electron Temperature
ax = fig.add_subplot(gs[2])
c = ax.pcolormesh(time_ac, rng_ac[:aidx_ac], te_ac[:,:aidx_ac].T, vmin=0., vmax=5.e3, cmap='inferno')
c = ax.pcolormesh(time_lp, rng_lp[aidx_lp:], te_lp[:,aidx_lp:].T, vmin=0., vmax=5.e3, cmap='inferno')
# ax.set_xlabel('Universal Time')
ax.set_ylabel('Range (m)')
fig.colorbar(c, label=r'Electron Temperature (K)')
# Plot Line-of-Site Velocity
ax = fig.add_subplot(gs[3])
c = ax.pcolormesh(time_ac, rng_ac[:aidx_ac], vlos_ac[:,:aidx_ac].T, vmin=-500., vmax=500., cmap='bwr')
c = ax.pcolormesh(time_lp, rng_lp[aidx_lp:], vlos_lp[:,aidx_lp:].T, vmin=-500., vmax=500., cmap='bwr')
ax.set_xlabel('Universal Time')
ax.set_ylabel('Range (m)')
fig.colorbar(c, label=r'Line-of-Site Velocity (m/s)')
<matplotlib.colorbar.Colorbar at 0x7fd6304eded0>
