Finished review work on code

686e249c · Leonie Pick · da98c4d7 · 686e249c · 686e249c · 686e249c
Commit 686e249c authored Aug 07, 2019 by Leonie Pick
5 changed files
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,7 @@ Dump/

 # Don't track the following files
 Input.nc
+Input_old.nc
 Readme.md

 # Don't track all files with these endings:

--- a/.ipynb_checkpoints/ClassifyStorms-checkpoint.ipynb
+++ b/.ipynb_checkpoints/ClassifyStorms-checkpoint.ipynb
--- a/ClassifyStorms.ipynb
+++ b/ClassifyStorms.ipynb
--- a/Modules.py
+++ b/Modules.py
@@ -55,7 +55,7 @@ def Search_TargetEvents(HMC, HMC11y, HMC5d, dHMC, HTime, DTime, grid, Plot, Save
        elif (IndexDiff[i] > 1):                                        # difference is greater than 1

            if j >= jmax:
-                jmax= j
+                jmax = j

            if j > mini_range:
                
@@ -135,7 +135,7 @@ def Search_TargetEvents(HMC, HMC11y, HMC5d, dHMC, HTime, DTime, grid, Plot, Save
    return Index_thres1, StormIndices, IndexMin1, IndexMin2
 ###
 ###
-def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Reference, Time, Date, YearStart, Plot, Save):
+def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Dst, SYM, SYMMinutes, Reference, Time, Date, YearStart, Mode, Plot, Save):

    YearsIndex = np.where(Time[:,0] >= YearStart)[0]
    TrClass = Reference[:,0]
@@ -172,7 +172,7 @@ def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Reference, T
                                     Time[YearsIndex,:], Date[YearsIndex], grid[i,:], Plot, Save)

        if Plot == True:
-            pl.Selection(Time[YearsIndex,:], Date[YearsIndex], HMC[YearsIndex], Index_thres, Storms_a, Storms_b, Storms, Kp_all, KpHours_all, Save)
+            pl.Selection(Time[YearsIndex,:], Date[YearsIndex], HMC[YearsIndex], Index_thres, Storms_a, Storms_b, Storms, Kp_all, KpHours_all, Dst[YearsIndex], SYM, SYMMinutes, Save)

        #Found_CIRs = np.where(np.in1d(TrTimeIndex[np.logical_or(TrClass==0,TrClass==2)],Storms+YearsIndex[0]))[0]
        #Found_CMEs = np.where(np.in1d(TrTimeIndex[TrClass==1],Storms+YearsIndex[0]))[0]
@@ -208,14 +208,14 @@ def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Reference, T
    ###
    TargetResults = {'No. target events': len(StormsWin[SelectStr]),'Fraction [%]':len(StormsWin[SelectStr])*100/len(HMC[YearsIndex])}
    TrainingResults = [[len(CIRs), len(CMEs), len(CIRs)+len(CMEs)],[len(FoundStorms['CIRs_'+SelectStr]),len(FoundStorms['CMEs_'+SelectStr]),len(FoundStorms['CIRs_'+SelectStr])+len(FoundStorms['CMEs_'+SelectStr])]]
-    SelectionParams = {r'$P_n [nT]$':grid[:,0],r'$p_{sc}$':grid[:,1], r'$\Delta t$ [h]':grid[:,2], r' $Hs$ [nT]':-grid[:,3]}
+    SelectionParams = {r'$n$':grid[:,0],r'$p_{sc}$':grid[:,1], r'$\Delta t$ [h]':grid[:,2]+1, r' $Hs$ [nT]':-grid[:,3]}
    display(pd.DataFrame(data=SelectionParams))
    display(pd.DataFrame(data=TargetResults,index=['Selection']))
    display(pd.DataFrame(data=TrainingResults,columns=['No. CIRs','No. CMEs', 'Total'], index=['Reference set', 'Training set']))
    ###
    
    if Plot == True:
