Skip to content

Commit

Permalink
functionalize plotting
Browse files Browse the repository at this point in the history
  • Loading branch information
scivision committed Mar 30, 2018
1 parent 99e94a2 commit eed36af
Show file tree
Hide file tree
Showing 4 changed files with 145 additions and 149 deletions.
31 changes: 5 additions & 26 deletions AltitudeProfile.py
@@ -1,38 +1,17 @@
#!/usr/bin/env python
""" Height Profile Example """
import pyiri2016
import pyiri2016 as iri
import pyiri2016.plots as piri
#
import numpy as np
from matplotlib.pyplot import figure, show
from matplotlib.pyplot import show

glat, glon = -11.95, -76.77

alt_km = np.arange(80,1000,20.)

iri = pyiri2016.IRI('2012-08-21T12', alt_km, glat, glon)
iono = iri.IRI('2012-08-21T12', alt_km, glat, glon)

fig = figure(figsize=(16,6))
axs = fig.subplots(1,2)

fig.suptitle(f'{str(iri.time[0].values)[:-13]}\n Glat, Glon: {iri.glat}, {iri.attrs["glon"]}')


pn = axs[0]
pn.plot(iri['ne'].squeeze(), iri.alt_km, label='N$_e$')
#pn.set_title(iri2016Obj.title1)
pn.set_xlabel('Density (m$^{-3}$)')
pn.set_ylabel('Altitude (km)')
pn.set_xscale('log')
pn.legend(loc='best')
pn.grid(True)

pn = axs[1]
pn.plot(iri['Ti'].squeeze(), iri.alt_km, label='T$_i$')
pn.plot(iri['Te'].squeeze(), iri.alt_km, label='T$_e$')
#pn.set_title(iri2016Obj.title2)
pn.set_xlabel('Temperature (K)')
pn.set_ylabel('Altitude (km)')
pn.legend(loc='best')
pn.grid(True)
piri.altprofile(iono)

show()
50 changes: 14 additions & 36 deletions LatitudeProfile.py
@@ -1,45 +1,23 @@
#!/usr/bin/env python
import pyiri2016
from numpy import arange
from matplotlib.pyplot import figure, show
import pyiri2016 as iri
import pyiri2016.plots as piri
from matplotlib.pyplot import show

""" Geog. Latitude Profile Example """

latlim = [-60, 60]
latstp = 2.
sim = pyiri2016.geoprofile(altkm=600, latlim=latlim, dlat=latstp, \
glon=-76.77, time='2004-01-01T17')

latbins = arange(latlim[0], latlim[1], latstp)
if __name__ == '__main__':
from argparse import ArgumentParser
p = ArgumentParser()
p.add_argument('-o','--outfn',help='write data to file')
p = p.parse_args()

latlim = [-60, 60]
latstp = 2.
iono = iri.geoprofile(altkm=600, latlim=latlim, dlat=latstp,
glon=-76.77, time='2004-01-01T17')

fig = figure(figsize=(8,12))
axs = fig.subplots(2,1, sharex=True)
piri.latprofile(iono)

pn = axs[0]

pn.plot(latbins, sim['NmF2'].squeeze(), label='N$_m$F$_2$')
pn.plot(latbins, sim['NmF1'].squeeze(), label='N$_m$F$_1$')
pn.plot(latbins, sim['NmE'].squeeze(), label='N$_m$E')
pn.set_title(str(sim.time[0].values)[:-13] + ' latitude'+str(latlim))
pn.set_xlim(latbins[[0, -1]])
pn.set_xlabel('Geog. Lat. ($^\circ$)')
pn.set_ylabel('(m$^{-3}$)')
pn.set_yscale('log')


pn = axs[1]
pn.plot(latbins, sim['hmF2'].squeeze(), label='h$_m$F$_2$')
pn.plot(latbins, sim['hmF1'].squeeze(), label='h$_m$F$_1$')
pn.plot(latbins, sim['hmE'].squeeze(), label='h$_m$E')
pn.set_xlim(latbins[[0, -1]])
pn.set_title(str(sim.time[0].values)[:-13] + ' latitude'+str(latlim))
pn.set_xlabel('Geog. Lat. ($^\circ$)')
pn.set_ylabel('(km)')

for a in axs:
a.legend(loc='best')
a.grid(True)


