plot_and_save_single_sat.py
plot_and_save_single_sat.pyΒΆ
NB: Do not use the download button above. Go to the repository on GitHub to access the files.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 26 11:09:38 2021
@author: blagau
"""
#SAVES INPUT AND OUTPUT TO ASCII FILES
#====================================
str_trange = dtime_beg.replace('-','')[:8]+'_'+ \
dtime_beg.replace(':','')[11:17] + '_'+dtime_end.replace(':','')[11:17]
fname_out = 'FAC_'+rez+'_sw'+sat[0]+'_atmuso_'+ str_trange +'.dat'
fname_fig = 'FAC_'+rez+'_sw'+sat[0]+'_atmuso_'+ str_trange +'.pdf'
bad_ang = np.less(np.abs(np.cos(ang_BR*np.pi/180.)), np.cos(np.deg2rad(60)))
if rez == 'LR':
info_points = '\nres.: LR missing points: '+str(len(ind_gaps)) + ' ' +\
'low-latitude points: ' +str(len(np.where(bad_ang)[0]))
else:
info_points = '\nres.: HR low-latitude points: ' +str(len(np.where(bad_ang)[0]))
add_text = 'Time interval: ' + dtime_beg[:10]+' ' + \
dtime_beg[11:19] + ' - '+dtime_end[11:19] + info_points
if tincl is not None:
add_text = add_text + '\nInclination interval: ' + \
np.datetime_as_string(tincl[0])[11:19] + ' - ' + \
np.datetime_as_string(tincl[1])[11:19]
# exports the results
out_df = j_df.copy()
text_trj = '# time [YYYY-mm-ddTHH:MM:SS.f],\t Rmid_X [km],\t Rmid_Y,\t\
Rmid_Z,\t Jfac [microA/m^2],\t IRC,\t errJfac,\t errJirc,\t '
if use_filter:
out_df = pd.concat([j_df[['Rmid X','Rmid_Y','Rmid Z','FAC','IRC','FAC_er','IRC_er']],\
jflt_df[['FAC_flt','IRC_flt','FAC_flt_er','IRC_flt_er','ang_BR','incl']]], axis=1)
text_trj = text_trj + 'Jfac_flt,\t Jirc_flt,\t errJfac_flt,\t errJirc_flt,\t '
text_header = text_trj + 'ang_BR [deg],\t incl \n'
with open(fname_out, 'w') as file:
file.write('# Swarm sat: ' + str(sat[0]) + '\n')
file.write('# Model: ' + Bmodel + '\n')
if rez == 'LR':
file.write('# Number of missing data points in L1b file: ' + str(len(ind_gaps)) + '\n')
if len(ind_gaps) > 0:
file.write('#' + '; '.join(timegaps.strftime('%H:%M:%S').values) + '\n')
file.write(text_header)
out_df.to_csv(fname_out, mode='a', sep=",", na_rep = 'NaN', date_format='%Y-%m-%dT%H:%M:%S.%f',\
float_format='%15.4f', header=False)
#PLOTS THE RESULTS
#====================================
# creates fig and axes objects
fig_size = (8.27,11.69)
fig = plt.figure(figsize=fig_size, frameon=True)
fig.suptitle('FAC density estimate with single-satellite method\n', size='xx-large',\
weight = 'bold', y=0.99)
plt.figtext(0.21, 0.96, add_text, fontsize = 'x-large', \
va='top', ha = 'left')
xle, xri, ybo, yto = 0.125, 0.92, 0.095, 0.105 # plot margins
ratp = np.array([1, 1, 1, 1, 0.5, 0.5]) #relative hight of each panel
hsep = 0.005 # vspace between panels
nrp = len(ratp)
hun = (1-yto - ybo - (nrp-1)*hsep)/ratp.sum()
yle = np.zeros(nrp) # y left for each panel
yri = np.zeros(nrp)
for ii in range(nrp):
yle[ii] = (1 - yto) - ratp[: ii+1].sum()*hun - ii*hsep
yri[ii] = yle[ii] + hun*ratp[ii]
# creates axes for each panel
ax = [0] * nrp
for ii in range(nrp):
ax[ii] = fig.add_axes([xle, yle[ii], xri-xle, hun*ratp[ii]])
ax[ii].set_xlim(pd.Timestamp(dtime_beg), pd.Timestamp(dtime_end))
for ii in range(nrp -1):
ax[ii].set_xticklabels([])
#Plots time-series quantities
ax[0].plot(dB_df)
ax[0].set_ylabel('$dB_{GEO}$ sw'+sat[0]+'\n[nT]', linespacing=1.7)
ax[0].legend(['dB_X', 'dB_Y', 'dB_Z' ], loc='upper right')
ax[1].plot(j_df['FAC'], label='FAC')
if use_filter:
ax[1].plot(jflt_df['FAC_flt'], label='FAC_flt')
ax[1].get_shared_x_axes().join(ax[0], ax[1], ax[2], ax[3], ax[4], ax[5])
ax[1].set_ylabel(r'$J_{FAC}$'+'\n'+r'$[\mu A/m^2]$', linespacing=1.7)
ax[1].legend(loc='upper right')
ax[2].plot(j_df['IRC'], label='IRC')
if use_filter:
ax[2].plot(jflt_df['IRC_flt'], label='IRC_flt')
ax[2].get_shared_x_axes().join(ax[0], ax[1], ax[2], ax[3], ax[4], ax[5])
ax[2].set_ylabel(r'$J_{IRC}$'+'\n'+r'$[\mu A/m^2]$', linespacing=1.7)
ax[2].legend(loc='upper right')
ax[3].plot(j_df['FAC'], label='FAC')
ax[3].plot(FAC_L2['FAC_L2'], label='FAC_L2')
ax[3].get_shared_x_axes().join(ax[0], ax[1], ax[2], ax[3], ax[4], ax[5])
ax[3].set_ylabel(r'$J_{FAC}$'+'\n'+r'$[\mu A/m^2]$', linespacing=1.7)
ax[3].legend(loc='upper right')
ax[4].plot(j_df['ang_BR'])
ax[4].get_shared_x_axes().join(ax[0], ax[1], ax[2], ax[3], ax[4], ax[5])
ax[4].set_ylabel('ang. BR'+'\n'+r'$[deg]$', linespacing=1.7)
ax[4].legend()
ax[5].plot(j_df['incl'])
ax[5].get_shared_x_axes().join(ax[0], ax[1], ax[2], ax[3], ax[4], ax[5])
ax[5].set_ylabel('incl.'+'\n'+r'$[deg]$', linespacing=1.7)
ax[5].legend()
if tincl is not None:
for ii in range(nrp):
ax[ii].axvline(tincl[0], ls='--', c='k')
ax[ii].axvline(tincl[1], ls='--', c='k')
latc = Rsph[:,0]
lonc = Rsph[:,1]
locx = ax[nrp-1].get_xticks()
latc_ipl = np.round(np.interp(locx, mdt.date2num(dat_df.index), \
latc), decimals=2).astype('str')
lonc_ipl = np.round(np.interp(locx, mdt.date2num(dat_df.index), \
lonc), decimals=2).astype('str')
qdlat_ipl = np.round(np.interp(locx, mdt.date2num(dat_df.index), \
dat_df['QDLat']), decimals=2).astype('str')
qdlon_ipl = np.round(np.interp(locx, mdt.date2num(dat_df.index), \
dat_df['QDLon']), decimals=2).astype('str')
mlt_ipl = np.round(np.interp(locx, mdt.date2num(dat_df.index), \
dat_df['MLT']), decimals=1).astype('str')
lab_fin = ['']*len(locx)
for ii in range(len(locx)):
lab_ini = mdt.num2date(locx[ii]).strftime('%H:%M:%S')
lab_fin[ii] = lab_ini + '\n' +latc_ipl[ii] + '\n' + lonc_ipl[ii]+ \
'\n'+ qdlat_ipl[ii] + '\n' +qdlon_ipl[ii] + '\n' + mlt_ipl[ii]
ax[nrp-1].set_xticklabels(lab_fin)
plt.figtext(0.01, 0.01, 'Time\nLat\nLon\nQLat\nQLon\nMLT')
plt.show()
fig.savefig(fname_fig)