-        pl.IndexDist(Time,YearsIndex,StormsWin[SelectStr],TrTimeIndex,Kp_all,KpHours_all,HMC,Save)
+        pl.IndexDist(Time,YearsIndex,StormsWin[SelectStr],TrTimeIndex,Kp_all,KpHours_all,HMC,Dst,Mode,Save)

    return StormsWin, TrFound, FoundStorms
 ###

--- a/Plots.py
+++ b/Plots.py
 from Imports import *
 ###
-def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin1,IndexMin2,Kp,KpHours,Save):
+def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin1,IndexMin2,Kp,KpHours,Dst,SYM,SYMMinutes,Save):
    
    fig, ax = plt.subplots(4,1,facecolor='w', edgecolor='k')
    fig.set_size_inches([6.69, 8.86], forward=True)
@@ -19,25 +19,36 @@ def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin1,IndexMin2,Kp,KpH

    dot_size=20
    # HMC
-    ax[0].plot(Time[start:end,4],HMC[start:end],color='black',linewidth=0.75,label='HMC',zorder=0)
-    ax[1].plot(Time[start:end,4],HMC[start:end],color='black',linewidth=0.75,zorder=0)
-    ax[2].plot(Time[start:end,4],HMC[start:end],color='black',linewidth=0.75,zorder=1)
+    ax[0].plot(Time[start:end,4],HMC[start:end],color='black',linewidth=0.75,label='HMC',zorder=1)
+    ax[1].plot(Time[start:end,4],HMC[start:end],color='black',linewidth=0.75,zorder=1)
+    ax[2].plot(Time[start:end,4],HMC[start:end],color='black',linewidth=0.75,zorder=2)

    # HMC selection
-    ax[0].plot(Time[start:end,4],Index_thres1[start:end],linestyle='--',color='black',label=r'$Hl_t\in$ ['+str(np.around(max(Index_thres1[start:end]),2))+' nT, '+str(np.around(min(Index_thres1[start:end]),2))+' nT]',zorder=2)
+    ax[0].plot(Time[start:end,4],Index_thres1[start:end],linewidth=0.75,linestyle='--',color='black',label=r'$Hl_t\in$ ['+str(np.around(max(Index_thres1[start:end]),2))+' nT, '+str(np.around(min(Index_thres1[start:end]),2))+' nT]',zorder=2)
    ax[0].scatter(Time[StormIndices,4],HMC[StormIndices],color='black',s=10,zorder=4)
    ax[1].scatter(Time[np.asarray(IndexMin1,dtype=int),4],HMC[np.asarray(IndexMin1,dtype=int)],color='black',s=dot_size,zorder=2)
    ax[2].scatter(Time[np.asarray(IndexMin2,dtype=int),4],HMC[np.asarray(IndexMin2,dtype=int)],color='black',s=dot_size,zorder=3)
    
    # Dst
    Dst_color='dimgray'
+    ax[0].plot(Time[start:end,4],Dst[start:end],color=Dst_color,linewidth=0.75,label='Dst',zorder=0)
+    ax[1].plot(Time[start:end,4],Dst[start:end],color=Dst_color,linewidth=0.75,zorder=0)
+    ax[2].plot(Time[start:end,4],Dst[start:end],color=Dst_color,linewidth=0.75,zorder=1)
+
+    # Dst storm scale
+    G1Dst= Rectangle((Time[start,4],-30),len(Time[start:end,4]),-20,edgecolor='yellow',facecolor='yellow',zorder=0,alpha=0.3,label='weak')
+    G2Dst = Rectangle((Time[start,4],-50),len(Time[start:end,4]),-50,edgecolor='gold',facecolor='gold',zorder=0,alpha=0.3,label='moderate')
+    #G3Dst = Rectangle((Time[start,4],-100),len(Time[start:end,4]),-100,edgecolor='orange',facecolor='orange',zorder=0,alpha=0.3,label='strong')
+    #G4Dst = Rectangle((Time[start,4],-200),len(Time[start:end,4]),-150,edgecolor='red',facecolor='red',zorder=0,alpha=0.3,label='severe')
+    #G5Dst = Rectangle((Time[start,4],-300),len(Time[start:end,4]),-600,edgecolor='maroon',facecolor='maroon',zorder=0,alpha=0.3,label='great')
+    ax[2].add_patch(G1Dst);ax[2].add_patch(G2Dst)#;ax[2].add_patch(G3Dst);ax[2].add_patch(G4Dst);ax[2].add_patch(G5Dst)

    # Kp 
    startKPT = 1981.665925; endKPT = 1981.747774
    startKP = np.where(KpHours==startKPT)[0][0]
    endKP = np.where(KpHours==endKPT)[0][0]