show()
show()
92 changes: 5 additions & 87 deletions TimeProfile.py
Expand Up @@ -2,13 +2,10 @@
""" Time Profile: IRI2016 """
import numpy as np
from datetime import timedelta
from matplotlib.pyplot import figure, show
try:
import ephem
except ImportError:
ephem = None
from matplotlib.pyplot import show
#
import pyiri2016
import pyiri2016 as iri
import pyiri2016.plots as piri


# %% user parameters
Expand All @@ -17,88 +14,9 @@
glat,glon = 65,-148
alt_km = np.arange(120, 180, 20)
# %% ru
sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
sim = iri.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
alt_km,glat,glon)

# %% Plots
Nplot=3

if Nplot>2:
fig = figure(figsize=(16,12))
axs = fig.subplots(3,1, sharex=True).ravel()
else:
fig = figure(figsize=(16,6))
axs = fig.subplots(1,2).ravel()

fig.suptitle(f'{str(sim.time[0].values)[:-13]} to {str(sim.time[-1].values)[:-13]}\n Glat, Glon: {sim.glat}, {sim.glon}')

ax = axs[0]
#NmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[2, :])
#NmE = sim.loc[4, :]
ax.plot(sim.time, sim['NmF2'].squeeze(), label='N$_m$F$_2$')
#ax.plot(sim.time, NmF1, label='N$_m$F$_1$')
#ax.plot(sim.time, NmE, label='N$_m$E')
ax.set_title('Maximum number densities vs. ionospheric layer')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(m$^{-3}$)')
ax.set_yscale('log')
ax.legend(loc='best')

ax = axs[1]
#hmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[3, :])
#hmE = sim.b[5, :]
ax.plot(sim.time, sim['hmF2'].squeeze(), label='h$_m$F$_2$')
#ax.plot(sim.time, hmF1, label='h$_m$F$_1$')
#ax.plot(sim.time, hmE, label='h$_m$E')
ax.set_title('Height of maximum density vs. ionospheric layer')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(km)')
ax.legend(loc='best')
# %%
if Nplot > 2:
ax = axs[2]

for a in sim.alt_km:
ax.plot(sim.time, sim['ne'].squeeze(), marker='.', label=f'{a.item()} km')
ax.set_xlabel('time UTC (hours)')
ax.set_ylabel('[m$^{-3}$]')
ax.set_title(f'$N_e$ vs. altitude and time')
ax.set_yscale('log')
ax.legend(loc='best')
# %%
if Nplot > 4:
ax = axs[4]
tec = sim.b[36, :]
ax.plot(sim.time, tec, label=r'TEC')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(m$^{-2}$)')
#ax.set_yscale('log')
ax.legend(loc='best')

ax = axs[5]
vy = sim.b[43, :]
ax.plot(sim.time, vy, label=r'V$_y$')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(m/s)')
ax.legend(loc='best')

for a in axs.ravel():
a.grid(True)


if __name__ == '__main__':
from argparse import ArgumentParser
p = ArgumentParser(description='IRI2016 time profile plot')
p.add_argument('-t','--trange',help='START STOP STEP (hours) time [UTC]',nargs=3,
default=('2012-08-21','2012-08-22',0.25))
p.add_argument('--alt',help='START STOP STEP altitude [km]',type=float, nargs=3,default=(120,180,20))
p.add_argument('-c','--latlon',help='geodetic coordinates of simulation',
type=float,default=(65,-147.5))
p.add_argument('--f107',type=float,default=200.)
p.add_argument('--f107a', type=float,default=200.)
p.add_argument('--ap', type=int, default=4)
p.add_argument('--species',help='species to plot',nargs='+',default=('ne'))
p = p.parse_args()

piri.timeprofile(sim)

show()
121 changes: 121 additions & 0 deletions pyiri2016/plots.py
@@ -0,0 +1,121 @@
import xarray
from matplotlib.pyplot import figure

def timeprofile(iono:xarray.Dataset):
# %% Plots
Nplot=3

if Nplot>2:
fig = figure(figsize=(16,12))
axs = fig.subplots(3,1, sharex=True).ravel()
else:
fig = figure(figsize=(16,6))
axs = fig.subplots(1,2).ravel()

fig.suptitle(f'{str(iono.time[0].values)[:-13]} to {str(iono.time[-1].values)[:-13]}\n Glat, Glon: {iono.glat}, {iono.glon}')

ax = axs[0]

