Module TIEGCM_Statistics.Plotter_AltProfComparison

Creates altitude profiles plots of the median of two datasets in order to compare them. Reads two data structures with results of statistical calculations and creates an altitude profiles plot where the two datasets are plotted together.
It creates one subfigure per Magnetic Local Time - Kp index combination.
Each subfigure contains the two altitude profiles of the 50th percentile (median) of each dataset.
In case the second dataset is not provided it automatically reads the data from the Tromso EISCAT radar.
X axis represents the variable studied and Y axis represents the altitude.

Expand source code
"""
Creates altitude profiles plots of the median of two datasets in order to compare them.
Reads two data structures with results of statistical calculations and creates an altitude profiles plot where the two datasets are plotted together.  
It creates one subfigure per Magnetic Local Time - Kp index combination.  
Each subfigure contains the two altitude profiles of the 50th percentile (median) of each dataset.  
In case the second dataset is not provided it automatically reads the data from the Tromso EISCAT radar.  
X axis represents the variable studied and Y axis represents the altitude.
"""

import Data as D

import numpy as np
import plotly
import chart_studio.plotly as py 
import plotly.graph_objects as go
import plotly.figure_factory as ff
from plotly.subplots import make_subplots
import scipy.io
import scipy
import math

def plotAltProf_MedianComparison( VariableName, Buckets, CurveColor="dodgerblue", Buckets2=None, CurveColor2="dodgerblue", SuperTitle="" ):
    '''
    Creates comparison plots of two data sets.  
    The median values of each dataset are plotted together as altitude profiles for each MLT-Kp bin.
    In case the second dataset (Buckets2) is None, then the function reads data produced by the Tromso EISCAT radar.
    Args:
        VariableName (string): The physical variable on which the calculation has been applied. 
        Buckets (dictionary): The data structure which contains the statistical calculation results of the 1st dataset. See the function Data.init_ResultDataStructure() for details.
        CurveColor (string): The 1st dataset will be plotted with this color.
        Buckets2 (dictionary): The data structure which contains the statistical calculation results of the 2nd dataset. See the function Data.init_ResultDataStructure() for details.
        CurveColor2 (string): The 2nd dataset will be plotted with this color.
        SuperTitle (string): This title will be displayed at the top of the plot.
    '''
    HEIGHT_INTEGRATED_RATIO_ALL_average = 0
    HEIGHT_INTEGRATED_RATIO_UPPER_average = 0
    HEIGHT_INTEGRATED_RATIO_LOWER_average = 0
    TIEGCMarea_Upper = 0
    TIEGCMarea_Lower = 0
    TIEGCMarea2_Upper = 0
    TIEGCMarea2_Lower = 0
    EISCATcolor = CurveColor2 
    print("------------------ TIEGCM info start ------------------\n")
    print( "ALT_distance_of_a_bucket:", D.ALT_distance_of_a_bucket )
    print( "ALTsequence:", D.ALTsequence )
    print("------------------ TIEGCM info finish ------------------\n\n")
    
    if VariableName == "Joule Heating" or "JH" in VariableName:
        if Buckets2 != None:
            x_axes_range=[0, 3]
        else:
            x_axes_range=[0, 20]
        MultiplicationFactor = 10**8 
        new_units = "10^-8 W/m3"
    elif VariableName == "Pedersen Conductivity":
        x_axes_range=[0, 0.4]
        MultiplicationFactor = 10**3 
        new_units = "mS/m"
    else:
        x_axes_range=[0, 10]
        MultiplicationFactor = 1
        new_units = "?"

    
    # alter visibleALTsequence so that the point is displayed in the middle of the sub-bin
    visibleALTsequence = D.ALTsequence.copy()
    for i in range(1, len(visibleALTsequence)-1):
        visibleALTsequence[i] += D.ALT_distance_of_a_bucket/2
    visibleALTsequence[0] = D.ALTsequence[0]
    visibleALTsequence[-1] = D.ALTsequence[-1] + D.ALT_distance_of_a_bucket
    
    # construct the column MLT titles #("0-3", "3-6", "6-9", "9-12", "12-15", "15-18", "18-21", "21-24")
    ColumnTitles = list()
    
    for i in range(0, len(D.MLTsequence)):
        MLTfrom = int(D.MLTsequence[i])
        if MLTfrom > 24: MLTfrom -=24
        MLTto = int(D.MLTsequence[i]+D.MLT_duration_of_a_bucket)
        if MLTto > 24: MLTto -=24
        ColumnTitles.append( "MLT " + str(MLTfrom) + "-"  + str(MLTto) )
    # define secondary y-axis at the right of the plot
    mySpecs = list()
    for row in range(0, len(D.KPsequence)):
        mySpecs.append( list() )
        for col in range(0, len(D.MLTsequence)):
            mySpecs[row].append( {"secondary_y": True} )

    #make plot
    if VariableName == "Joule Heating": 
        XXtitle = 'Joule heating (10<sup>-8</sup> W/m<sup>3</sup>)'
    elif VariableName == "Pedersen Conductivity":
        XXtitle = 'Pedersen conductivity (mS/m)'
    else:
        XXtitle = VariableName
    fig = make_subplots(rows=len(D.KPsequence), cols=len(D.MLTsequence), x_title=XXtitle, shared_xaxes=True, shared_yaxes=True, vertical_spacing=0.035, horizontal_spacing=0.02, subplot_titles=ColumnTitles, specs=mySpecs)
    
    fig.update_layout( font=dict( family="arial black", size=24 ) )
    fig.update_annotations( font=dict( family="arial black", size=24) )
    #fig.update_xaxes(title_font_family="Arial black", title_font_size=20)
    #fig.update_yaxes(title_font_family="Arial black", title_font_size=20)
    fig.update_xaxes(tickfont_size=22)
    fig.update_yaxes(tickfont_size=22)
    fig.layout.annotations[4]["font"] = {'size': 30}  # this is the XXtitle at the bottom
    
    for aKP in D.KPsequence:
        for aMLT in D.MLTsequence:
            #Means = list()
            TIEGCMmedian = list()
            TIEGCMmedian2 = list()
            hits  = 0

            # compute TIEGCM 2ND RESULT percentiles
            if Buckets2 != None:
                TIEGCMarea_Upper = 0
                TIEGCMarea_Lower = 0
                TIEGCMarea2_Upper = 0
                TIEGCMarea2_Lower = 0
    
                for anALT in D.ALTsequence:
                    TIEGCMmedian2.append( Buckets2[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor )
                #
                for anALT in D.ALTsequence:
                    if anALT >= 120:
                        TIEGCMarea2_Upper += Buckets2[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                    else:
                        TIEGCMarea2_Lower += Buckets2[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                # plot TIEGCM 2ND RESULT median
                if CurveColor==CurveColor2:
                    linetype = 'dot'
                else:
                    linetype = 'solid'
                fig.add_trace( go.Scatter(x=TIEGCMmedian2, y=visibleALTsequence, mode='lines', fill=None, fillcolor=None, line=dict(color=CurveColor,width=4,dash=linetype,), showlegend=False), row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1 )
                # CALCULATE the Height_integration-Vaue for TIEGCM 2ND RESULT= area under median curve
                TIEGCMarea2 = 0
                for i in range(0, len(TIEGCMmedian2)):
                    if math.isnan(TIEGCMmedian2[i]) == False: 
                        TIEGCMarea2 += TIEGCMmedian2[i]*D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;

            # compute TIEGCM percentiles
            for anALT in D.ALTsequence:
                TIEGCMmedian.append( Buckets[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor )
            #
            for anALT in D.ALTsequence:
                if anALT >= 120:
                    TIEGCMarea_Upper += Buckets[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                else:
                    TIEGCMarea_Lower += Buckets[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
            # plot TIEGCM median
            fig.add_trace( go.Scatter(x=TIEGCMmedian, y=visibleALTsequence, mode='lines', fill=None, fillcolor=None, line=dict(color=CurveColor,width=4,), showlegend=False), row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1 )
            
            # CALCULATE the Height_integration-Vaue for TIEGCM = area under median curve
            TIEGCMarea = 0
            for i in range(0, len(TIEGCMmedian)):
                if math.isnan(TIEGCMmedian[i]) == False: 
                    if VariableName == "Joule Heating" or "JH" in VariableName: 
                        TIEGCMarea += TIEGCMmedian[i]*D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                    elif VariableName == "Pedersen Conductivity":
                        TIEGCMarea += TIEGCMmedian[i]*D.ALT_distance_of_a_bucket # area*1000 * math.pow(10,-8) * 1000;

            # read the median curve of EISCAT
            [EISCATmedian, EISCATmedianTHIN] = getEISCAT_MedianCurve(VariableName, aKP, aMLT)
            
            # CALCULATE the Height_integration-Vaue for EISCAT = area under median curve
            EISCATarea = 0.0
            for i in range(0, len(EISCATmedian)):
                if math.isnan(EISCATmedian[i]) == False: 
                    if VariableName == "Joule Heating" or "JH" in VariableName: 
                        EISCATarea += EISCATmedian[i]*0.01  #area += EISCATmedian[i]*1000 * math.pow(10,-8) * 1000;
                    elif VariableName == "Pedersen Conductivity":
                        EISCATarea += EISCATmedian[i]

            # Calculate the Percentage Difference
            try:
                SimilarityFactor_eiscat = (TIEGCMarea-EISCATarea) / ((TIEGCMarea+EISCATarea)/2)
                SimilarityFactor_eiscat = int(round(100*SimilarityFactor_eiscat,   0)) # %
            except:
                SimilarityFactor_eiscat = 0
            if Buckets2 != None:
                try:
                    SimilarityFactor_winds =  (TIEGCMarea2-TIEGCMarea) / TIEGCMarea 
                    SimilarityFactor_winds = int(round(100*SimilarityFactor_winds, 0)) # %    
                    #print ( "HEIGHT_INTEGRATED_RATIO_ALL", aMLT, aKP, "\t", round(TIEGCMarea/TIEGCMarea2,2) )
                    #print ( "HEIGHT_INTEGRATED_RATIO_UPPER", aMLT, aKP, "\t", round(TIEGCMarea_Upper/TIEGCMarea2_Upper ,2) )
                    #print ( "HEIGHT_INTEGRATED_RATIO_LOWER", aMLT, aKP, "\t", round(TIEGCMarea_Lower/TIEGCMarea2_Lower ,2) )
                    HEIGHT_INTEGRATED_RATIO_ALL_average += TIEGCMarea/TIEGCMarea2
                    HEIGHT_INTEGRATED_RATIO_UPPER_average += TIEGCMarea_Upper/TIEGCMarea2_Upper
                    HEIGHT_INTEGRATED_RATIO_LOWER_average += TIEGCMarea_Lower/TIEGCMarea2_Lower
                except:
                    pass
                #
                
            
                
            sim_factor_color = "purple"    

            # add annotations
            if VariableName=="Joule Heating":
                if Buckets2 != None:
                    fig.add_annotation(xref='x domain', yref='y domain', x=0.99, y=1, text=F"<b>{SimilarityFactor_winds}%</b>", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color=CurveColor) )
                    #fig.add_annotation(xref='x domain',yref='y domain', x=0.5, y=1, text=F"{round(TIEGCMarea_Upper/TIEGCMarea2_Upper,2)}", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color='black') )
                    #fig.add_annotation(xref='x domain',yref='y domain', x=0.5, y=0.5, text=F"{round(TIEGCMarea/TIEGCMarea2,2)}", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color='black') )
                    #fig.add_annotation(xref='x domain',yref='y domain', x=0.5, y=0, text=F"{round(TIEGCMarea_Lower/TIEGCMarea2_Lower,2)}", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color='black') )
                
                
            # add a trace in order to display secondary y-axis at the right
            fig.add_trace( go.Scatter(x=[-1000], y=[-1000], line=dict(color=CurveColor,width=1), showlegend=False), row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, secondary_y=True )
            
            
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    
    print ( "HEIGHT_INTEGRATED_RATIO_ALL", "average", "\t", round(HEIGHT_INTEGRATED_RATIO_ALL_average/12 ,2) )
    print ( "HEIGHT_INTEGRATED_RATIO_UPPER", "average", "\t", round(HEIGHT_INTEGRATED_RATIO_UPPER_average/12 ,2) )
    print ( "HEIGHT_INTEGRATED_RATIO_LOWER", "average", "\t", round(HEIGHT_INTEGRATED_RATIO_LOWER_average/12 ,2) )
    
    fig.update_xaxes( range=x_axes_range, row=1, col=1)
    fig.update_xaxes( range=x_axes_range, row=1, col=2)
    fig.update_xaxes( range=x_axes_range, row=1, col=3)
    fig.update_xaxes( range=x_axes_range, row=1, col=4)
    
    fig.update_xaxes( range=x_axes_range, row=2, col=1)
    fig.update_xaxes( range=x_axes_range, row=2, col=2)
    fig.update_xaxes( range=x_axes_range, row=2, col=3)
    fig.update_xaxes( range=x_axes_range, row=2, col=4)
    
    fig.update_xaxes( range=x_axes_range, row=3, col=1)
    fig.update_xaxes( range=x_axes_range, row=3, col=2)
    fig.update_xaxes( range=x_axes_range, row=3, col=3)
    fig.update_xaxes( range=x_axes_range, row=3, col=4)
    
    for aKP in D.KPsequence:
        fig.update_yaxes( title_text="Altitude(km)", row=D.KPsequence.index(aKP)+1, col=1, side='left', secondary_y=False)
        row_title = "Kp " + str(aKP) + " - "
        if aKP == 0:
            row_title +=  "2"
        elif aKP == 2:
            row_title +=  "4"
        else:
            row_title +=  "9"
        fig.update_yaxes( title_text=row_title, row=D.KPsequence.index(aKP)+1, col=len(D.MLTsequence),  side='right', secondary_y=True, showticklabels=False )
        for aMLT in D.MLTsequence:
            fig.update_yaxes( row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, secondary_y=True, showticklabels=False )
    #fig.update_xaxes( range=x_axes_range )
    fig.update_yaxes( range=[80, 150], tick0=90, dtick=20 )  
    fig.update_layout( title = SuperTitle,
                       width=400+len(D.MLTsequence)*250, height=200+200*len(D.KPsequence), showlegend=False, legend_orientation="h", legend_y=-0.04) 

    if Buckets2 == None:
        plotEISCAT( VariableName, fig )
    
    plotly.offline.init_notebook_mode(connected=True)
    plotly.offline.iplot(fig) 

    # plot more zoom versions
    '''
    new_x_axes_range = [x * (2/3) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (1/2) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (3/2) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (2.5) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (10) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    '''

    
    
    

def getEISCAT_MedianCurve( VariableName, aKP, aMLT ):
    '''
    Args:
        VariableName (string): The physical variable on which the calculation has been applied. 
        aKP (float): A Kp index value (0-9)
        aMLT (float): A Magnetic Local Time value 
    Returns:
        A list of points representing an altitude profile median curve for the desired KP and MLT combination
    '''
    if aMLT > 24: aMLT -= 24
    
    Values = None
    matlabStruct = scipy.io.loadmat('./EISCAT_DATA/data_2009_2019_TS.mat')
    
    allALTs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][0] ).flatten()
    allKPs  = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][1][0] ) )
    allMLTs = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][2][0] )[:-1] )
    allJHs  = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][3] )
    allPEDs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][4] )
    
    if VariableName == "Pedersen Conductivity":
        Values = allPEDs
        MultiplicationFactor = 10**3 
        new_units = "mS/m"
    else:
        Values = allJHs
        MultiplicationFactor = 10**8 
        new_units = "10^-8 W/m3"
    
    ALTsequence =  allALTs
    MLTsequence = allMLTs
    KPsequence = [ 0, 2, 4 ]
    MLT_duration_of_a_profile = 6
    
    # alter visibleALTsequence so that the point is displayed in the middle of the sub-bin
    visibleALTsequence = ALTsequence.copy()
    for i in range(1, len(visibleALTsequence)-1):
        visibleALTsequence[i] += 0.5
    
    MedianCurve = Values[KPsequence.index(aKP), MLTsequence.index(aMLT), :, 2] * MultiplicationFactor
    
    #print( "~~~~~~~~~~~~Thinning EISCAT median to compare with TIEGCM median", len(ALTsequence) )
    EISCATmedianTHIN = []
    for i in range( 0, len(MedianCurve) ):
        if ALTsequence[i] in D.ALTsequence:
            EISCATmedianTHIN.append( MedianCurve[i] )
    
    return [ MedianCurve, EISCATmedianTHIN ]

    