-    # Kp
+    # Kp storm scale
    Kp_step = 0.333333
    ax3b = ax[3].twinx()
    ax3b.bar(Time[start:end,4][::3],Kp[startKP:endKP+1],3/(24*365),align='edge',edgecolor=Dst_color,color=Dst_color,zorder=1)
@@ -51,6 +62,7 @@ def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin1,IndexMin2,Kp,KpH
    ax[3].patch.set_visible(False)

    # Sym-H
+    ax[3].plot(SYMMinutes,SYM,color='black',linewidth=.75)

    # Axes
    x_unique, x_indices = np.unique(Time[start:end,2],return_index=True)
@@ -89,62 +101,40 @@ def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin1,IndexMin2,Kp,KpH
    plt.show()

    ####################################################################
-    # SYM = np.loadtxt(Path_Input+'SYM_ASY_09_1981.txt',skiprows=15,usecols=[3,4,5,6])
-    # SYM_Time = np.linspace(HTime['dyear'][start],HTime['dyear'][end],len(SYM[30:len(SYM)-29,3]))
-    
-    # # Dst
-    # Dst_color='dimgray'
-    # ax[0].plot(HTime['dyear'][start:end],Dst[start:end],color=Dst_color,linewidth=0.75,label='Dst',zorder=1)
-    # ax[1].plot(HTime['dyear'][start:end],Dst[start:end],color=Dst_color,zorder=1)
-    # ax[2].plot(HTime['dyear'][start:end],Dst[start:end],color=Dst_color,zorder=2)
    # # Dst selection
    # ax[0].plot(HTime['dyear'][start:end],np.ones(len(HTime['dyear'][start:end]))-30,linestyle='--',color=Dst_color,label=r'$Hl = $-30 nT',zorder=3)
    # ax[0].scatter(HTime['dyear'][np.asarray(Dst_StormIndices)],Dst[np.asarray(Dst_StormIndices)],color=Dst_color,s=10,zorder=5)
    # ax[1].scatter(HTime['dyear'][np.asarray(Dst_IndexMin1)],Dst[np.asarray(Dst_IndexMin1)],color=Dst_color,s=dot_size,zorder=3)
    # ax[2].scatter(HTime['dyear'][np.asarray(Dst_IndexMin2)],Dst[np.asarray(Dst_IndexMin2)],color=Dst_color,s=dot_size,zorder=4)