ax.plot(iono.time, iono['NmF2'].squeeze(), label='N$_m$F$_2$')
ax.plot(iono.time, iono['NmF1'].squeeze(), label='N$_m$F$_1$')
ax.plot(iono.time, iono['NmE'].squeeze(), label='N$_m$E')
ax.set_title('Maximum number densities vs. ionospheric layer')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(m$^{-3}$)')
ax.set_yscale('log')
ax.legend(loc='best')

ax = axs[1]
ax.plot(iono.time, iono['hmF2'].squeeze(), label='h$_m$F$_2$')
ax.plot(iono.time, iono['hmF1'].squeeze(), label='h$_m$F$_1$')
ax.plot(iono.time, iono['hmE'].squeeze(), label='h$_m$E')
ax.set_title('Height of maximum density vs. ionospheric layer')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(km)')
ax.legend(loc='best')
# %%
if Nplot > 2:
ax = axs[2]

for a in iono.alt_km:
ax.plot(iono.time, iono['ne'].squeeze(), marker='.', label=f'{a.item()} km')
ax.set_xlabel('time UTC (hours)')
ax.set_ylabel('[m$^{-3}$]')
ax.set_title(f'$N_e$ vs. altitude and time')
ax.set_yscale('log')
ax.legend(loc='best')
# %%
if Nplot > 4:
ax = axs[4]
tec = iono.b[36, :]
ax.plot(iono.time, tec, label=r'TEC')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(m$^{-2}$)')
#ax.set_yscale('log')
ax.legend(loc='best')

ax = axs[5]
vy = iono.b[43, :]
ax.plot(iono.time, vy, label=r'V$_y$')
ax.set_xlabel('Hour (UT)')
ax.set_ylabel('(m/s)')
ax.legend(loc='best')

for a in axs.ravel():
a.grid(True)


def altprofile(iono:xarray.Dataset):
fig = figure(figsize=(16,6))
axs = fig.subplots(1,2)

fig.suptitle(f'{str(iono.time[0].values)[:-13]}\n Glat, Glon: {iono.glat}, {iono.attrs["glon"]}')


pn = axs[0]
pn.plot(iono['ne'].squeeze(), iono.alt_km, label='N$_e$')
#pn.set_title(iri2016Obj.title1)
pn.set_xlabel('Density (m$^{-3}$)')
pn.set_ylabel('Altitude (km)')
pn.set_xscale('log')
pn.legend(loc='best')
pn.grid(True)

pn = axs[1]
pn.plot(iono['Ti'].squeeze(), iono.alt_km, label='T$_i$')
pn.plot(iono['Te'].squeeze(), iono.alt_km, label='T$_e$')
#pn.set_title(iri2016Obj.title2)
pn.set_xlabel('Temperature (K)')
pn.set_ylabel('Altitude (km)')
pn.legend(loc='best')
pn.grid(True)

def latprofile(iono:xarray.Dataset):

fig = figure(figsize=(8,12))
axs = fig.subplots(2,1, sharex=True)

ax = axs[0]

ax.plot(iono.lat, iono['NmF2'].squeeze(), label='N$_m$F$_2$')
ax.plot(iono.lat, iono['NmF1'].squeeze(), label='N$_m$F$_1$')
ax.plot(iono.lat, iono['NmE'].squeeze(), label='N$_m$E')
ax.set_title(str(iono.time[0].values)[:-13] + f' latitude {iono.lat[[0, -1]].values}')
#ax.set_xlim(iono.lat[[0, -1]])
ax.set_xlabel('Geog. Lat. ($^\circ$)')
ax.set_ylabel('(m$^{-3}$)')
ax.set_yscale('log')


ax = axs[1]
ax.plot(iono.lat, iono['hmF2'].squeeze(), label='h$_m$F$_2$')
ax.plot(iono.lat, iono['hmF1'].squeeze(), label='h$_m$F$_1$')
ax.plot(iono.lat, iono['hmE'].squeeze(), label='h$_m$E')
ax.set_xlim(iono.lat[[0, -1]])
ax.set_title(str(iono.time[0].values)[:-13] + f' latitude {iono.lat[[0, -1]].values}')
ax.set_xlabel('Geog. Lat. ($^\circ$)')
ax.set_ylabel('(km)')

for a in axs:
a.legend(loc='best')
a.grid(True)

0 comments on commit eed36af

Please sign in to comment.