def plotEISCAT( VariableName, fig ):
    '''
    Adds altitude profile curves of the median value of a variable as calculated by EISCAT
    Args:
        VariableName (string): The physical variable on which the calculation has been applied. 
        fig (plotly object): the plotly figure upon which the EISCAT altitude profiles of the median value will be plotted.
    '''
    EISCATcolor = "limegreen"
    
    matlabStruct = scipy.io.loadmat('./EISCAT_DATA/data_2009_2019_TS.mat')
    
    allALTs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][0] ).flatten()
    allKPs  = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][1][0] ) )
    allMLTs = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][2][0] )[:-1] )
    allJHs  = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][3] )
    allPEDs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][4] )
    
    print("------------------ EISCAT info start ------------------")
    print( "Altitudes:", allALTs[0], allALTs[1], "...", allALTs[-1] )
    print( "KPs:", allKPs  )
    print( "MLTs:", allMLTs )
    print( "JHs shape:", allJHs.shape )
    print( "PEDs shape:", allPEDs.shape )    
    print("------------------ EISCAT info finish ------------------\n\n")
    
    if VariableName == "Pedersen Conductivity":
        Values = allPEDs
        x_axes_range=[0, 0.4]
        MultiplicationFactor = 10**3 
        new_units = "mS/m"
    else:
        Values = allJHs
        x_axes_range=[0, 20]
        MultiplicationFactor = 10**8 
        new_units = "10^-8 W/m3"
    
    ALTsequence =  allALTs
    MLTsequence = allMLTs
    KPsequence = [ 0, 2, 4 ]  #list( mat_medians[ 'jouleMedians' ][0][0][3] )
    MLT_duration_of_a_profile = 6
    
    # alter visibleALTsequence so that the point is displayed in the middle of the sub-bin
    visibleALTsequence = ALTsequence.copy()
    for i in range(1, len(visibleALTsequence)-1):
        visibleALTsequence[i] += 0.5
    
    for aKP in KPsequence:
        for aMLT in MLTsequence:
            #Means = list()
            EISCATmedian = list()            
            hits  = 0
            
            # compute percentiles
            EISCATmedian = Values[KPsequence.index(aKP), MLTsequence.index(aMLT), :, 2] * MultiplicationFactor #EISCATmedian = JHmedians[1,1,:] * MultiplicationFactor 

            fig.add_trace( go.Scatter(x=EISCATmedian, y=visibleALTsequence, mode='lines', fill=None, fillcolor=EISCATcolor, line=dict(color=EISCATcolor,width=4,), showlegend=False), row=KPsequence.index(aKP)+1, col=MLTsequence.index(aMLT)+1 )
            
      