-
-    # # Dst storm scale
-    # G1Dst= Rectangle((HTime['dyear'][start],-30),len(HTime['dyear'][start:end]),-20,edgecolor='yellow',facecolor='yellow',zorder=0,alpha=0.3,label='weak')
-    # G2Dst = Rectangle((HTime['dyear'][start],-50),len(HTime['dyear'][start:end]),-50,edgecolor='gold',facecolor='gold',zorder=0,alpha=0.3,label='moderate')
-    # #G3Dst = Rectangle((HTime['dyear'][start],-100),len(HTime['dyear'][start:end]),-100,edgecolor='orange',facecolor='orange',zorder=0,alpha=0.3,label='strong')
-    # #G4Dst = Rectangle((HTime['dyear'][start],-200),len(HTime['dyear'][start:end]),-150,edgecolor='red',facecolor='red',zorder=0,alpha=0.3,label='severe')
-    # #G5Dst = Rectangle((HTime['dyear'][start],-300),len(HTime['dyear'][start:end]),-600,edgecolor='maroon',facecolor='maroon',zorder=0,alpha=0.3,label='great')
-    # ax[2].add_patch(G1Dst);ax[2].add_patch(G2Dst)#;ax[2].add_patch(G3Dst);ax[2].add_patch(G4Dst);ax[2].add_patch(G5Dst);
-
-    # # Sym-H (only for this specific time)
-    # ax[3].plot(SYM_Time,SYM[30:len(SYM)-29,3],color='black',alpha=0.75)
-    # #ax[3].plot(SYM_Time,-SYM[30:len(SYM)-29,0],color='blue')
   
 ###
 ###
-def IndexDist(Time,YearsIndex,Storms,RefStorms,Kp_all,KpHours_all,HMC,Save):
+def IndexDist(Time,YearsIndex,Storms,RefStorms,Kp_all,KpHours_all,HMC,Dst,mode,Save):
    
-    mode = 'Kp'
    TimeDecYear = Time[:,4]

    # Overlap of storm times with Kp
    Kp_indices = np.where(np.logical_and(KpHours_all >= min(TimeDecYear[YearsIndex]),KpHours_all < max(TimeDecYear[YearsIndex])))[0]
    Kp_red = Kp_all[Kp_indices]; KpHours_red = KpHours_all[Kp_indices]
    Kp_step = 0.3333333
-    #Kp_range = np.around(np.arange(0,9,Kp_step),3)
    Kp_range = np.around(np.arange(0-Kp_step/2,9+Kp_step/2,Kp_step),3)

    Storms1932 = np.where(TimeDecYear[Storms] >= 1932.000171)[0]
-    HMC_step = 10#7.78
+    Storms1957 = np.where(TimeDecYear[Storms] >= 1957.000171)[0]
+    HMC_step = 10
    HMC_range = np.arange(0,480+HMC_step,HMC_step)
    
-    if mode == 'HMC':
+    if mode == 'HMC' or mode == 'Dst':
        Index_range = HMC_range
        step = HMC_step
    else:
        Index_range = Kp_range
        step = Kp_step 
        
-    years = np.arange(1932,2016,1)
-    Percentages = np.zeros((len(years),len(Index_range)))
+    years = np.arange(1930,2016,1)
    Hist = np.zeros((len(years),len(Index_range)-1))
    RHist =  np.zeros((len(years),len(Index_range)-1))
-    KpValuesAll = []
+    ValuesAll = []
    RValuesAll = []

    for l in range(len(years)):
@@ -186,33 +176,44 @@ def IndexDist(Time,YearsIndex,Storms,RefStorms,Kp_all,KpHours_all,HMC,Save):
        #KpValues = Kp_hours[KpStorms,1]
        #RKpValues = Kp_hours[RKpStorms,1]

-        HMCValues = -HMC[Storms[StormsBlock]] 
-        RHMCValues = -HMC[RefStorms[RBlockHMC]]
+        if mode == 'HMC':
+            Values = -HMC[Storms[StormsBlock]] 
+            RValues = -HMC[RefStorms[RBlockHMC]]
+        elif mode == 'Dst':
+            Values = -Dst[Storms[StormsBlock]] 
+            RValues = -Dst[RefStorms[RBlockHMC]]
+        else:
+            Values = KpValues
+            RValues = RKpValues
+            
+        ValuesAll.extend(Values.tolist())
+        RValuesAll.extend(RValues.tolist())
+        Hist[l,:], bin_edges = np.histogram(Values, bins = Index_range)
+        RHist[l,:], bin_edges = np.histogram(RValues, bins = Index_range)

