diff --git a/scripts/CMakeLists.txt b/scripts/CMakeLists.txt index 5f7cad3..414059b 100644 --- a/scripts/CMakeLists.txt +++ b/scripts/CMakeLists.txt @@ -13,4 +13,5 @@ install_scripts (cv_header_list.sh) install_scripts (varbcdiagnose.sh plotvarbccoeff.py) install_scripts (rundiacov.sh diag_diacov.sh plotdiacov.py) install_scripts (dfscomp.sh plotdfs.py) +install_scripts (dep_stats.py) #install_scripts (diacov.sh festat.sh) diff --git a/scripts/JOdiag3DVar.py b/scripts/JOdiag3DVar.py new file mode 100755 index 0000000..9fbf79a --- /dev/null +++ b/scripts/JOdiag3DVar.py @@ -0,0 +1,305 @@ +### +import datetime +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +from matplotlib.dates import DayLocator, HourLocator +from matplotlib.ticker import (MultipleLocator, + FormatStrFormatter, + AutoMinorLocator) +import numpy as np +import os +import pandas as pd +import re +######################################### +years = ['2024'] +months = ['03'] +days = ['06'] +cycles = ['00','03','06','09','12','15','18', '21'] +path = '/scratch/nkim/logFilesFromDINI' +out_path = '/scratch/nkim/logFilesFromDINI/JOstatistics' +keep_JOtables = True +figs = True +Variable = 'Z' +Obstype = 'SYNOP' +################################# +def main (): + + print (">>>> Diagnostics using cost function information from log files <<<<<") + + Codetype = SetCodetype ( Obstype, Variable) + +# Sanity checks + if (np.size(Codetype) == 0 or np.size(Variable) == 0): + print ("ERROR: check your settings for Codetype and Variable") + print("Codetype: ", Codetype) + print("Variable : ", Variable) + exit() + + + num_code = 0 + + for cty in Codetype: + num_code = num_code + 1 + JoGlobalAll = [] + JoCtypeAll = [] + date_dfAll = [] + + for year in years: + for month in months : + + for day in days: + for cy in cycles: + fileLog = '{}/HM_Date_{}{}{}{}.html'.format(path,year,month,day,cy) + fileJOtable = '{}/JOtable_{}{}{}{}.txt'.format(out_path,year,month,day,cy) + + print(">>>>>>> Reading the log file in: ", fileLog) + extractJOtable(fileLog, fileJOtable) + + print(">>>>>>> Reading the JoGlobal") + Global = readGlobal(fileJOtable) + df1 = pd.DataFrame(Global['JoGlobal']) + df2 = pd.DataFrame(Global['Ntotal']) + dataListG = createDataList(df1, df2) + JoGlobal = pd.DataFrame(dataListG) + + d = datetime.datetime(int(year),int(month), + int(day),int(cy)).strftime("%Y-%m-%d %H") + date_df = pd.DataFrame ({"datetime": [d]} ) + JoGlobalAll.append(JoGlobal) + date_dfAll.append(date_df) + + print(">>>>>>> Reading the Jo{}: ".format(cty)) + Observations = readObservations(fileJOtable, cty, Variable) + df3 = pd.DataFrame(Observations['JoOBS']) + df4 = pd.DataFrame(Observations['NOBS']) + dataListC =createDataList(df3, df4) + JoCTy = pd.DataFrame(dataListC) + JoCtypeAll.append(JoCTy) + + if keep_JOtables is False: + print(">>>>>>> deleting: ", fileJOtable) + os.remove (fileJOtable) + + print(">>>>>>> Creating the list of Jo(s)") + JoGlobalAll = pd.concat(JoGlobalAll, ignore_index=True, axis=0) + JoCtypeAll = pd.concat(JoCtypeAll, ignore_index=True, axis=0) + dateAll = pd.concat(date_dfAll, ignore_index= True, axis=0) + + + JoGlobalAll = pd.concat([dateAll, JoGlobalAll], axis=1) + JoCtypeAll = pd.concat([dateAll, JoCtypeAll], axis=1) + + print(">>>>>>> writing Jo(s) in :", out_path) + Jodf = writeJoGlobal(out_path, JoGlobalAll,"Global") + JoCTydf = writeJoGlobal(out_path, JoCtypeAll, "Observations_{}_{}".format(Obstype,num_code, Variable)) + print(">>>>>>> saving Jo(s) plots in :", out_path) + + if (figs): plotJo('JoCtype', out_path, JoCTydf, Obstype, num_code, cty, Variable) + + if (figs): plotJo('JoGlobal', out_path, Jodf,"Global","", "", "") + +####### +def SetCodetype ( Obstype, Variable): + + if (Obstype == 'SYNOP') : + Codetype = ['Codetype 14', 'Codetype 24', \ + 'Codetype 170', 'Codetype 182'] + if ( Variable != 'Z' ): + print ("{} variable not supported for SYNOP (only Z)".format(Variable)) + exit() + elif (Obstype == 'Aircraft'): + Codetype = ['Codetype 141', 'Codetype 146', 'Codetype 147'] + if ( Variable != 'U' and Variable != 'T'): + print ("{} variable not supported for Aircraft (only U and T)".format(Variable)) + exit() + + elif (Obstype == 'Radiosondes'): + Codetype = ['Codetype 109 === BUFR Land TEMP' ] + if ( Variable != 'U' and Variable != 'T' and Variable != 'Q'): + print ("{} variable not supported for Radiosondes (only U, T and Q)".format(Variable)) + exit() + elif (Obstype == 'AMSUA'): + Codetype = ['Codetype 210 === metop 1 3 SENSOR=amsua', \ + 'Codetype 210 === metop 3 5 SENSOR=amsua', \ + 'Codetype 210 === noaa 18 209 SENSOR=amsua', \ + 'Codetype 210 === noaa 19 223 SENSOR=amsua'] + if ( Variable != 'RAD'): + print ("{} variable not supported for AMSUA (only RAD)".format(Variable)) + exit() + elif (Obstype == 'MHS'): + Codetype = ['Codetype 210 === metop 1 3 SENSOR=mhs', \ + 'Codetype 210 === metop 3 5 SENSOR=mhs', \ + 'Codetype 210 === noaa 19 223 SENSOR=mhs'] + if ( Variable != 'RAD'): + print ("{} variable not supported for MHS (only RAD)".format(Variable)) + exit() + elif (Obstype == 'ATMS'): + Codetype = ['Codetype 210 === jpss 0 224 SENSOR=atms', \ + 'Codetype 210 === noaa 20 225 SENSOR=atms'] + if ( Variable != 'RAD'): + print ("{} variable not supported for ATMS (only RAD)".format(Variable)) + exit() + elif (Obstype == 'MWHS2'): + Codetype = ['Codetype 210 === fy3 4'] + if ( Variable != 'RAD'): + print ("{} variable not supported for MWHS2 (only RAD)".format(Variable)) + exit() + elif (Obstype == 'IASI' ): + Codetype = ['Codetype 210 === metop 1 5 SENSOR=iasi',\ + 'Codetype 210 === metop 3 5 SENSOR=iasi'] + if ( Variable != 'RAD'): + print ("{} variable not supported for IASI (only RAD)".format(Variable)) + exit() + elif (Obstype == 'ASCAT'): + Codetype = ['Codetype 139 === METOP-B ASCAT', 'Codetype 139 === METOP-C ASCAT'] + if ( Variable != 'U10'): + print ("{} variable not supported for ASCAT (only U10)".format(Variable)) + exit() + elif (Obstype == 'LIMB'): + Codetype = ['Codetype 250 === METOP-1', 'Codetype 250 === METOP-3',\ + 'Codetype 250 === SPIRE'] + if ( Variable != 'RO'): + print ("{} variable not supported for LIMB (only RO)".format(Variable)) + exit() + elif (Obstype == 'AMV'): + Codetype = ['Codetype 90 === METOP 5 METHOD=IR',\ + 'Codetype 90 === METEOSAT 56 METHOD=WVCL1',\ + 'Codetype 90 === METEOSAT 56 METHOD=WVCL2', \ + 'Codetype 90 === METEOSAT 56 METHOD=IR3',\ + 'Codetype 90 === METEOSAT 57 METHOD=WVCL1',\ + 'Codetype 90 === METEOSAT 57 METHOD=WVCL2',\ + 'Codetype 90 === METEOSAT 57 METHOD=IR3',\ + 'Codetype 90 === dual-MTOP 852 METHOD=IR'] + if ( Variable != 'U' and Variable != 'T'): + print ("{} variable not supported for AMV (only U and T)".format(Variable)) + exit() + else: + print ("{} Obstype not supported.".format(Obstype)) + + + return Codetype +#### +def extractJOtable(fileLog, fileJOtable): + + + in_JOtable = False + mylines = [] + fh = open(fileJOtable, 'w') + with open (fileLog, 'rt') as myfile: + for myline in myfile: + # add its contents to mylines. + if ("Diagnostic JO-table (JOT)" and "MINIMISATION JOB" and "NSIM4D= 999 " in myline): + in_JOtable = True + elif("End of JO-table (JOT)" in myline) : + in_JOtable = False + if in_JOtable: + fh.write(myline) + fh.close() +#### +def readObservations (fileJOtable, code, var): + + in_JOtable = True + mylines = [] + rec1 = [] + rec2 = [] + nextLine = False + # + var = var + " " + + with open (fileJOtable, 'rt') as myfile: + # searches for codetypes + for mylines in myfile: + if (in_JOtable): + if ( re.search(code, mylines) ): + nextLine = True + if (nextLine): + if ( re.search (var , mylines)): + record = mylines.split() + rec1.append(int(record[1])) + rec2.append(float(record[3])) + nextLine = False + + if ("Jo Global :" in mylines ): + in_JOtable = False + + in_JOtable = True + + datalist = ({ + 'JoOBS' : [rec2], + 'NOBS' : [rec1] + + }) + + print(datalist) + return datalist + + +#### +def readGlobal(fileJOtable): + + mylines = [] + rec1 = [] + rec2 = [] + + with open(fileJOtable,'rt') as myfile: + for mylines in myfile: + if re.search("Jo Global", mylines): + record = mylines.split() + rec1.append(int(record[3])) + rec2.append(float(record[5])) + + datalist = ({ + 'JoGlobal' : [rec2], + 'Ntotal' : [rec1] + + }) + return datalist +####### +def createDataList(data1, data2): + + # sanity test + if (np.size(data2) == 0) : + data2[0] = 0 + data1[0] = 0 + + datalist = ({ + 'endMinim': [data1[0]], + 'N': [data2[0]], + }) + + return datalist +##### +def writeJoGlobal(out_path,JoGlobal,name): + + filename = '{}/Jo{}.csv'.format(out_path,name) + df = pd.DataFrame(JoGlobal) + df.to_csv(filename, index=False, sep =' ') + print (">>>>> df: ", df) + return df + +def plotJo(label,path, Jo, obs, num, cty, var): + + Jo = Jo.set_index("datetime") + print ("========= JO ========") + png = "{}/Jo{}{}{}.png".format(path,obs,num,var) + plt.close("all") + fig, ax = plt.subplots(figsize=(12, 6)) + ax.plot(Jo["endMinim"].astype(float),label='end Minim',linestyle='solid',color='blue') + ax.set_ylim(0,1.5) + ax.set_ylabel("Jo/n") + ax2 = ax.twinx() + ax2.plot(Jo["N"].astype(int),label='N',linestyle='dashed',color='grey') + ax2.set_ylabel ("N total") + legend = ax.legend(loc='upper right', fontsize='x-small') + plt.gcf().autofmt_xdate() + plt.xlabel("Datetime") + plt.title("Cost Function {} {}".format(obs, cty, var)) + plt.savefig(png) + plt.show() + + print ("figure saved in: ", png) + + + +if __name__=="__main__": + main() diff --git a/scripts/JOdiag4DVar.py b/scripts/JOdiag4DVar.py new file mode 100755 index 0000000..0ae45e7 --- /dev/null +++ b/scripts/JOdiag4DVar.py @@ -0,0 +1,345 @@ + +import datetime +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +from matplotlib.dates import DayLocator, HourLocator +from matplotlib.ticker import (MultipleLocator, + FormatStrFormatter, + AutoMinorLocator) +import numpy as np +import os +import pandas as pd +import re +######################################### +years = ['2023'] +months = ['05'] +#days = ['01','02','03','04','05','06','07', '08', '09','10', '11', '12', '13', '14', '15'] +days = ['15', '16'] +#days = ['15'] +cycles = ['00','03','06','09','12','15','18','21'] +#cycles = ['00','03','06','12','15','18'] +#cycles = ['06','09','12','21'] +path = '/scratch/nkim/logFilesFrom4dvar_v22' +#out_path = '/scratch/nkim/logFilesFrom4dvar_v22/JOstatistics' +out_path = '/scratch/nkim/' +keep_JOtables = True +figs = True +Variable = 'U' +Obstype = 'AMV' +################################# +def main (): + + print (">>>> Diagnostics using cost function information from log files <<<<<") + + Codetype = SetCodetype ( Obstype, Variable) + +# Sanity checks + if (np.size(Codetype) == 0 or np.size(Variable) == 0): + print ("ERROR: check your settings for Codetype and Variable") + print("Codetype: ", Codetype) + print("Variable : ", Variable) + exit() + + + num_code = 0 + + for cty in Codetype: + num_code = num_code + 1 + JoGlobalAll = [] + JoCtypeAll = [] + date_dfAll = [] + + for year in years: + for month in months : + + for day in days: + for cy in cycles: + fileLog = '{}/HM_Date_{}{}{}{}.html'.format(path,year,month,day,cy) + fileJOtable = '{}/JOtable_{}{}{}{}.txt'.format(out_path,year,month,day,cy) + + print(">>>>>>> Reading the log file in: ", fileLog) + #extractJOtable(fileLog, fileJOtable) + + print(">>>>>>> Reading the JoGlobal") + Global = readGlobal(fileJOtable) + df1 = pd.DataFrame(Global['JoGlobal']) + df2 = pd.DataFrame(Global['Ntotal']) + dataListG = createDataList(df1, df2) + JoGlobal = pd.DataFrame(dataListG) + + d = datetime.datetime(int(year),int(month), + int(day),int(cy)).strftime("%Y-%m-%d %H") + date_df = pd.DataFrame ({"datetime": [d]} ) + JoGlobalAll.append(JoGlobal) + date_dfAll.append(date_df) + + print(">>>>>>> Reading the Jo{}: ".format(cty)) + Observations = readObservations(fileJOtable, cty, Variable) + df3 = pd.DataFrame(Observations['JoOBS']) + df4 = pd.DataFrame(Observations['NOBS']) + dataListC =createDataList(df3, df4) + JoCTy = pd.DataFrame(dataListC) + JoCtypeAll.append(JoCTy) + + if keep_JOtables is False: + print(">>>>>>> deleting: ", fileJOtable) + os.remove (fileJOtable) + + print(">>>>>>> Creating the list of Jo(s)") + JoGlobalAll = pd.concat(JoGlobalAll, ignore_index=True, axis=0) + JoCtypeAll = pd.concat(JoCtypeAll, ignore_index=True, axis=0) + dateAll = pd.concat(date_dfAll, ignore_index= True, axis=0) + + + JoGlobalAll = pd.concat([dateAll, JoGlobalAll], axis=1) + JoCtypeAll = pd.concat([dateAll, JoCtypeAll], axis=1) + + print(">>>>>>> writing Jo(s) in :", out_path) + Jodf = writeJoGlobal(out_path, JoGlobalAll,"Global") + JoCTydf = writeJoGlobal(out_path, JoCtypeAll, "Observations_{}_{}".format(Obstype,num_code, Variable)) + print(">>>>>>> saving Jo(s) plots in :", out_path) + + if (figs): plotJo('JoCtype', out_path, JoCTydf, Obstype, num_code, cty, Variable) + + if (figs): plotJo('JoGlobal', out_path, Jodf,"Global","", "", "") + +####### +def SetCodetype ( Obstype, Variable): + + if (Obstype == 'SYNOP') : + Codetype = ['Codetype 14', 'Codetype 24', \ + 'Codetype 170', 'Codetype 182'] + if ( Variable != 'Z' ): + print ("{} variable not supported for SYNOP (only Z)".format(Variable)) + exit() + elif (Obstype == 'Aircraft'): + Codetype = ['Codetype 141', 'Codetype 146', 'Codetype 147'] + if ( Variable != 'U' and Variable != 'T'): + print ("{} variable not supported for Aircraft (only U and T)".format(Variable)) + exit() + + elif (Obstype == 'Radiosondes'): + Codetype = ['Codetype 109 === BUFR Land TEMP' ] + if ( Variable != 'U' and Variable != 'T' and Variable != 'Q'): + print ("{} variable not supported for Radiosondes (only U, T and Q)".format(Variable)) + exit() + elif (Obstype == 'AMSUA'): + Codetype = ['Codetype 210 === metop 1 3 SENSOR=amsua', \ + 'Codetype 210 === metop 3 5 SENSOR=amsua', \ + 'Codetype 210 === noaa 18 209 SENSOR=amsua', \ + 'Codetype 210 === noaa 19 223 SENSOR=amsua'] + if ( Variable != 'RAD'): + print ("{} variable not supported for AMSUA (only RAD)".format(Variable)) + exit() + elif (Obstype == 'MHS'): + Codetype = ['Codetype 210 === metop 1 3 SENSOR=mhs', \ + 'Codetype 210 === metop 3 5 SENSOR=mhs', \ + 'Codetype 210 === noaa 19 223 SENSOR=mhs'] + if ( Variable != 'RAD'): + print ("{} variable not supported for MHS (only RAD)".format(Variable)) + exit() + elif (Obstype == 'ATMS'): + Codetype = ['Codetype 210 === jpss 0 224 SENSOR=atms', \ + 'Codetype 210 === noaa 20 225 SENSOR=atms'] + if ( Variable != 'RAD'): + print ("{} variable not supported for ATMS (only RAD)".format(Variable)) + exit() + elif (Obstype == 'MWHS2'): + Codetype = ['Codetype 210 === fy3 4'] + if ( Variable != 'RAD'): + print ("{} variable not supported for MWHS2 (only RAD)".format(Variable)) + exit() + elif (Obstype == 'IASI' ): + Codetype = ['Codetype 210 === metop 1 5 SENSOR=iasi',\ + 'Codetype 210 === metop 3 5 SENSOR=iasi'] + if ( Variable != 'RAD'): + print ("{} variable not supported for IASI (only RAD)".format(Variable)) + exit() + elif (Obstype == 'ASCAT'): + Codetype = ['Codetype 139 === METOP-B ASCAT', 'Codetype 139 === METOP-C ASCAT'] + if ( Variable != 'U10'): + print ("{} variable not supported for ASCAT (only U10)".format(Variable)) + exit() + elif (Obstype == 'LIMB'): + Codetype = ['Codetype 250 === METOP-1', 'Codetype 250 === METOP-3',\ + 'Codetype 250 === SPIRE'] + if ( Variable != 'RO'): + print ("{} variable not supported for LIMB (only RO)".format(Variable)) + exit() + elif (Obstype == 'AMV'): + Codetype = ['Codetype 90 === METOP 5 METHOD=IR',\ + 'Codetype 90 === METEOSAT 56 METHOD=WVCL1',\ + 'Codetype 90 === METEOSAT 56 METHOD=WVCL2', \ + 'Codetype 90 === METEOSAT 56 METHOD=IR3',\ + 'Codetype 90 === METEOSAT 57 METHOD=WVCL1',\ + 'Codetype 90 === METEOSAT 57 METHOD=WVCL2',\ + 'Codetype 90 === METEOSAT 57 METHOD=IR3',\ + 'Codetype 90 === dual-MTOP 852 METHOD=IR'] + if ( Variable != 'U' and Variable != 'T'): + print ("{} variable not supported for AMV (only U and T)".format(Variable)) + exit() + else: + print ("{} Obstype not supported.".format(Obstype)) + + + return Codetype +#### +def extractJOtable(fileLog, fileJOtable): + + + in_JOtable = False + mylines = [] + fh = open(fileJOtable, 'w') + with open (fileLog, 'rt') as myfile: + for myline in myfile: + # add its contents to mylines. + # looks for the start of the 1st 4DVminim + if( "Diagnostic JO-table (JOT) MINIMISATION JOB T0089 NCONF= 131 NSIM4D= 0 NUPTRA= 0" in myline): + in_JOtable = True + elif("End of JO-table (JOT)" in myline) : + in_JOtable = False + + # looks for the end of the 1st 4DVminim + elif ("Diagnostic JO-table (JOT) MINIMISATION JOB T0089 NCONF= 131 NSIM4D= 999 NUPTRA= 0" in myline): + in_JOtable = True + elif("End of JO-table (JOT)" in myline) : + in_JOtable = False + + # looks for the 1st trajectory run + if ("Diagnostic JO-table (JOT) TRAJECTORY JOB T0539 NCONF= 1 NSIM4D= 0 NUPTRA= 1" in myline): + in_JOtable = True + elif("End of JO-table (JOT)" in myline) : + in_JOtable = False + + # looks for the end of the 2nd 4DVminim + elif ("Diagnostic JO-table (JOT) MINIMISATION JOB T0179 NCONF= 131 NSIM4D= 999 NUPTRA= 1" in myline): + in_JOtable = True + elif("End of JO-table (JOT)" in myline) : + in_JOtable = False + + # looks for the 2nd trajectory run + if ("Diagnostic JO-table (JOT) TRAJECTORY JOB T0539 NCONF= 1 NSIM4D= 0 NUPTRA= 2" in myline): + in_JOtable = True + elif("End of JO-table (JOT)" in myline) : + in_JOtable = False + + #### add all the relevant JOT in a txt file (can be saved or not + if in_JOtable: + fh.write(myline) + fh.close() +#### +def readObservations (fileJOtable, code, var): + + in_JOtable = True + mylines = [] + rec1 = [] + rec2 = [] + nextLine = False + # + var = var + " " + + with open (fileJOtable, 'rt') as myfile: + # searches for codetypes + for mylines in myfile: + if (in_JOtable): + if ( re.search(code, mylines) ): + nextLine = True + if (nextLine): + if ( re.search (var , mylines)): + record = mylines.split() + rec1.append(int(record[1])) + rec2.append(float(record[3])) + nextLine = False + + if ("Jo Global :" in mylines ): + in_JOtable = False + + in_JOtable = True + + datalist = ({ + 'JoOBS' : [rec2], + 'NOBS' : [rec1] + + }) + + print(datalist) + return datalist + + +#### +def readGlobal(fileJOtable): + + mylines = [] + rec1 = [] + rec2 = [] + + with open(fileJOtable,'rt') as myfile: + for mylines in myfile: + if re.search("Jo Global", mylines): + record = mylines.split() + rec1.append(int(record[3])) + rec2.append(float(record[5])) + + datalist = ({ + 'JoGlobal' : [rec2], + 'Ntotal' : [rec1] + + }) + return datalist +####### +def createDataList(data1, data2): + + # sanity test + if (np.size(data2) == 0) : + data2[0] = 0 + data1[0] = 0 + data1[1] = 0 + data1[2] = 0 + data1[3] = 0 + data1[4] = 0 + + datalist = ({ + 'start4DVmn1': [data1[0]], + '4DVmn1': [data1[1]], + '4DVtr1': [data1[2]], + '4DVmn2': [data1[3]], + '4DVtr2': [data1[4]], + 'N' : [data2[0]] + }) + + return datalist +##### +def writeJoGlobal(out_path,JoGlobal,name): + + filename = '{}/Jo{}.csv'.format(out_path,name) + df = pd.DataFrame(JoGlobal) + df.to_csv(filename, index=False, sep =' ') + return df + +def plotJo(label,path, Jo, obs, num, cty, var): + + Jo = Jo.set_index("datetime") + print ("========= JO ========") + png = "{}/Jo{}{}{}.png".format(path,obs,num,var) + plt.close("all") + fig, ax = plt.subplots(figsize=(12, 6)) + ax.plot(Jo["start4DVmn1"].astype(float),label='start4DVmn1',linestyle='solid',color='blue') + ax.plot(Jo["4DVmn2"].astype(float), label = 'end4DVmn2',linestyle='solid',color='orange') + ax.plot(Jo["4DVtr2"].astype(float), label= '4DVtr2',linestyle='solid',color='red') + ax.set_ylim(0,1.5) + ax.set_ylabel("Jo/n") + ax2 = ax.twinx() + ax2.plot(Jo["N"].astype(int),label='N',linestyle='dashed',color='grey') + ax2.set_ylabel ("N total") + legend = ax.legend(loc='upper right', fontsize='x-small') + plt.gcf().autofmt_xdate() + plt.xlabel("Datetime") + plt.title("Cost Function {} {}".format(obs, cty, var)) + plt.savefig(png) + plt.show() + + print ("figure saved in: ", png) + + +if __name__=="__main__": + main() diff --git a/scripts/dep_stats.py b/scripts/dep_stats.py new file mode 100755 index 0000000..052da8f --- /dev/null +++ b/scripts/dep_stats.py @@ -0,0 +1,231 @@ +### +import argparse +import sys +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.colors as colors +from mpl_toolkits.basemap import Basemap + + +def main(Loop1, Loop2, Loop3): + INPUT1="{}".format(Loop1) + INPUT2="{}".format(Loop2) + INPUT3="{}".format(Loop3) +#load data from ASCII from odb +## From loop1: +### date , time, statid@hdr, lat, lon, endtime, timeslot, fg_depar, lores +### 0 1 2,3 4 5 6 7 8 9 +### where, loop1 start: fg_depar=FirstGuessDep and loop1 end: lores=dep after 4DVminim1 +##### +## From loop2: +### date , time, statid@hdr, lat, lon, endtime, timeslot, lores, hires +### 0 1 2,3 4 5 6 7 8 9 +### where, loop2 start hires = depa after 4DVtraj1 and loop2 end lores=dep after 4DVminim2 +### +### +### + print('++++++++++++++++++++++++++++ Reading ODB files : ', INPUT1) + df1 = pd.read_csv(INPUT1, sep=" ") + + df1.replace('NULL',np.NaN) + lat1 = pd.to_numeric(df1['lat@hdr'], errors='coerce') + lon1 = pd.to_numeric(df1['lon@hdr'], errors='coerce') +# + print('++++++++++++++++++++++++++++ Reading ODB files : ', INPUT2) + df2 = pd.read_csv(INPUT2, sep=" ") + + df2.replace('NULL',np.NaN) + lat2 = pd.to_numeric(df2['lat@hdr'], errors='coerce') + lon2 = pd.to_numeric(df2['lon@hdr'], errors='coerce') + +# + print('++++++++++++++++++++++++++++ Reading ODB files : ', INPUT3) + df3 = pd.read_csv(INPUT3, sep=" ") + + df3.replace('NULL',np.NaN) + lat3 = pd.to_numeric(df3['lat@hdr'], errors='coerce') + lon3 = pd.to_numeric(df3['lon@hdr'], errors='coerce') + +###################### 2D departures + timeslot1 = pd.to_numeric(df1['timeslot@timeslot_index'], errors='coerce') + fg_dep = pd.to_numeric(df1['fg_depar@body'], errors='coerce') + minim1 = pd.to_numeric(df1['lores@update_1'], errors='coerce') +## time series + + print('++++++++++++++++++++++++++++ Statistics by TSlot, loop1 ') + FG_dep, FG_std = getStatsTimeSlot(fg_dep, timeslot1) + MINI1_dep, MINI1_std = getStatsTimeSlot(minim1, timeslot1) + +###################### 2D departures + timeslot2 = pd.to_numeric(df2['timeslot@timeslot_index'], errors='coerce') + traj1 = pd.to_numeric(df2['hires@update_2'], errors='coerce') + minim2 = pd.to_numeric(df2['lores@update_2'], errors='coerce') +## time series + print('++++++++++++++++++++++++++++ Statistics by TSlot, loop2 ') + TRAJ1_dep, TRAJ1_std = getStatsTimeSlot(traj1, timeslot2) + MINI2_dep, MINI2_std = getStatsTimeSlot(minim2, timeslot2) + +###################### 2D departures + timeslot3 = pd.to_numeric(df3['timeslot@timeslot_index'], errors='coerce') + traj2 = pd.to_numeric(df3['hires@update_3'], errors='coerce') +## time series + print('++++++++++++++++++++++++++++ Statistics by TSlot, loop3 (end of loop2) ') + TRAJ2_dep, TRAJ2_std = getStatsTimeSlot(traj2, timeslot3) +# +############################################################# +### figures with statistics + + print('++++++++++++++++++++++++++++ figures ') + + + createPlots (lat1, lon1, fg_dep, minim1, + lat2, lon2, traj1, minim2, lat3, lon3, traj2, FG_dep, + MINI1_dep, TRAJ1_dep, MINI2_dep, TRAJ2_dep, FG_std, + MINI1_std, TRAJ1_std, MINI2_std, TRAJ2_std) + + +############################################################# +### ASCII files with time series + + print('++++++++++++++++++++++++++++ writes departures as a function of TSlot in an ASCII file ') + + createASCII (FG_dep,MINI1_dep, TRAJ1_dep, MINI2_dep, TRAJ2_dep, + FG_std,MINI1_std, TRAJ1_std, MINI2_std, TRAJ2_std) + +############################################################# +def getStatsTimeSlot(field, timeslot): + + Function_mean = np.empty(10) + Function_std = np.empty(10) + for i in range (0, 9): + Function_mean[i] = field.loc[timeslot == i + 1 ].mean() + Function_std[i] = field.loc[timeslot == i + 1 ].std() + + return Function_mean, Function_std +# +############################################################# +def createPlots (lat1, lon1, fg_dep, minim1, + lat2, lon2, traj1, minim2, lat3, lon3, traj2, + FG_dep, MINI1_dep, TRAJ1_dep, MINI2_dep, TRAJ2_dep, + FG_std, MINI1_std, TRAJ1_std, MINI2_std, TRAJ2_std): + +## +## Maps: +## minlim = colorbar min +## maxlim = colorbar max +## cm = colorbar chhosen +## Name = filename produced +## +########### user settings: +## + minlim = -2 + maxlim = 2 + cm = plt.cm.get_cmap('RdBu_r') + Name = "teste.png" + + + +## To basemap: + ulat = np.max(lat1) + llat = np.min(lat1) + rlon = np.max(lon1) + llon = np.min(lon1) + + + plt.close("all") + + fig = plt.figure(1,figsize=(12,8)) + + Tsl = np.array([el for el in range(1,11, 1)]) + + ax0 = plt.subplot2grid((2, 4), (0, 0), colspan=2) + ax0.set_title("departures") + ax0.set_xlim((0, 11)) + ax0.set_xticks((1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,10)) + ax0.plot(Tsl, FG_dep, color='blue', linestyle='--', marker='^', fillstyle='none', markersize=8, label='fg') + ax0.plot(Tsl, MINI1_dep, color='dodgerblue', linestyle='--', marker='x', fillstyle='none', markersize=8, label='minim1') + ax0.plot(Tsl, TRAJ1_dep,color='red', linestyle='--', marker='o', fillstyle='none', markersize=8,label='traj1') + ax0.plot(Tsl, MINI2_dep,color='darkred', linestyle='--', marker='d', fillstyle='none', markersize=8, label='minim2') + ax0.plot(Tsl, TRAJ2_dep,color='green', linestyle='--', marker='+', fillstyle='none', markersize=8, label='traj2') + ax0.grid(True) + ax0.legend(loc='upper right', borderaxespad=0.) + ax0.set_xlabel(" time slot number") + ax0.set_ylabel("mean") +## + ax1 = plt.subplot2grid((2, 4), (1, 0), colspan=2) + ax1.set_xlim((0, 11)) + ax1.set_xticks((1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,10)) + ax1.plot(Tsl, FG_std, color='blue', linestyle='--', marker='^', fillstyle='none', markersize=8, label='fg') + ax1.plot(Tsl, MINI1_std, color='dodgerblue', linestyle='--', marker='x', fillstyle='none', markersize=8, label='minim1') + ax1.plot(Tsl, TRAJ1_std,color='red', linestyle='--', marker='o', fillstyle='none', markersize=8,label='traj1') + ax1.plot(Tsl, MINI2_std,color='darkred', linestyle='--', marker='d', fillstyle='none', markersize=8, label='minim2') + ax1.plot(Tsl, TRAJ2_std,color='green', linestyle='--', marker='+', fillstyle='none', markersize=8, label='traj2') + ax1.grid(True) + ax1.legend(loc='lower left', borderaxespad=0.) + ax1.set_xlabel(" time slot number") + ax1.set_ylabel(" std") + + ax2 = plt.subplot2grid((2, 4), (0, 2), colspan=2) + plt.title('departures loop1: Fg-4DVtraj1') + m=Basemap(projection='cyl', llcrnrlon=llon, urcrnrlon=rlon, + llcrnrlat=llat, urcrnrlat=ulat, lat_0=llat, lon_0=llon, resolution='h') + x, y = m(lon1,lat1) + m.drawcoastlines() + m.drawparallels(np.arange(0.,360.,5.),labels=[1,0,0,0], fontsize=7) + m.drawmeridians(np.arange(-180.,180.,5.),labels=[0,0,0,1], fontsize=7) + im=m.scatter(x,y,c=fg_dep-traj1, marker='s',vmin = minlim , vmax = maxlim, s= 1,cmap=cm) + cbar=m.colorbar(im, location='bottom',pad="15%") + +## + ax3 = plt.subplot2grid((2, 4), (1, 2), colspan=2) + plt.title('departures loop2 : 4DVminim2-4DVtraj2') + m=Basemap(projection='cyl', llcrnrlon=llon, urcrnrlon=rlon, llcrnrlat=llat, urcrnrlat=ulat, lat_0=llat, lon_0=llon +, resolution='h') + x, y = m(lon2,lat2) + m.drawcoastlines() + m.drawparallels(np.arange(0.,360.,5.),labels=[1,0,0,0], fontsize=7) + m.drawmeridians(np.arange(-180.,180.,5.),labels=[0,0,0,1], fontsize=7) + im=m.scatter(x,y,c=minim2-traj2, marker='s',vmin = minlim, vmax = maxlim, s= 1,cmap=cm) + cbar=m.colorbar(im, location='bottom',pad="15%") + + plt.savefig(Name) + plt.show() + +############ +def createASCII (FG_dep, MINI1_dep, TRAJ1_dep, MINI2_dep, TRAJ2_dep, + FG_std, MINI1_std, TRAJ1_std, MINI2_std, TRAJ2_std): + +## fName = ASCII file name where statistics should be written +## +########### user settings: + fName = "teste.txt" + + Tsl = np.array([el for el in range(1,11, 1)]) + + a = np.array([Tsl, FG_dep, MINI1_dep, TRAJ1_dep, MINI2_dep, TRAJ2_dep, FG_std,MINI1_std, TRAJ1_std, MINI2_std, TRAJ2_std]) + b = np.transpose(a) + HEADER = 'Tsl FG 4DVminim1 4DVtraj1 4DVminim2 4DVtraj2 FG_std 4DVminim1_std 4DVtraj1_std 4DVmini2_std 4DVtraj2_std' + + np.savetxt(fName, b, fmt=['%0.2d','%2.4f', '%2.4f','%2.4f','%2.4f','%2.4f','%2.4f','%2.4f','%2.4f','%2.4f','%2.4f'], header=HEADER, delimiter=" ", newline='\n') + + return + +############################################################# +## +def create_parser(): + + parser = argparse.ArgumentParser() + parser.add_argument("Filename1", help = "name and path 1st loop ASCII file ($SCRATCH/)", type = str) + parser.add_argument("Filename2", help = "name and path 2nd loop ASCII file ($SCRATCH/)", type = str) + parser.add_argument("Filename3", help = "name and path 3rd loop ASCII file ($SCRATCH/)", type = str) + + return parser + +if __name__ == "__main__": +# + parser = create_parser() + args = parser.parse_args() + main(args.Filename1, args.Filename2, args.Filename3) + +