Functions

def getEISCAT_MedianCurve(VariableName, aKP, aMLT)

Args

VariableName : string
The physical variable on which the calculation has been applied.
aKP : float
A Kp index value (0-9)
aMLT : float
A Magnetic Local Time value

Returns

A list of points representing an altitude profile median curve for the desired KP and MLT combination

Expand source code
def getEISCAT_MedianCurve( VariableName, aKP, aMLT ):
    '''
    Args:
        VariableName (string): The physical variable on which the calculation has been applied. 
        aKP (float): A Kp index value (0-9)
        aMLT (float): A Magnetic Local Time value 
    Returns:
        A list of points representing an altitude profile median curve for the desired KP and MLT combination
    '''
    if aMLT > 24: aMLT -= 24
    
    Values = None
    matlabStruct = scipy.io.loadmat('./EISCAT_DATA/data_2009_2019_TS.mat')
    
    allALTs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][0] ).flatten()
    allKPs  = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][1][0] ) )
    allMLTs = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][2][0] )[:-1] )
    allJHs  = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][3] )
    allPEDs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][4] )
    
    if VariableName == "Pedersen Conductivity":
        Values = allPEDs
        MultiplicationFactor = 10**3 
        new_units = "mS/m"
    else:
        Values = allJHs
        MultiplicationFactor = 10**8 
        new_units = "10^-8 W/m3"
    
    ALTsequence =  allALTs
    MLTsequence = allMLTs
    KPsequence = [ 0, 2, 4 ]
    MLT_duration_of_a_profile = 6
    
    # alter visibleALTsequence so that the point is displayed in the middle of the sub-bin
    visibleALTsequence = ALTsequence.copy()
    for i in range(1, len(visibleALTsequence)-1):
        visibleALTsequence[i] += 0.5
    
    MedianCurve = Values[KPsequence.index(aKP), MLTsequence.index(aMLT), :, 2] * MultiplicationFactor
    
    #print( "~~~~~~~~~~~~Thinning EISCAT median to compare with TIEGCM median", len(ALTsequence) )
    EISCATmedianTHIN = []
    for i in range( 0, len(MedianCurve) ):
        if ALTsequence[i] in D.ALTsequence:
            EISCATmedianTHIN.append( MedianCurve[i] )
    
    return [ MedianCurve, EISCATmedianTHIN ]