+    # HMC / Kp / Dst
+    if mode == 'HMC' or mode == 'Dst':
        if mode == 'HMC':
-            KpValuesAll.extend(HMCValues.tolist())
-            RValuesAll.extend(RHMCValues.tolist())
-            Hist[l,:], bin_edges = np.histogram(HMCValues, bins = Index_range)
-            RHist[l,:], bin_edges = np.histogram(RHMCValues, bins = Index_range)
+            SumHist = np.sum(Hist*100/len(Storms),axis=0)
+            SumHist1957 = np.sum(Hist[27:,:]*100/len(Storms1957),axis=0)
+            bound25_1957 = np.percentile(np.asarray(ValuesAll[27:]).flatten(),25)
+            bound75_1957 = np.percentile(np.asarray(ValuesAll[27:]).flatten(),75)
        else:
-            KpValuesAll.extend(KpValues.tolist())
-            RValuesAll.extend(RKpValues.tolist())
-            Hist[l,:], bin_edges = np.histogram(KpValues, bins = Index_range)
-            RHist[l,:], bin_edges = np.histogram(RKpValues, bins = Index_range) 
-
-    # HMC / Kp
-    SumHist = np.sum(Hist*100/len(Storms1932),axis=0)
-    CumSum = np.cumsum(SumHist)
-    CumSum2 = np.cumsum(np.insert(SumHist[0:len(SumHist)-1],0,0))
-    bound25 = np.percentile(np.asarray(KpValuesAll).flatten(),25)
-    bound75 = np.percentile(np.asarray(KpValuesAll).flatten(),75)
-    bound50 = np.percentile(np.asarray(KpValuesAll).flatten(),50)
-    #print(bound25,bound50, bound75)
+            SumHist = np.sum(Hist*100/len(Storms1957),axis=0)
+    else:
+        SumHist = np.sum(Hist*100/len(Storms1932),axis=0)
+
+    bound25 = np.percentile(np.asarray(ValuesAll).flatten(),25)
+    bound75 = np.percentile(np.asarray(ValuesAll).flatten(),75)
+    bound50 = np.percentile(np.asarray(ValuesAll).flatten(),50)
+    print(bound25,bound50, bound75)

    # Reference HMC / Kp
    SumRHist = np.sum(RHist*100/len(RefStorms),axis=0)
    bound25_R = np.nanpercentile(np.asarray(RValuesAll).flatten(),25)
    bound75_R = np.nanpercentile(np.asarray(RValuesAll).flatten(),75)
+    bound50_R = np.nanpercentile(np.asarray(RValuesAll).flatten(),50)
+    print(bound25_R, bound50_R, bound75_R)

    fig = plt.figure()
    fig.set_size_inches([4, 8], forward=True)
@@ -261,8 +262,12 @@ def IndexDist(Time,YearsIndex,Storms,RefStorms,Kp_all,KpHours_all,HMC,Save):
    axins.set_yticklabels(['0','5','10','15','20'])

    ####################################################################
-
    # Dst-scale
+    G1Dst = Rectangle(((30-Index_range[0])/step,0),20/step,100,edgecolor='yellow',facecolor='yellow',zorder=0,alpha=0.3,label='weak')
+    G2Dst = Rectangle(((50-Index_range[0])/step,0),50/step,100,edgecolor='gold',facecolor='gold',zorder=0,alpha=0.3,label='moderate')
+    G3Dst = Rectangle(((100-Index_range[0])/step,0),100/step,100,edgecolor='orange',facecolor='orange',zorder=0,alpha=0.3,label='strong')
+    G4Dst = Rectangle(((200-Index_range[0])/step,0),150/step,100,edgecolor='red',facecolor='red',zorder=0,alpha=0.3,label='severe')
+    G5Dst = Rectangle(((350-Index_range[0])/step,0),600/step,100,edgecolor='maroon',facecolor='maroon',zorder=0,alpha=0.3)#,label='great')

    # NOAA G-scale
    G1 = Rectangle(((5-1.5*step-Index_range[0])/step-.5,0),3,100,edgecolor='yellow',facecolor='yellow',zorder=0,alpha=0.3,label='G1')
