first time density profile plot workign

space-physics · Mar 5, 2018 · e4dc4e3 · e4dc4e3
1 parent adc1d9f
commit e4dc4e3
Show file tree

Hide file tree

Showing 2 changed files with 56 additions and 50 deletions.
diff --git a/examples/TimeProfile.py b/examples/TimeProfile.py
@@ -1,27 +1,29 @@
 #!/usr/bin/env python
 """ Time Profile: IRI2016 """
-import pyiri2016
-from numpy import arange
+import numpy as np
 from datetime import timedelta
 from matplotlib.pyplot import figure, show
 try:
     import ephem
 except ImportError:
     ephem = None
+#
+import pyiri2016
+
+
 # %% user parameters
 hrlim = [0, 24] # [0,24] does whole day(s)
 hrstp = 0.25 # time step [decimal hours]
 #lat = -11.95; lon = -76.77
 glat = 0
 glon = 0
-alt_km = arange(120, 180, 10)
+alt_km = np.arange(120, 180, 10)
 # %% ru
 sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
                             alt_km,glat,glon)
 
 # %% Plots
-Nplot=2
-
+Nplot=4
 
 if Nplot>2:
     fig = figure(figsize=(16,12))
@@ -30,58 +32,59 @@
     fig = figure(figsize=(16,6))
     axs = fig.subplots(1,2).ravel()
 
-pn = axs[0]
+fig.suptitle(f'{np.datetime_as_string(sim.time[0])[:-13]} to {np.datetime_as_string(sim.time[-1])[:-13]}\n Glat, Glon: {sim.attrs["glat"]}, {sim.attrs["glon"]}')
+
+ax = axs[0]
 #NmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[2, :])
 #NmE = sim.loc[4, :]
-pn.plot(sim.time, sim.attrs['NmF2'], label='N$_m$F$_2$')
-#pn.plot(sim.time, NmF1, label='N$_m$F$_1$')
-#pn.plot(sim.time, NmE, label='N$_m$E')
-pn.set_title(sim.title1)
-pn.set_xlabel('Hour (UT)')
-pn.set_ylabel('(m$^{-3}$)')
-pn.set_yscale('log')
-pn.legend(loc='best')
-
-pn = axs[1]
-hmF2 = sim.b[1, :]
-hmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[3, :])
-hmE = sim.b[5, :]
-pn.plot(sim.time, hmF2, label='h$_m$F$_2$')
-pn.plot(sim.time, hmF1, label='h$_m$F$_1$')
-pn.plot(sim.time, hmE, label='h$_m$E')
-pn.set_title(sim.title2)
-pn.set_xlabel('Hour (UT)')
-pn.set_ylabel('(km)')
-pn.legend(loc='best')
+ax.plot(sim.time, sim.attrs['NmF2'], 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.attrs['hmF2'], 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:
-    pn = axs[2]
-
-    for alt in alt_km:
-        sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
-                                    alt,glat,glon)
+    ax = axs[2]
 
-        pn.plot(sim.time, sim.loc[:,'ne'], marker='.', label=alt)
-    pn.set_xlabel('time UTC (hours)')
-    pn.set_ylabel('[m$^{-3}$]')
-    pn.set_title(f'$N_e$ vs. altitude and time')
-    pn.legend(loc='best')
+    sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
+                                alt_km, glat,glon)
 
+    for a in sim.alt_km:
+        ax.plot(sim.time, sim.loc[:,a,'ne'], 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.legend(loc='best')
+# %%
 if Nplot > 4:
-    pn = axs[4]
+    ax = axs[4]
     tec = sim.b[36, :]
-    pn.plot(sim.time, tec, label=r'TEC')
-    pn.set_xlabel('Hour (UT)')
-    pn.set_ylabel('(m$^{-2}$)')
-    #pn.set_yscale('log')
-    pn.legend(loc='best')
+    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')
 
-    pn = axs[5]
+    ax = axs[5]
     vy = sim.b[43, :]
-    pn.plot(sim.time, vy, label=r'V$_y$')
-    pn.set_xlabel('Hour (UT)')
-    pn.set_ylabel('(m/s)')
-    pn.legend(loc='best')
+    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)

diff --git a/pyiri2016/__init__.py b/pyiri2016/__init__.py
@@ -118,7 +118,8 @@ def IRI( time, altkm, glat, glon, ap=None, f107=None, ssn=None, var=None):
 #     #------------------------------------------------------------------------------
 
     mmdd = 100*time.month + time.day               # month and dom (MMDD)
-
+    # hour + 25 denotes UTC time
+    dhour = (time.hour + 25) + time.minute/60.
 # %% more inputs
     jmag = 0            #  0: geographic; 1: geomagnetic
   #  iut = 0             #  0: for LT;     1: for UT
@@ -130,10 +131,10 @@ def IRI( time, altkm, glat, glon, ap=None, f107=None, ssn=None, var=None):
 #        a, b = iriwebg(jmag, jf, glat, glon, int(time.year), mmdd, iut, time.hour,
 #            height, h_tec_max, ivar, ivbeg, ivend, ivstp, addinp, self.iriDataFolder)
 
-    # hour + 25 denotes UTC time
+
 
     outf,oarr = iri2016.iri_sub(jf, jmag, glat, glon,
-                                time.year, mmdd, time.hour+25, altkm,
+                                time.year, mmdd, dhour, altkm,
                                 proot/'data/')
 
 # %% collect output
@@ -169,6 +170,8 @@ def timeprofile(tlim:tuple, dt:timedelta,
 
     T = datetimerange(tlim[0], tlim[1], dt)
 
+    altkm = np.atleast_1d(altkm)
+
     iono = xarray.DataArray(np.empty((len(T),altkm.size,9)),
                             coords={'time':T,'alt_km':altkm, 'sim':simout},
                             dims=['time','alt_km','sim'],