def plotAltProf_MedianComparison(VariableName, Buckets, CurveColor='dodgerblue', Buckets2=None, CurveColor2='dodgerblue', SuperTitle='')

Creates comparison plots of two data sets.
The median values of each dataset are plotted together as altitude profiles for each MLT-Kp bin. In case the second dataset (Buckets2) is None, then the function reads data produced by the Tromso EISCAT radar.

Args

VariableName : string
The physical variable on which the calculation has been applied.
Buckets : dictionary
The data structure which contains the statistical calculation results of the 1st dataset. See the function Data.init_ResultDataStructure() for details.
CurveColor : string
The 1st dataset will be plotted with this color.
Buckets2 : dictionary
The data structure which contains the statistical calculation results of the 2nd dataset. See the function Data.init_ResultDataStructure() for details.
CurveColor2 : string
The 2nd dataset will be plotted with this color.
SuperTitle : string
This title will be displayed at the top of the plot.
Expand source code
def plotAltProf_MedianComparison( VariableName, Buckets, CurveColor="dodgerblue", Buckets2=None, CurveColor2="dodgerblue", SuperTitle="" ):
    '''
    Creates comparison plots of two data sets.  
    The median values of each dataset are plotted together as altitude profiles for each MLT-Kp bin.
    In case the second dataset (Buckets2) is None, then the function reads data produced by the Tromso EISCAT radar.
    Args:
        VariableName (string): The physical variable on which the calculation has been applied. 
        Buckets (dictionary): The data structure which contains the statistical calculation results of the 1st dataset. See the function Data.init_ResultDataStructure() for details.
        CurveColor (string): The 1st dataset will be plotted with this color.
        Buckets2 (dictionary): The data structure which contains the statistical calculation results of the 2nd dataset. See the function Data.init_ResultDataStructure() for details.
        CurveColor2 (string): The 2nd dataset will be plotted with this color.
        SuperTitle (string): This title will be displayed at the top of the plot.
    '''
    HEIGHT_INTEGRATED_RATIO_ALL_average = 0
    HEIGHT_INTEGRATED_RATIO_UPPER_average = 0
    HEIGHT_INTEGRATED_RATIO_LOWER_average = 0
    TIEGCMarea_Upper = 0
    TIEGCMarea_Lower = 0
    TIEGCMarea2_Upper = 0
    TIEGCMarea2_Lower = 0
    EISCATcolor = CurveColor2 
    print("------------------ TIEGCM info start ------------------\n")
    print( "ALT_distance_of_a_bucket:", D.ALT_distance_of_a_bucket )
    print( "ALTsequence:", D.ALTsequence )
    print("------------------ TIEGCM info finish ------------------\n\n")
    
    if VariableName == "Joule Heating" or "JH" in VariableName:
        if Buckets2 != None:
            x_axes_range=[0, 3]
        else:
            x_axes_range=[0, 20]
        MultiplicationFactor = 10**8 
        new_units = "10^-8 W/m3"
    elif VariableName == "Pedersen Conductivity":
        x_axes_range=[0, 0.4]
        MultiplicationFactor = 10**3 
        new_units = "mS/m"
    else:
        x_axes_range=[0, 10]
        MultiplicationFactor = 1
        new_units = "?"

    
    # alter visibleALTsequence so that the point is displayed in the middle of the sub-bin
    visibleALTsequence = D.ALTsequence.copy()
    for i in range(1, len(visibleALTsequence)-1):
        visibleALTsequence[i] += D.ALT_distance_of_a_bucket/2
    visibleALTsequence[0] = D.ALTsequence[0]
    visibleALTsequence[-1] = D.ALTsequence[-1] + D.ALT_distance_of_a_bucket
    
    # construct the column MLT titles #("0-3", "3-6", "6-9", "9-12", "12-15", "15-18", "18-21", "21-24")
    ColumnTitles = list()
    
    for i in range(0, len(D.MLTsequence)):
        MLTfrom = int(D.MLTsequence[i])
        if MLTfrom > 24: MLTfrom -=24
        MLTto = int(D.MLTsequence[i]+D.MLT_duration_of_a_bucket)
        if MLTto > 24: MLTto -=24
        ColumnTitles.append( "MLT " + str(MLTfrom) + "-"  + str(MLTto) )
    # define secondary y-axis at the right of the plot
    mySpecs = list()
    for row in range(0, len(D.KPsequence)):
        mySpecs.append( list() )
        for col in range(0, len(D.MLTsequence)):
            mySpecs[row].append( {"secondary_y": True} )

    #make plot
    if VariableName == "Joule Heating": 
        XXtitle = 'Joule heating (10<sup>-8</sup> W/m<sup>3</sup>)'
    elif VariableName == "Pedersen Conductivity":
        XXtitle = 'Pedersen conductivity (mS/m)'
    else:
        XXtitle = VariableName
    fig = make_subplots(rows=len(D.KPsequence), cols=len(D.MLTsequence), x_title=XXtitle, shared_xaxes=True, shared_yaxes=True, vertical_spacing=0.035, horizontal_spacing=0.02, subplot_titles=ColumnTitles, specs=mySpecs)
    
    fig.update_layout( font=dict( family="arial black", size=24 ) )
    fig.update_annotations( font=dict( family="arial black", size=24) )
    #fig.update_xaxes(title_font_family="Arial black", title_font_size=20)
    #fig.update_yaxes(title_font_family="Arial black", title_font_size=20)
    fig.update_xaxes(tickfont_size=22)
    fig.update_yaxes(tickfont_size=22)
    fig.layout.annotations[4]["font"] = {'size': 30}  # this is the XXtitle at the bottom
    
    for aKP in D.KPsequence:
        for aMLT in D.MLTsequence:
            #Means = list()
            TIEGCMmedian = list()
            TIEGCMmedian2 = list()
            hits  = 0

            # compute TIEGCM 2ND RESULT percentiles
            if Buckets2 != None:
                TIEGCMarea_Upper = 0
                TIEGCMarea_Lower = 0
                TIEGCMarea2_Upper = 0
                TIEGCMarea2_Lower = 0
    
                for anALT in D.ALTsequence:
                    TIEGCMmedian2.append( Buckets2[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor )
                #
                for anALT in D.ALTsequence:
                    if anALT >= 120:
                        TIEGCMarea2_Upper += Buckets2[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                    else:
                        TIEGCMarea2_Lower += Buckets2[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                # plot TIEGCM 2ND RESULT median
                if CurveColor==CurveColor2:
                    linetype = 'dot'
                else:
                    linetype = 'solid'
                fig.add_trace( go.Scatter(x=TIEGCMmedian2, y=visibleALTsequence, mode='lines', fill=None, fillcolor=None, line=dict(color=CurveColor,width=4,dash=linetype,), showlegend=False), row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1 )
                # CALCULATE the Height_integration-Vaue for TIEGCM 2ND RESULT= area under median curve
                TIEGCMarea2 = 0
                for i in range(0, len(TIEGCMmedian2)):
                    if math.isnan(TIEGCMmedian2[i]) == False: 
                        TIEGCMarea2 += TIEGCMmedian2[i]*D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;

            # compute TIEGCM percentiles
            for anALT in D.ALTsequence:
                TIEGCMmedian.append( Buckets[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor )
            #
            for anALT in D.ALTsequence:
                if anALT >= 120:
                    TIEGCMarea_Upper += Buckets[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                else:
                    TIEGCMarea_Lower += Buckets[aKP, anALT, D.LAT_min, aMLT, "Percentile50"] * MultiplicationFactor * D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
            # plot TIEGCM median
            fig.add_trace( go.Scatter(x=TIEGCMmedian, y=visibleALTsequence, mode='lines', fill=None, fillcolor=None, line=dict(color=CurveColor,width=4,), showlegend=False), row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1 )
            
            # CALCULATE the Height_integration-Vaue for TIEGCM = area under median curve
            TIEGCMarea = 0
            for i in range(0, len(TIEGCMmedian)):
                if math.isnan(TIEGCMmedian[i]) == False: 
                    if VariableName == "Joule Heating" or "JH" in VariableName: 
                        TIEGCMarea += TIEGCMmedian[i]*D.ALT_distance_of_a_bucket * 0.01 # area*1000 * math.pow(10,-8) * 1000;
                    elif VariableName == "Pedersen Conductivity":
                        TIEGCMarea += TIEGCMmedian[i]*D.ALT_distance_of_a_bucket # area*1000 * math.pow(10,-8) * 1000;

            # read the median curve of EISCAT
            [EISCATmedian, EISCATmedianTHIN] = getEISCAT_MedianCurve(VariableName, aKP, aMLT)
            
            # CALCULATE the Height_integration-Vaue for EISCAT = area under median curve
            EISCATarea = 0.0
            for i in range(0, len(EISCATmedian)):
                if math.isnan(EISCATmedian[i]) == False: 
                    if VariableName == "Joule Heating" or "JH" in VariableName: 
                        EISCATarea += EISCATmedian[i]*0.01  #area += EISCATmedian[i]*1000 * math.pow(10,-8) * 1000;
                    elif VariableName == "Pedersen Conductivity":
                        EISCATarea += EISCATmedian[i]

            # Calculate the Percentage Difference
            try:
                SimilarityFactor_eiscat = (TIEGCMarea-EISCATarea) / ((TIEGCMarea+EISCATarea)/2)
                SimilarityFactor_eiscat = int(round(100*SimilarityFactor_eiscat,   0)) # %
            except:
                SimilarityFactor_eiscat = 0
            if Buckets2 != None:
                try:
                    SimilarityFactor_winds =  (TIEGCMarea2-TIEGCMarea) / TIEGCMarea 
                    SimilarityFactor_winds = int(round(100*SimilarityFactor_winds, 0)) # %    
                    #print ( "HEIGHT_INTEGRATED_RATIO_ALL", aMLT, aKP, "\t", round(TIEGCMarea/TIEGCMarea2,2) )
                    #print ( "HEIGHT_INTEGRATED_RATIO_UPPER", aMLT, aKP, "\t", round(TIEGCMarea_Upper/TIEGCMarea2_Upper ,2) )
                    #print ( "HEIGHT_INTEGRATED_RATIO_LOWER", aMLT, aKP, "\t", round(TIEGCMarea_Lower/TIEGCMarea2_Lower ,2) )
                    HEIGHT_INTEGRATED_RATIO_ALL_average += TIEGCMarea/TIEGCMarea2
                    HEIGHT_INTEGRATED_RATIO_UPPER_average += TIEGCMarea_Upper/TIEGCMarea2_Upper
                    HEIGHT_INTEGRATED_RATIO_LOWER_average += TIEGCMarea_Lower/TIEGCMarea2_Lower
                except:
                    pass
                #
                
            
                
            sim_factor_color = "purple"    

            # add annotations
            if VariableName=="Joule Heating":
                if Buckets2 != None:
                    fig.add_annotation(xref='x domain', yref='y domain', x=0.99, y=1, text=F"<b>{SimilarityFactor_winds}%</b>", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color=CurveColor) )
                    #fig.add_annotation(xref='x domain',yref='y domain', x=0.5, y=1, text=F"{round(TIEGCMarea_Upper/TIEGCMarea2_Upper,2)}", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color='black') )
                    #fig.add_annotation(xref='x domain',yref='y domain', x=0.5, y=0.5, text=F"{round(TIEGCMarea/TIEGCMarea2,2)}", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color='black') )
                    #fig.add_annotation(xref='x domain',yref='y domain', x=0.5, y=0, text=F"{round(TIEGCMarea_Lower/TIEGCMarea2_Lower,2)}", showarrow=False, row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, font=dict(color='black') )
                
                
            # add a trace in order to display secondary y-axis at the right
            fig.add_trace( go.Scatter(x=[-1000], y=[-1000], line=dict(color=CurveColor,width=1), showlegend=False), row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, secondary_y=True )
            
            
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    
    print ( "HEIGHT_INTEGRATED_RATIO_ALL", "average", "\t", round(HEIGHT_INTEGRATED_RATIO_ALL_average/12 ,2) )
    print ( "HEIGHT_INTEGRATED_RATIO_UPPER", "average", "\t", round(HEIGHT_INTEGRATED_RATIO_UPPER_average/12 ,2) )
    print ( "HEIGHT_INTEGRATED_RATIO_LOWER", "average", "\t", round(HEIGHT_INTEGRATED_RATIO_LOWER_average/12 ,2) )
    
    fig.update_xaxes( range=x_axes_range, row=1, col=1)
    fig.update_xaxes( range=x_axes_range, row=1, col=2)
    fig.update_xaxes( range=x_axes_range, row=1, col=3)
    fig.update_xaxes( range=x_axes_range, row=1, col=4)
    
    fig.update_xaxes( range=x_axes_range, row=2, col=1)
    fig.update_xaxes( range=x_axes_range, row=2, col=2)
    fig.update_xaxes( range=x_axes_range, row=2, col=3)
    fig.update_xaxes( range=x_axes_range, row=2, col=4)
    
    fig.update_xaxes( range=x_axes_range, row=3, col=1)
    fig.update_xaxes( range=x_axes_range, row=3, col=2)
    fig.update_xaxes( range=x_axes_range, row=3, col=3)
    fig.update_xaxes( range=x_axes_range, row=3, col=4)
    
    for aKP in D.KPsequence:
        fig.update_yaxes( title_text="Altitude(km)", row=D.KPsequence.index(aKP)+1, col=1, side='left', secondary_y=False)
        row_title = "Kp " + str(aKP) + " - "
        if aKP == 0:
            row_title +=  "2"
        elif aKP == 2:
            row_title +=  "4"
        else:
            row_title +=  "9"
        fig.update_yaxes( title_text=row_title, row=D.KPsequence.index(aKP)+1, col=len(D.MLTsequence),  side='right', secondary_y=True, showticklabels=False )
        for aMLT in D.MLTsequence:
            fig.update_yaxes( row=D.KPsequence.index(aKP)+1, col=D.MLTsequence.index(aMLT)+1, secondary_y=True, showticklabels=False )
    #fig.update_xaxes( range=x_axes_range )
    fig.update_yaxes( range=[80, 150], tick0=90, dtick=20 )  
    fig.update_layout( title = SuperTitle,
                       width=400+len(D.MLTsequence)*250, height=200+200*len(D.KPsequence), showlegend=False, legend_orientation="h", legend_y=-0.04) 

    if Buckets2 == None:
        plotEISCAT( VariableName, fig )
    
    plotly.offline.init_notebook_mode(connected=True)
    plotly.offline.iplot(fig) 

    # plot more zoom versions
    '''
    new_x_axes_range = [x * (2/3) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (1/2) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (3/2) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (2.5) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    new_x_axes_range = [x * (10) for x in x_axes_range]
    fig.update_xaxes( range=new_x_axes_range )
    plotly.offline.iplot(fig) 
    '''
def plotEISCAT(VariableName, fig)

Adds altitude profile curves of the median value of a variable as calculated by EISCAT

Args

VariableName : string
The physical variable on which the calculation has been applied.
fig : plotly object
the plotly figure upon which the EISCAT altitude profiles of the median value will be plotted.
Expand source code
def plotEISCAT( VariableName, fig ):
    '''
    Adds altitude profile curves of the median value of a variable as calculated by EISCAT
    Args:
        VariableName (string): The physical variable on which the calculation has been applied. 
        fig (plotly object): the plotly figure upon which the EISCAT altitude profiles of the median value will be plotted.
    '''
    EISCATcolor = "limegreen"
    
    matlabStruct = scipy.io.loadmat('./EISCAT_DATA/data_2009_2019_TS.mat')
    
    allALTs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][0] ).flatten()
    allKPs  = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][1][0] ) )
    allMLTs = list( np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][2][0] )[:-1] )
    allJHs  = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][3] )
    allPEDs = np.array( matlabStruct[ 'data_2009_2019_TS' ][0][0][4] )
    
    print("------------------ EISCAT info start ------------------")
    print( "Altitudes:", allALTs[0], allALTs[1], "...", allALTs[-1] )
    print( "KPs:", allKPs  )
    print( "MLTs:", allMLTs )
    print( "JHs shape:", allJHs.shape )
    print( "PEDs shape:", allPEDs.shape )    
    print("------------------ EISCAT info finish ------------------\n\n")
    
    if VariableName == "Pedersen Conductivity":
        Values = allPEDs
        x_axes_range=[0, 0.4]
        MultiplicationFactor = 10**3 
        new_units = "mS/m"
    else:
        Values = allJHs
        x_axes_range=[0, 20]
        MultiplicationFactor = 10**8 
        new_units = "10^-8 W/m3"
    
    ALTsequence =  allALTs
    MLTsequence = allMLTs
    KPsequence = [ 0, 2, 4 ]  #list( mat_medians[ 'jouleMedians' ][0][0][3] )
    MLT_duration_of_a_profile = 6
    
    # alter visibleALTsequence so that the point is displayed in the middle of the sub-bin
    visibleALTsequence = ALTsequence.copy()
    for i in range(1, len(visibleALTsequence)-1):
        visibleALTsequence[i] += 0.5
    
    for aKP in KPsequence:
        for aMLT in MLTsequence:
            #Means = list()
            EISCATmedian = list()            
            hits  = 0
            
            # compute percentiles
            EISCATmedian = Values[KPsequence.index(aKP), MLTsequence.index(aMLT), :, 2] * MultiplicationFactor #EISCATmedian = JHmedians[1,1,:] * MultiplicationFactor 

            fig.add_trace( go.Scatter(x=EISCATmedian, y=visibleALTsequence, mode='lines', fill=None, fillcolor=EISCATcolor, line=dict(color=EISCATcolor,width=4,), showlegend=False), row=KPsequence.index(aKP)+1, col=MLTsequence.index(aMLT)+1 )