@@ -271,18 +276,23 @@ def IndexDist(Time,YearsIndex,Storms,RefStorms,Kp_all,KpHours_all,HMC,Save):
    G4 = Rectangle(((8-1.5*step-Index_range[0])/step-.5,0),4,100,edgecolor='red',facecolor='red',zorder=0,alpha=0.3,label='G4')
    G5 = Rectangle(((9-0.5*step-Index_range[0])/step-.5,0),1,100,edgecolor='maroon',facecolor='maroon',zorder=0,alpha=0.3,label='G5')
  
-    if mode == 'HMC': 
+    if mode == 'HMC' or mode == 'Dst': 

        ax2.bar(np.arange(0,len(Index_range)-1,1),SumHist,width=1.0,align='edge',color='darkgray',label='Target')
        ax2.bar(np.arange(0,len(Index_range)-1,1),SumRHist,width=1.0,align='edge',edgecolor='midnightblue',fill=False,label='Reference')
-        #ax2.bar(np.arange(0,len(Index_range)-1,1),SumHist2,width=1.0,align='edge',edgecolor='maroon',fill=False,label='Target 1957')
        ax2.axvline((bound25-Index_range[0])/step,color='black',linewidth=.75,linestyle='--'); ax2.axvline((bound75-Index_range[0])/step,color='black',linewidth=.75,linestyle='--',label='Q1, Q3')
        ax2.axvline((bound25_R-Index_range[0])/step,color='midnightblue',linewidth=.75,linestyle='--'); ax2.axvline((bound75_R-Index_range[0])/step,color='midnightblue',linewidth=.75,linestyle='--')    
-        #ax2.axvline((bound25_2-Index_range[0])/step,color='maroon',linestyle='--'); ax2.axvline((bound75_2-Index_range[0])/step,color='maroon',linestyle='--')
-        
+
+        if mode == 'HMC':
+            ax2.bar(np.arange(0,len(Index_range)-1,1),SumHist1957,width=1.0,align='edge',edgecolor='maroon',fill=False,label='Target 1957')
+            ax2.axvline((bound25_1957-Index_range[0])/step,color='maroon',linewidth=.75,linestyle='--'); ax2.axvline((bound75_1957-Index_range[0])/step,color='maroon',linewidth=.75,linestyle='--')
+            ax2.set_xlabel('-HMC [nT]')
+        else:
+            ax2.add_patch(G1Dst);ax2.add_patch(G2Dst);ax2.add_patch(G3Dst);ax2.add_patch(G4Dst);ax2.add_patch(G5Dst)
+            ax2.set_xlabel('-Dst [nT]')
+   
        ax2.set_xlim([0,xlimit])
        ax2.set_xticklabels(Index_range[::3], fontsize=9)
-        ax2.set_xlabel('-HMC [nT]')
        ax2.legend(loc=1,ncol=1,frameon=False,fontsize=9)

    else:
@@ -304,8 +314,8 @@ def IndexDist(Time,YearsIndex,Storms,RefStorms,Kp_all,KpHours_all,HMC,Save):
    ax2.set_ylabel('Occurrence [%]')

    if Save == True:
-        fig.savefig('./Dump/Fig/development/KpDist_test.pdf',format='pdf',dpi=200,transparent=True)
-        fig.savefig('./Dump/Fig/development/KpDist_test.png',format='png',dpi=200,transparent=True)
+        fig.savefig('./Dump/Fig/development/'+mode+'Dist.pdf',format='pdf',dpi=200,transparent=True)
+        fig.savefig('./Dump/Fig/development/'+mode+'Dist.png',format='png',dpi=200,transparent=True)
    plt.show()
 ###
 ###