Preparing files to upload to pypi, predstorm_l5 now runs with persistence model (pickled NOAA real-time files) when STEREO-A is not available.

Author: R Bailey

predstorm/data.py

@@ -222,7 +222,7 @@ def __len__(self):
 
 
     def __str__(self):
-        """Print string describing object."""
+        """Return string describing object."""
 
         ostr = "Length of data:\t\t{}\n".format(len(self))
         ostr += "Keys in data:\t\t{}\n".format(self.vars)
@@ -1111,7 +1111,7 @@ def make_kp_prediction(self):
     def get_state(self):
         """Finds state of wind and fills self.state attribute."""
 
-        print("Coming soon.")
+        logger.info("Coming soon.")
 
     # -----------------------------------------------------------------------------------
     # Data archiving
@@ -1120,7 +1120,7 @@ def get_state(self):
     def archive(self):
         """Make archive of long-term data."""
 
-        print("Not yet implemented.")
+        logger.info("Not yet implemented.")
 
 
 class PositionData():
@@ -1907,7 +1907,6 @@ def get_omni_data(starttime=None, endtime=None, filepath='', download=False, dld
     #check how many rows exist in this file
     f=open(filepath)
     dataset= len(f.readlines())
-    #print(dataset)
     #global Variables
     spot=np.zeros(dataset) 
     btot=np.zeros(dataset) #floating points
@@ -1939,7 +1938,7 @@ def get_omni_data(starttime=None, endtime=None, filepath='', download=False, dld
     with open('data/omni2_all_years.dat') as f:
         for line in f:
             line = line.split() # to deal with blank 
-            #print line #41 is Dst index, in nT
+            #41 is Dst index, in nT
             dst[j]=line[40]
             ae[j]=line[41]
             kp[j]=line[38]
@@ -2475,7 +2474,6 @@ def merge_Data(satdata1, satdata2, keys=None):
 
 def getpositions(filename):  
     pos=scipy.io.readsav(filename)  
-    print
     return pos
 
 
@@ -2505,15 +2503,12 @@ def converttime():
     http://matplotlib.org/examples/pylab_examples/date_demo2.html
     """
 
-    print('convert time start')
     for index in range(0,dataset):
         #first to datetimeobject 
         timedum=datetime(int(year[index]), 1, 1) + timedelta(day[index] - 1) +timedelta(hours=hour[index])
         #then to matlibplot dateformat:
         times1[index] = matplotlib.dates.date2num(timedum)
-        #print time
-        #print year[index], day[index], hour[index]
-    print('convert time done')   #for time conversion
+    logger.info('convert time done')   #for time conversion
 
 
 def round_to_hour(dt):

predstorm/plot.py

@@ -141,9 +141,14 @@ def plot_solarwind_and_dst_prediction(DSCOVR_data, STEREOA_data, DST_data, DSTPR
 
     # STEREO-A minute resolution data with timeshift
     plt.plot_date(stam['time'][sta_index_future], stam['btot'][sta_index_future],
-                  '-', c=c_sta, linewidth=lw, label='B STEREO-A')
-    plt.plot_date(stam['time'][sta_index_future], stam['bn'][sta_index_future],
-                  '-', c=c_sta, alpha=0.5, linewidth=lw, label='Bn RTN STEREO-A')
+                  '-', c=c_sta, linewidth=lw, label='B {}'.format(stam.source))
+    if 'stereo' in stam.source.lower():
+        plt.plot_date(stam['time'][sta_index_future], stam['bn'][sta_index_future],
+                      '-', c=c_sta, alpha=0.5, linewidth=lw, label='Bn RTN {}'.format(stam.source))
+    else:
+        plt.plot_date(stam['time'][sta_index_future], stam['bz'][sta_index_future],
+                      '-', c=c_sta, alpha=0.5, linewidth=lw, label='Bz GSM {}'.format(stam.source))
+
 
     # Indicate 0 level for Bz
     plt.plot_date([plotstart,plotend], [0,0],'--k', alpha=0.5, linewidth=1)
@@ -157,8 +162,8 @@ def plot_solarwind_and_dst_prediction(DSCOVR_data, STEREOA_data, DST_data, DSTPR
 
     if 'stereo' in stam.source.lower():
         pred_source = 'STEREO-Ahead Beacon'
-    elif 'omni' in stam.source.lower():
-        pred_source = '27-day SW-Recurrence Model (OMNI)'
+    elif 'dscovr' in stam.source.lower() or 'noaa' in stam.source.lower():
+        pred_source = '27-day SW-Recurrence Model (NOAA)'
     plt.title('L1 real time solar wind from NOAA SWPC for '+ datetime.strftime(timestamp, "%Y-%m-%d %H:%M")+ ' UT & {}'.format(pred_source), fontsize=fs_title)
 
     # SUBPLOT 2: Solar wind speed
@@ -212,17 +217,17 @@ def plot_solarwind_and_dst_prediction(DSCOVR_data, STEREOA_data, DST_data, DSTPR
     dstplotmin = np.nanmin(np.concatenate((dst['dst'], dst_pred['dst'])))-20
 
     if dstplotmin > -100:       # Low activity (normal)
-        plt.ylim([-100,np.nanmax(dst['dst'])+20])
+        plt.ylim([-100, dstplotmax + 30])
     else:                       # High activity
         plt.ylim([dstplotmin, dstplotmax])
 
     # Plot predicted Dst
     if not verification_mode:
         plt.plot_date(dst_pred['time'], dst_pred['dst'], '-', c=c_sta_dst, label=dst_label, markersize=3, linewidth=1)
         # Add generic error bars of +/-15 nT:
-        error=15
+        error=8
         plt.fill_between(dst_pred['time'], dst_pred['dst']-error, dst_pred['dst']+error, alpha=0.2,
-                         label='Error for high speed streams')
+                         label='prediction accuracy +/- 8 nT')
     else:
         #load saved data l prefix is for loaded - WARNING This will crash if called right now
         [timenowb, sta_ptime, sta_vr, sta_btime, sta_btot, sta_br,sta_bt, sta_bn, rbtime_num, rbtot, rbzgsm, rptime_num, rpv, rpn, lrdst_time, lrdst, lcom_time, ldst_burton, ldst_obrien,ldst_temerin_li]=pickle.load(open(verify_filename,'rb') )
@@ -240,10 +245,10 @@ def plot_solarwind_and_dst_prediction(DSCOVR_data, STEREOA_data, DST_data, DSTPR
 
         # Plot solar wind density:
         avg_newell_coupling = newell_coupling.get_weighted_average('ec')
-        plt.plot_date(newell_coupling['time'], avg_newell_coupling/4421., '-', color=c_ec, label='Newell coupling',linewidth=1.5)
+        plt.plot_date(newell_coupling['time'], avg_newell_coupling/4421., '-', color=c_ec, label='Newell coupling 4h weighted mean',linewidth=1.5)
         plt.ylabel('Newell Coupling / 4421\n$\mathregular{[(km/s)^{4/3} nT^{2/3}]}$',fontsize=fs_ylabel)
         # For y limits check where the maximum and minimum are for DSCOVR and STEREO taken together:
-        plt.ylim([0,np.nanmax(avg_newell_coupling/4421.)*1.1])
+        plt.ylim([0,np.nanmax(avg_newell_coupling/4421.)*1.2])
 
         # Indicate level of interest (Ec/4421 = 1.0)
         plt.plot_date([plotstart,plotend], [1,1],'--k', alpha=0.5, linewidth=1)
@@ -288,7 +293,7 @@ def plot_solarwind_and_dst_prediction(DSCOVR_data, STEREOA_data, DST_data, DSTPR
     logger.info('Plot saved as png:\n'+ filename)
 
 
-def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, timestamp=None, past_days=7, future_days=7, outfile='predstorm_science.png', **kwargs):
+def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, timestamp=None, past_days=7, future_days=7, plot_step=20, outfile='predstorm_science.png', **kwargs):
     """
     Plots solar wind variables, past from DSCOVR and future/predicted from STEREO-A.
     Total B-field and Bz (top), solar wind speed (second), particle density (third)
@@ -355,6 +360,7 @@ def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, t
     timeutc = mdates.date2num(timestamp)
 
     n_plots = 3
+    plst = plot_step
 
     plotstart = timeutc - past_days
     plotend = timeutc + future_days
@@ -365,10 +371,10 @@ def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, t
     axes.append(ax1)
 
     # Total B-field and Bz (DSCOVR)
-    plt.plot_date(dism['time'], dism['btot'],'-', c='black', label='B', linewidth=lw)
-    plt.plot_date(dism['time'], dism['bx'],'-', c='teal', label='Bx', linewidth=lw)
-    plt.plot_date(dism['time'], dism['by'],'-', c='orange', label='By', linewidth=lw)
-    plt.plot_date(dism['time'], dism['bz'],'-', c='purple', label='Bz', linewidth=lw)
+    plt.plot_date(dism['time'][::plst], dism['btot'][::plst],'-', c='black', label='B', linewidth=lw)
+    plt.plot_date(dism['time'][::plst], dism['bx'][::plst],'-', c='teal', label='Bx', linewidth=lw)
+    plt.plot_date(dism['time'][::plst], dism['by'][::plst],'-', c='orange', label='By', linewidth=lw)
+    plt.plot_date(dism['time'][::plst], dism['bz'][::plst],'-', c='purple', label='Bz', linewidth=lw)
 
     # STEREO-A minute resolution data with timeshift
     plt.plot_date(stam['time'][sta_index_future], stam['btot'][sta_index_future],
@@ -392,8 +398,8 @@ def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, t
 
     if 'stereo' in stam.source.lower():
         pred_source = 'STEREO-Ahead Beacon'
-    elif 'omni' in stam.source.lower():
-        pred_source = '27-day SW-Recurrence Model (OMNI)'
+    elif 'dscovr' in stam.source.lower() or 'noaa' in stam.source.lower():
+        pred_source = '27-day SW-Recurrence Model (NOAA)'
     plt.title('L1 real time solar wind from NOAA SWPC for '+ datetime.strftime(timestamp, "%Y-%m-%d %H:%M")+ ' UT & {}'.format(pred_source), fontsize=fs_title)
 
     # SUBPLOT 2: Solar wind speed
@@ -402,7 +408,7 @@ def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, t
     axes.append(ax2)
 
     # Plot solar wind speed (DSCOVR):
-    plt.plot_date(dism['time'], dism['speed'],'-', c='black', label='speed',linewidth=lw)
+    plt.plot_date(dism['time'][::plst], dism['speed'][::plst],'-', c='black', label='speed',linewidth=lw)
     plt.ylabel('Speed $\mathregular{[km \\ s^{-1}]}$', fontsize=fs_ylabel)
 
     # Plot STEREO-A data with timeshift and savgol filter
@@ -425,7 +431,7 @@ def plot_solarwind_science(DSCOVR_data, STEREOA_data, verification_mode=False, t
     axes.append(ax3)
 
     # Plot solar wind density:
-    plt.plot_date(dism['time'], dism['density'],'-k', label='density',linewidth=lw)
+    plt.plot_date(dism['time'][::plst], dism['density'][::plst],'-k', label='density',linewidth=lw)
     plt.ylabel('Density $\mathregular{[ccm^{-3}]}$',fontsize=fs_ylabel)
     # For y limits check where the maximum and minimum are for DSCOVR and STEREO taken together:
     plt.ylim([0,np.nanmax(np.nanmax(np.concatenate((dism['density'],stam['density'][sta_index_future])))+10)])

predstorm/predict.py

@@ -412,7 +412,6 @@ def _jit_calc_dst_temerin_li_2002(time, btot, bx, by, bz, speed, speedx, density
         bp = (by[i]**2 + bx[i]**2)**0.5  
         #contains t1, but in cos and sin 
         dh = bp*np.cos(np.arctan2(bx[i],by[i])+6.10) * (3.59e-2 * np.cos(2*np.pi*t1/yearli + 0.04) - 2.18e-2*np.sin(2*np.pi*t1-1.60))
-        #print(i, bx[i], by[i], bz[i])
         theta_li = -(np.arccos(-bz[i]/bt)-np.pi)/2
         exx = 1e-3 * abs(speedx[i]) * bt * np.sin(theta_li)**6.1
         #t1 and dt are in days

predstorm_l5.py

@@ -173,6 +173,7 @@ def main():
             rec_new = ps.merge_Data(rec_27days, rec_dism)
             rec_new.cut(starttime=num2date(timenow-365))  # don't keep more than 1 year
             rec_new.source = 'NOAA-L1'
+            logger.info("Pickled real-time data contains data from {} till {}".format(num2date(rec_new['time'][0]), num2date(rec_new['time'][-1])))
             with open(rec_model_path, 'wb') as f:
                 pickle.dump(rec_new, f)
     elif run_mode == 'historic':
@@ -181,7 +182,6 @@ def main():
             logger.info("Using historic mode with current DSCOVR data and {} timestamp".format(timestamp))
             dism = ps.get_dscovr_realtime_data()
             dism = dism.cut(endtime=timestamp)
-            dis = dism.make_hourly_data()
         else:
             logger.info("Using historic mode with archive DSCOVR data")
             dism = ps.get_dscovr_data(starttime=timestamp-timedelta(days=plot_past_days+1),
@@ -190,6 +190,7 @@ def main():
         timenow = timeutc
     dism.interp_nans()
     dis = dism.make_hourly_data()
+    sw_past_min = dism
     sw_past = dis
 
     timeutcstr = datetime.strftime(timestamp, tstr_format)
@@ -204,35 +205,36 @@ def main():
     #------------------------ (1b) Get real-time STEREO-A beacon data -----------------------
 
     logger.info("(2) Getting STEREO-A data...")
+    use_recurrence_model = False
     if run_mode == 'normal':
-        stam = ps.get_stereo_beacon_data(starttime=startread, endtime=timestamp-timedelta(minutes=1))
+        try:
+            stam = ps.get_stereo_beacon_data(starttime=startread, endtime=timestamp-timedelta(minutes=1))
+            stam.load_positions()
+            sta_details = stam.return_position_details(timestamp)
+            if stam.h['PlasmaDataIntegrity'] == 0: # very low quality data
+                use_recurrence_model = True
+        except Exception as e:
+            logger.info("STEREO-A read failed for reason: {}".format(e))
+            use_recurrence_model = True
     elif run_mode == 'historic':
         lag_L1, lag_r = ps.get_time_lag_wrt_earth(timestamp=timestamp, satname='STEREO-A')
         est_timelag = lag_L1 + lag_r
         stam = ps.get_stereo_beacon_data(starttime=timestamp-timedelta(days=plot_future_days+est_timelag),
                                          endtime=timestamp+timedelta(days=2))
+        stam.load_positions()
+        sta_details = stam.return_position_details(timestamp)
 
-    logger.info('UTC time of last datapoint in STEREO-A beacon data:')
-    logger.info('\t{}'.format(num2date(stam['time'][-1])))
-    logger.info('Time lag in minutes: {}'.format(int(round((timeutc-stam['time'][-1])*24*60))))
-
-    # Prepare data:
-    stam.load_positions()
-    sta_details = stam.return_position_details(timestamp)
-
-    use_recurrence_model = False
-    if stam.h['PlasmaDataIntegrity'] == 0: # very low quality data
+    if run_mode == 'normal' and use_recurrence_model:
         use_recurrence_model = True
         logger.info("STEREO-A plasma data is missing/corrupted, using 27-day recurrence model for plasma data instead.")
         rec_start = timestamp - timedelta(days=27)
         rec_end = timestamp - timedelta(days=27-plot_future_days)
-        sw_future = ps.get_omni_data_new(starttime=rec_start, endtime=rec_end)
-        sw_future['time'] = sw_future['time'] + 27.
-        sw_future.h['DataSource'] += ' t+27days'
-        sw_future.source += '+27days'
-        # Remove OMNI indices, don't need them:
-        sw_future['kp'], sw_future['ae'], sw_future['dst'] = 0., 0., 0.
-        sw_future.vars = [v for v in sw_future.vars if v not in ['kp', 'ae', 'dst']]
+        with open(rec_model_path, 'rb') as f:
+            sw_future_min = pickle.load(f)
+        sw_future_min.cut(starttime=rec_start, endtime=rec_end)
+        sw_future_min['time'] += 27.
+        sw_future_min.h['DataSource'] += ' t+27days'
+        sw_future_min.source += '+27days'
 
     #------------------------- (1c) Load NOAA Dst for comparison ----------------------------
 
@@ -252,39 +254,38 @@ def main():
 
     #========================== (2) PREDICTION CALCULATIONS ==================================
 
-    logger.info("\n-------------------------\nL5-to-L1 MAPPING\n-------------------------")
-    #------------------------ (2a) Corrections to time-shifted STEREO-A data ----------------
-    logger.info("Doing corrections to STEREO-A data...")
-
-    logger.info("(1) Shift time at STEREO-A according to solar wind rotation")
-    stam.shift_time_to_L1()
-    logger.info("STA-to-L1 adjusted time of last datapoint in STEREO-A:")
-    logger.info("\t{}".format(num2date(stam['time'][-1])))
-
-    logger.info("(2) Make correction for difference in heliocentric distance")
-    stam.shift_wind_to_L1()
-
-    logger.info("(3) Conversion from RTN to GSE and then to GSM as if STEREO was on Sun-Earth line")
-    stam['bx'], stam['by'], stam['bz'] = stam['br'], -stam['bt'], stam['bn']    # RTN to quasi-GSE
-    stam.convert_GSE_to_GSM()
-
-    use_recurrence_model = False
     if not use_recurrence_model:
+        logger.info('UTC time of last datapoint in STEREO-A beacon data:')
+        logger.info('\t{}'.format(num2date(stam['time'][-1])))
+        logger.info('Time lag in minutes: {}'.format(int(round((timeutc-stam['time'][-1])*24*60))))
+        logger.info("\n-------------------------\nL5-to-L1 MAPPING\n-------------------------")
+        #------------------------ (2a) Corrections to time-shifted STEREO-A data ----------------
+        logger.info("Doing corrections to STEREO-A data...")
+
+        logger.info("(1) Shift time at STEREO-A according to solar wind rotation")
+        stam.shift_time_to_L1()
+        logger.info("STA-to-L1 adjusted time of last datapoint in STEREO-A:")
+        logger.info("\t{}".format(num2date(stam['time'][-1])))
+
+        logger.info("(2) Make correction for difference in heliocentric distance")
+        stam.shift_wind_to_L1()
+
+        logger.info("(3) Conversion from RTN to GSE and then to GSM as if STEREO was on Sun-Earth line")
+        stam['bx'], stam['by'], stam['bz'] = stam['br'], -stam['bt'], stam['bn']    # RTN to quasi-GSE
+        stam.convert_GSE_to_GSM()
         sw_future = stam.make_hourly_data()
         sw_future_min = stam
     else:
-        sta = stam.interp_to_time(sw_future['time'])
-        sw_future.vars += ['br', 'bt', 'bn']
-        for bvar in ['bx', 'by', 'bz', 'btot', 'br', 'bt', 'bn']:
-            sw_future[bvar] = sta[bvar] 
-        sw_future.interp_nans()
-        sw_future_min = copy.deepcopy(sw_future)
+        sw_future_min.vars += ['br', 'bt', 'bn']
+        sw_future_min['br'], sw_future_min['bt'], sw_future_min['bn'] = sw_future_min['bx'], sw_future_min['by'], sw_future_min['bz']
+        sw_future_min.interp_nans()
+        sw_future = sw_future_min.make_hourly_data()
 
     #------------------- (2b) COMBINE DSCOVR and time-shifted STEREO-A data -----------------
 
     sw_merged = ps.merge_Data(sw_past, sw_future)
     try:
-        sw_merged_min = ps.merge_Data(dism, stam)
+        sw_merged_min = ps.merge_Data(sw_past_min, sw_future_min)
         savemindata = True
     except:
         logger.warning("No minute data available.")
@@ -334,15 +335,15 @@ def main():
     logger.info("\n-------------------------\nPLOTTING\n-------------------------")
     # ********************************************************************
     logger.info("Creating output plot...")
-    plot_solarwind_and_dst_prediction([dism, sw_past], [sw_future_min, sw_future], 
+    plot_solarwind_and_dst_prediction([sw_past_min, sw_past], [sw_future_min, sw_future], 
                                       dst, dst_pred,
                                       newell_coupling=newell_coupling,
                                       past_days=plot_past_days,
                                       future_days=plot_future_days,
                                       dst_label=dst_label,
                                       timestamp=timestamp)
     plt.close()
-    plot_solarwind_science([dism, sw_past], [sw_future_min, sw_future], 
+    plot_solarwind_science([sw_past_min, sw_past], [sw_future_min, sw_future], 
                                       timestamp=timestamp)
     # ********************************************************************
 
@@ -373,7 +374,8 @@ def main():
 
     #----------------------------- (4b) CALCULATE FORECAST RESULTS --------------------------
 
-    logger.info("\n\nSATELLITE POSITION DETAILS\n--------------------------\n"+sta_details)
+    if not use_recurrence_model:
+        logger.info("\n\nSATELLITE POSITION DETAILS\n--------------------------\n"+sta_details)
 
     future_times = np.where(sw_merged['time'] > date2num(timestamp))
     past_times = np.where(sw_merged['time'] < date2num(timestamp))

setup.cfg

@@ -0,0 +1,2 @@
+[metadata]
+description-file = README.md