diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b9eda9..86c95c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,15 +19,21 @@ endif() set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} ${bounds}") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${dialect}") + # -# Compile. +# .Compile # +# compilation of readgsa.F90 routine as a python extention +add_subdirectory(src) + add_executable(jbdiagnose.x src/jbdiagnose.F90) add_executable(cv_header_list.x src/cv_header_list.F90) add_executable(prepdiacov.x src/prepdiacov.F90) add_executable(clddet_analyzer.x src/clddet_analyzer.F90) + + # install executables and scripts install (TARGETS jbdiagnose.x PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ @@ -46,6 +52,8 @@ install (TARGETS clddet_analyzer.x DESTINATION libexec) + add_subdirectory (scripts) +add_subdirectory (modules) add_subdirectory (share) add_subdirectory (tests) diff --git a/modules/CMakeLists.txt b/modules/CMakeLists.txt new file mode 100644 index 0000000..b8da4aa --- /dev/null +++ b/modules/CMakeLists.txt @@ -0,0 +1,30 @@ +function (install_modules ) + foreach (file ${ARGV}) + get_filename_component (name_without_extension ${file} NAME_WE) + install (FILES ${file} + PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ + DESTINATION modules ) + endforeach () +endfunction() + + +# __init__ file +install_modules (__init__.py) +# Py modules For tuneBR +install_modules (readgsa.so + gsacov.py + read_odb.py + config_env.py + sigma_bo.py ) + +# Py modules for obstool +install_modules (build_sql.py + conv_stats.py + handle_df.py + io_base.py + setting.py + utils.py + obstype_info.py) + + + diff --git a/modules/__init__.py b/modules/__init__.py new file mode 100755 index 0000000..ed693ac --- /dev/null +++ b/modules/__init__.py @@ -0,0 +1,17 @@ +#-*- cvoding: utf-8 -*- + +# TuneBR modules +from . import gsacov +from . import config_env +from . import sigma_bo +from . import read_odb + + +# Obstool modules +from . import conv_stats +from . import build_sql +from . import obstype_info +from . import io_base +from . import setting +from . import utils + diff --git a/modules/build_sql.py b/modules/build_sql.py new file mode 100755 index 0000000..126925b --- /dev/null +++ b/modules/build_sql.py @@ -0,0 +1,153 @@ +# -*- coding:utf-8 -*- +import os ,sys +from odb4py.utils import SqlParser + + + +class SqlHandler: + """ + Class : @ Parse the different observation dictionaries and build the + SQL query according to obstype, codetype, varno , sensor and level range + @ Checks the SQL statement before sending. + Returns : The sql statement according to variable and varno + + Methods : BuildQuery + CheckQuery + + + """ + def __init__(self ): + + # For the moment ! + self.obstool_select=None + + #self.derozier_select=None + #self.jarvinen_select=None + return None + + + def BuildQuery(self , **kwarg ): + self.cols =",".join(kwarg["columns" ])[1:] + self.tabs =",".join(kwarg["tables" ]) + self.obs_ =kwarg["obs"] + self.varno =kwarg["obsvano"] + self.ctype =kwarg["codetype"] + self.obst =kwarg["obstype"] + self.levels =kwarg["lrange"] + self.vtype =kwarg["vertco_type"] + self.vert_coord=kwarg["vertco"] + self.sensor =kwarg["sensor" ] + self.other =kwarg["remaining_sql"] + + + # Obstool columns and tables + self.obstool_select="SELECT "+self.cols+" FROM " +self.tabs + + # Check obstype + type_list=[] + if self.obst is not None: + if isinstance (self.obst, list ): + for tp in self.obst: + if tp is not None: + type_list.append (" obstype =="+str( tp ) ) + elif isinstance (self.obst, int ): + type_list.append (" obstype =="+str( self.obst ) ) + else: + print(self.obs_, " : obstype must be integer or a list" ) + sys.exit() + + # Check varno + varno_list=[] + if self.varno is not None: + if isinstance (self.varno , list ) : + for vr in self.varno: + if vr is not None: + varno_list.append (" varno =="+str( vr ) ) + elif isinstance( self.varno , int ): + varno_list.insert(0, " varno =="+ str(self.varno) ) + + else: + print(self.obs_ , ": varno must be integer ot list" ) + sys.exit() + + # Check codetype + """ + Not needed for the moment ! + ctype_list=[""] + if self.ctype is not None: + if isinstance( self.ctype, list ): + for ct in self.ctype: + if ct is not None: + ctype_list.append (" codetype =="+str( ct ) ) + elif isinstance( self.ctype , int ): + ctype_list.append( " codetype =="+str(self.ctype ) ) + else: + print(self.obs_, ": codetype must be integer or list " ) + sys.exit()""" + + + # Check sensor + sensor_list =[] + if self.sensor != None: + if isinstance (self.sensor , list ): + for s in self.sensor: + if s is not None: + sensor_list.append( " sensor=="+str(s) ) + elif isinstance (self.sensor , int ): + sensor_list.append( " sensor =="+str( self.sensor ) ) + else: + print(self.obs_, ": sensor must be integer or list" ) + sys.exit() + + + # Check level range + level_list=[] + if self.levels != None : + if len (self.levels) <2 or len( self.levels ) > 2: + print ( "Level range must have two limites [l1 , l2], but got argument:", self.levels ) + sys.exit(0) + elif not isinstance (self.levels[0], int ) or not isinstance ( self.levels [1], int): + print( "One or the both level limites in the list is/are not intger(s)" ) + sys.exit(0) + elif self.levels[1] < self.levels[0]: + print( "Bottom level limit can't be greater than top level limit" ) + sys.exit(0) + else: + l1 =str(self.levels[0] ) + l2 =str(self.levels[1] ) + level_cond="vertco_reference_1 >="+l1+" AND vertco_reference_1 <= "+l2 + level_list.append(level_cond) + + + where_cond_list=[] + for tp in type_list: + if len( type_list) !=0 : + where_cond_list.append( "".join( ( tp ) ) ) + for vr in varno_list: + if len( varno_list) !=0: + where_cond_list.append( " AND ".join( (tp, vr )) ) + if len(level_list) !=0: + for lv in level_list: + where_cond_list.append( " AND ".join( (tp, vr , lv )) ) + + query=" " + if self.other !=None and len(where_cond_list) == 0: + query= self.obstool_select +" WHERE "+self.other + + elif self.other !=None and len(where_cond_list) > 0: + for condition in where_cond_list: + query=self.obstool_select +" WHERE " + condition + " AND "+self.other + + elif self.other ==None and len(where_cond) > 0: + for condition in where_cond_list: + query= self.obstool_select +" WHERE " + condition + + return query + + + + def CheckQuery(self , query ): + p =SqlParser() + nfunc =p.get_nfunc ( query ) # N Columns=N pure columns - N functions in the query + sql =p.clean_string ( query ) # Remove any unprintable character ( \r\t ....etc) + return nfunc , sql diff --git a/modules/config_env.py b/modules/config_env.py new file mode 100755 index 0000000..992a59e --- /dev/null +++ b/modules/config_env.py @@ -0,0 +1,49 @@ +# -*- coding: utf-8 -*- +import os +import configparser + + +class TuneEnv: + """ + Class: Parse the sections and get the items as a dict + """ + def __init__(self, config ): + self.BeginDate= config.get('DATES', 'DATESTART') + self.EndDate = config.get('DATES', 'DATEEND' ) + #self.ntaks = config.get('DATES', 'NSLICE' , fallback=2) + self.outfile = config.get('OPTIONS', 'OUTFILE' , fallback = "sigmabo_ratios" ) + self.dca_cpu = config.get('OPTIONS' ,'DCA_CPU', fallback= 4 ) + self.llverb = config.getboolean('OPTIONS' ,'LLVERB' , fallback="false") + self.lwrite = config.getboolean('OPTIONS' ,'LWRITE' , fallback="false") + self.lplot = config.getboolean('OPTIONS' ,'LPLOT' , fallback="false") + self.prog_bar = config.getboolean('OPTIONS' ,'PROGBAR', fallback="false") + + + + + pp = config.items("PATHS") + self.path_opt = { str(opt[0]):str(opt[1]) for opt in pp } + + self.basedir =self.path_opt["BASEDIR"] + self.stabal =self.path_opt["STATFILE"] + self.odbpath =self.path_opt["ODBPATH"] + self.odb_template=self.path_opt["ODB_TEMPLATE"] + self.rows_path=self.path_opt["ROWS_PATH"] + self.tmp_dir =self.path_opt["WORKDIR"] + + mm = config.items("MODEL") + self.model_opt= { str(opt[0]):str(opt[1]) for opt in mm } + + self.cycle_inc =self.model_opt["CYCLE_INC"] + self.deltax =self.model_opt["DELTAX"] + self.nsmax =self.model_opt["NSMAX"] + self.nflev =self.model_opt["NFLEV"] + self.rednmc =self.model_opt["REDNMC"] + return None + + def __Dicts__(self ): + PathAttr = self.path_opt + ModelAttr= self.model_opt + return PathAttr ,ModelAttr + + diff --git a/modules/conv_stats.py b/modules/conv_stats.py new file mode 100755 index 0000000..179d142 --- /dev/null +++ b/modules/conv_stats.py @@ -0,0 +1,267 @@ +#-*-coding:utf-8 -*- +import os +import pandas as pd +import numpy as np +from collections import defaultdict + +# obstool modules +from .handle_df import GroupDf + + + +class DHLStat: + """ + Class : Compute the diffrent statistics + covariance, correlation and standard deviations + from FG1, FG2,OA1 and OA2 departures in observation space + + Following the methods : + Desroziers , Hollingsworth,Lonneberg ( DHL ) + + functions : + getCov + getSig + getCor + """ + def __init__(self, df, max_dist =100, bin_dist =10 , delta_t =60 ): + + # Stats + # Get the values in the middle of the bins (the maximum distance and bin interval could be reset !!) + # DEFAULT VALUES + self.dist_max=100 #Km + self.bin_int =10 #Km + + # Overwrite if different values + if max_dist != 100: + self.dist_max = max_dist + print("New value for maximum distance has been set. Maximum distance value = {} Km".format( self.dist_max ) ) + else: + print("Default value for maximum distance is used : {} Km".format( self.dist_max ) ) + + if bin_dist != 10 : + self.bin_int = bin_dist + print("New binning distance has been set. bin_interval= {} Km".format( self.bin_int) ) + else: + print("Default value for binning interval is used : {} Km".format( self.bin_int ) ) + + self.gp =GroupDf () + self.merged_df=df + return None + + + + def getCov(self , var , inplace=None ): + d1,d2, d3, d4, d5, d6, d7, d8 , dobs , _ , dt1, dt2 , var =self.gp.GroupByBins (self.merged_df, self.dist_max, self.bin_int ) + + # Bins distances + dist_list = d1.index.to_list() + + # Var list and period + lvar = [var]* len(dist_list) + lperiod = [ str(dt1) +"_"+str(dt2) for _ in range(len( dist_list )) ] + + # Nobs + nobs=list(dobs ) + + # HL (Holingsworth-LÖnnberg ) + t1 = np.divide (d5 ,dobs) + t2 = np.divide ( np.multiply(d2 ,d3 ) , np.power( dobs, 2 )) + cov_hl= np.subtract(t1, t2 ) + + # Desroziers B + tb1 = np.divide( d4, dobs ) + tb2 = np.divide( np.multiply(d1, d3 ) , np.power(dobs, 2 ) ) + drb_ =np.subtract( tb1, tb2 ) + cov_drB =np.subtract(cov_hl , drb_ ) + + # Desroziers R + tr1 = np.divide ( d4 , dobs ) + tr2 = np.divide (np.multiply( d1 , d3 ) , np.power(dobs , 2 ) ) + cov_drR = np.subtract(tr1 , tr2 ) + + # Covariances + df_cov = pd.DataFrame ({ + "dist" :dist_list, + "nobs" :nobs , + "var" :lvar , + "period" :lperiod , + "COV_HL" :cov_hl , + "COV_DR-B" :cov_drB , + "COV_DR-R" :cov_drR } ).astype( { + "nobs" :"int32" , + "var" :"category", + "period" :"category", + "COV_HL" :"float64" , + "COV_DR-B":"float64" , + "COV_DR-R":"float64" + }) + + df_cov = df_cov.reset_index(drop=True).dropna(inplace=False ).reset_index() + if inplace ==True : + return cov_hl , cov_drB, cov_drR # inside the class + else: + return df_cov + + + def getSig ( self , var , inplace =None ): + + d1,d2, d3, d4, d5, d6, d7, d8 , dobs , _ , dt1, dt2 , var =self.gp.GroupByBins (self.merged_df , self.dist_max ,self.bin_int ) + + # Bins distances + dist_list = d1.index.to_list() + + # Var list and period + lvar = [var]* len(dist_list) + lperiod = [ str(dt1) +"_"+str(dt2) for _ in range(len( dist_list )) ] + + # Nobs + nobs=list(dobs ) + + # SIGMA FG1 + st1=np.divide( d6 , dobs) + st2=np.divide( d2 , dobs)**2 + sigma_fg1=np.sqrt( np.subtract(st1, st2 )) + del st1 , st2 + + # SIGMA FG2 + st1=np.divide( d7 ,dobs ) + st2=np.divide( d3 ,dobs )**2 + sigma_fg2=np.sqrt( np.subtract(st1, st2 )) + del st1, st2 + + # SIGMA A1 + st1=np.divide( d8 ,dobs ) + st2=np.divide( d1 ,dobs )**2 + sigma_a1=np.sqrt( np.subtract(st1, st2 )) + df_sig =pd.DataFrame({ + "dist" :dist_list , + "nobs" :nobs , + "var" :lvar , + "period" :lperiod , + "sigma_fg1": sigma_fg1 , + "sigma_fg2": sigma_fg2 , + "sigma_a1" : sigma_a1 + }).astype( { "nobs" :"int32" , + "var" :"category", + "period" :"category", + "sigma_fg1":"float64" , + "sigma_fg2":"float64" , + "sigma_a1" :"float64" + }) + + df_sig = df_sig.reset_index(drop=True).dropna(inplace=False ).reset_index() + + + if inplace ==True: + return sigma_fg1 , sigma_fg2, sigma_a1 # used inside the class + else: + return df_sig + + + def getCor (self,var , inplace=None ): + d1 ,_, _ , _ , _ , _ , _ , _ , dobs , _ , dt1, dt2 , var =self.gp.GroupByBins (self.merged_df , self.dist_max ,self.bin_int ) + + # Bins distances + dist_list = d1.index.to_list() + + # Var list and period + lvar = [var]* len(dist_list) + lperiod = [ str(dt1) +"_"+str(dt2) for _ in range(len( dist_list )) ] + + # Nobs + nobs=list(dobs ) + + # COV & Sigma + cov_hl, cov_drB, cov_drR =self.getCov ( var , inplace =True ) + sfg1 , sfg2 , sa1 =self.getSig ( var , inplace =True ) + + + # H.L and DRZ CORRELATIONS + cor_hl =np.divide( cov_hl , np.multiply( sfg1 , sfg2 ) ) + cor_drB =np.divide( cov_drB, np.multiply( sfg1 , sfg2 ) ) + cor_drR =np.divide( cov_drR, np.multiply( sa1 , sfg2 ) ) + + + df_cor=pd.DataFrame({ + "dist" : dist_list , + "nobs" : nobs , + "var" : lvar , + "period" : lperiod , + "COR_HL" : cor_hl , + "COR_DR-R" : cor_drR , + "COR_DR-B" : cor_drB + } ).astype( { "nobs":"int32" , + "var" :"category", + "period":"category", + "COR_HL" :"float64" , + "COR_DR-R" :"float64" , + "COR_DR-B" :"float64" + }) + + + df_cor = df_cor.reset_index(drop=True).dropna(inplace=False ).reset_index() + + if inplace ==True: + return cor_hl , cor_drB , cor_drR + else: + return df_cor + + + def getStatFrame (self , var ): + d1 ,_, _ , _ , _ , _ , _ , _ , dobs , _ , dt1, dt2 , var =self.gp.GroupByBins (self.merged_df, self.dist_max ,self.bin_int ) + + # Bins distances + dist_list = d1.index.to_list() + + # Var list and period + lvar = [var]* len(dist_list) + lperiod = [ str(dt1) +"_"+str(dt2) for _ in range(len( dist_list )) ] + + # Nobs + nobs=list(dobs ) + + # Gather all Statistics + cov_hl, cov_drB, cov_drR =self.getCov ( var , inplace =True ) + cor_hl, cor_drB, cor_drR =self.getCor ( var , inplace =True ) + sigma_fg1,sigma_fg2,sigma_a1 =self.getSig ( var , inplace =True ) + + # DataFrame: Contains the Desroziers , Hollingsworth/Lonnberg Cov, Sigma and Corr + drhl_frame={ + "dist" :dist_list , + "nobs" :nobs , + "var" :lvar , + "period" :lperiod , + "COV_HL" :cov_hl , + "COV_DR-B" :cov_drB , + "COV_DR-R" :cov_drR , + "sigma_FG1":sigma_fg1 , + "sigma_FG2":sigma_fg2 , + "sigma_a1" :sigma_a1 , + "COR_HL" :cor_hl , + "COR_DR-R" :cor_drR , + "COR_DR-B" :cor_drB + } + + # Return statistics DF + # Remark that the dataframe is rounded to 6 decimal digits + # It's done to allow a consistent comparison with the values given by R code + stat_frame =pd.DataFrame ( drhl_frame ).astype({ + "nobs":"int32" , + "var" :"category", + "period":"category", + "COV_HL" :"float64" , + "COV_DR-R" :"float64" , + "COV_DR-B" :"float64" , + "sigma_FG1":"float64" , + "sigma_FG2":"float64" , + "sigma_a1" :"float64" , + "COR_HL" :"float64" , + "COR_DR-R" :"float64" , + "COR_DR-B" :"float64" + }).round(6) + + + + stat_frame = stat_frame.dropna(inplace=False ).reset_index(drop=True) + return stat_frame + # To be plotted !!! diff --git a/modules/gsacov.py b/modules/gsacov.py new file mode 100755 index 0000000..288152a --- /dev/null +++ b/modules/gsacov.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +import os, sys +from pathlib import Path +from math import sqrt, pi +import numpy as np + +# Py/Fortran module (Reading B matrix ) +from .readgsa import readcov + + + +class GSA: + """ + Class : CALL A FORTRAN ROUTINE readgsa.so + READ THE BACKGROUND COVARIANCES OF : + TEMPERATURE ,Q SPECIFIC ,VORTICIY AND DIVERGENCE + + : RETURNS VERTICAL AVERAGED SIGMA AND MEAN OF SIGMA + OF T,Q,VOR,DIV and Wind + """ + + dpar={2:"Temperature" , + 3:"Specific humidity", + 4:"Vorticity" , + 5:"Divergence" } + + def __init__ (self,paths , cfile , nsmax , nflev , deltax , lverb ,lwrite): + self.path =paths + self.cfile =cfile + self.nsmax =int (nsmax) + self.nlev =int (nflev) + self.deltax=float(deltax) + self.lverb =bool(lverb) + self.write =bool(lwrite) + + return None + + def Write2File (self , varname , pname , value ): + # OUT PATH ( BASEDIR/out) + outdir = Path(self.path["BASEDIR"]) / "out" + outdir.mkdir(parents=True, exist_ok=True) + + file_ = outdir / f"{varname}_{pname}" + + outfile=open(file_ , "w") + if isinstance ( value , list): + for i,j in enumerate (value): + outfile.write( str(i+1)+" "+str(j)+"\n") + outfile.close() + else: + outfile.write(str(value) ) + + + + def ComputeStd(self, pcov): + var=np.zeros((self.nlev ,self.nlev)) + std=[] + for i in range(0,self.nlev): + for j in range(0,self.nlev): + var[i,j]= sum( pcov[i,j,:] ) + + for i in range(self.nlev): + std.append(sqrt(var[i,i])) + mstd=sum( std )/len(std) + return std , mstd + + def GetGSA( self, kpar ): + nlev=None + if kpar == 2: + nlev=self.nlev + 1 + else: nlev=self.nlev + cov,kret = readcov( nlev, self.nsmax , kpar,self.cfile,self.lverb) + if kret != 0: + print( "FAILED TO READ MATRICES FOR PARAMETER :" , dpar[kpar] ) + else: + return cov + + def GetSigmaB(self, ipar): + tcov=self.GetGSA (2) # TEMPERATURE + qcov=self.GetGSA (3) # SPECIFIC Q + vcov=self.GetGSA (4) # VORTICITY + dcov=self.GetGSA (5) # DIVERGENCE + + if ipar == 2: + tstd_ver , mean_av_t=self.ComputeStd( tcov ) + if self.write == True: + self.Write2File ( "sigmab_profile", "t" , tstd_ver ) + self.Write2File ( "sigmab_mean" , "t" , mean_av_t ) + return tstd_ver , mean_av_t + if ipar == 3: + qstd_ver , mean_av_q=self.ComputeStd( qcov ) + if self.write == True: + self.Write2File ( "sigmab_profile", "q" , qstd_ver ) + self.Write2File ( "sigmab_mean" , "q" , mean_av_q ) + return qstd_ver , mean_av_q + if ipar == 4: + vstd_ver , mean_av_v=self.ComputeStd( vcov ) + if self.write == True: + self.Write2File ( "sigmab_profile", "vor" , vstd_ver ) + self.Write2File ( "sigmab_mean" , "vor" , mean_av_v) + return vstd_ver , mean_av_v + if ipar == 5: + dstd_ver , mean_av_d=self.ComputeStd( dcov ) + if self.write == True: + self.Write2File ( "sigmab_profile", "div" , dstd_ver ) + self.Write2File ( "sigmab_mean" , "div" , mean_av_d) + return dstd_ver , mean_av_d + + # COMPUTE WIND SPEED FROM VORTICITY AND DIV + # FOTRAN k is from 1 ..NSMAX + # PYTHON k is from 0 ..NSMAX-1 + + # INIT + covuv=np.zeros((vcov.shape[0],vcov.shape[1],vcov.shape[2])) + + # PREPARE LAGGED INDICES + zidx =[i for i in range(1,self.nsmax +1) ] + + for i in range(0,self.nlev): + for j in range(0,self.nlev): + for k in range(0,self.nsmax): + ik=zidx[k] + z=-(pi*ik /((2*self.nsmax+3)*self.deltax*1000.0))**2 + if i ==j: + covuv[i,j,k] = ((dcov[i,j,ik]+vcov[i,j,ik])/ ( -2.0*z )) + + uvstd_ver , mean_av_uv=self.ComputeStd( covuv ) + if ipar == 999: # SET 999 FOR UV (NO KPAR VALUE IN stabal FILE ) + if self.write == True: + self.Write2File ( "sigmab_profile", "uv" , uvstd_ver ) + self.Write2File ( "sigmab_mean" , "uv" , mean_av_uv ) + return uvstd_ver , mean_av_uv diff --git a/modules/handle_df.py b/modules/handle_df.py new file mode 100755 index 0000000..2e2bb97 --- /dev/null +++ b/modules/handle_df.py @@ -0,0 +1,115 @@ +#-*- coding:utf-8 -*- +import gc +import pandas as pd +import numpy as np +from itertools import repeat + + +from .io_base import DataIO + + + +class SplitDf: + """ + Class :Split and subset DF + """ + + def __init__ (self, df , var , cdtg , max_dist =None , bin_dist =None , time_int=None ): + self.max_dist = max_dist # Maxumum distance for binning in [Km] + self.bin_int = bin_dist # Binning interval in [Km] + self.time_int = time_int # Time between OmG/OmA pairs in [+/- min] + self.dist_df = df + self.var = var + self.cdtg = cdtg + return None + + def SubsetDf (self , bin_dist , max_dist ): + dbin = [0,1]+list(np.arange(bin_dist, max_dist + bin_dist , bin_dist )) + dlabel = [0 ]+list(np.arange(bin_dist, max_dist + bin_dist , bin_dist )) + ldist_bin= pd.cut( self.dist_df ["dist"], bins=dbin,labels=dlabel, right=True,include_lowest=True ) + + # DIVIDE BY DIST INTERVALS + df_dist = self.dist_df[self.dist_df ["dist"] <= max_dist ] + df_dist ["ldist"] = ldist_bin.astype("int32") + + # Prepare Desroziers /H.L statistics (departres cross product ) + stat_df = df_dist.assign( + AFGsqr = df_dist.OA1 * df_dist.FG2, + FGsqr = df_dist.FG1 * df_dist.FG2, + FGsqr1 = df_dist.FG1 * df_dist.FG1, + FGsqr2 = df_dist.FG2 * df_dist.FG2, + Asqr1 = df_dist.OA1 * df_dist.OA1, + A1F1 = df_dist.OA1 * df_dist.FG1 ) + + # Binning + spdf = ( stat_df.groupby("ldist", observed=True).agg( + Asum1 = ("OA1" , "sum"), + FGsum1 = ("FG1" , "sum"), + FGsum2 = ("FG2" , "sum"), + AFGsqr = ("AFGsqr", "sum"), + FGsqr = ("FGsqr" , "sum"), + num = ("dist" , "count"), + FGsqr1 = ("FGsqr1", "sum"), + FGsqr2 = ("FGsqr2", "sum"), + Asqr1 = ("Asqr1" , "sum") ).reset_index() ).round(6) + + # Add var and date + spdf["var" ]=[ self.var for v in range(len( spdf["num"] ) ) ] + spdf["date"]=[ pd.to_datetime( str(self.cdtg), format="%Y%m%d%H", errors="raise") for v in range(len( spdf["num"] ) ) ] + return spdf + + +class ConcatDf: + """ + Class :Concat DF + Returns a concatenated dataframes for one variable and the whole period + """ + + def __init__(self ): + return None + + def ConcatFromListe (self, sub_df ): + merged_dict={} + merged_df =pd.DataFrame () + for k , v in sub_df.items(): + if len(v) !=0 : + merged_df =pd.concat ( v ) + merged_dict[k]=merged_df.reset_index().drop(columns =["index"]) + return merged_dict + + +class GroupDf: + """ + Class : Group DF and produce the pre-stat for cov, cor and sigma computation + """ + def __init__(self ): + return None + + def GroupByBins(self, merged_df, max_dist, bin_int): + # Set dist and max dist + self.max_dist = max_dist + self.bin_int = bin_int + + # Binning at the middles of the intervals + d_bins = np.arange(0, max_dist + bin_int, bin_int) + d_label = np.arange(bin_int/2, max_dist, bin_int) + merged_df["DIST"] = pd.cut(merged_df["ldist"], bins=d_bins, labels=d_label, include_lowest=False ) + + # Group/sum by dist bins + g = merged_df.groupby("DIST", observed=True) + d1 = g["Asum1" ].sum() + d2 = g["FGsum1"].sum() + d3 = g["FGsum2"].sum() + d4 = g["AFGsqr"].sum() + d5 = g["FGsqr" ].sum() + d6 = g["FGsqr1"].sum() + d7 = g["FGsqr2"].sum() + d8 = g["Asqr1" ].sum() + dobs = g["num" ].sum() + + # Get var and period + var = merged_df["var" ].iloc[ 0] + dt1 = merged_df["date"].iloc[ 0] + dt2 = merged_df["date"].iloc[-1] + return d1, d2, d3, d4, d5, d6, d7, d8, dobs, d1.index, dt1, dt2, var + diff --git a/modules/io_base.py b/modules/io_base.py new file mode 100755 index 0000000..02525a3 --- /dev/null +++ b/modules/io_base.py @@ -0,0 +1,89 @@ +# -*-coding :utf-8 -*- +import os , sys +import pandas as pd +import numpy as np +from pathlib import Path + + + +class DataIO: + """ + Class: DataIO : Contains methods to write and read the statstics dataframes + for each cycle (in .csv files) + Methods : + FlushFrame + ReadFrame + + """ + def __init__(self, compression="lz4"): + # If binary format is used , feather is the best for compression + # Write frames into feather format file with either "lz4" or "zstd" compression + # We use ".csv" for the moment + self.compression = compression + if self.compression not in ["ls4", "zstd"]: + self.compression = "ls4" # Fallback to default + + def FlushFrame(self, df, subdir, cdtg, var, fid=None, verbose =None, filetype=None ): + """ + Flush a DataFrame to Feather format with lz4 compression algo'. + Naming based on var, cdtg + """ + vrb = verbose + if df is None or df.empty: + print("Warning : Empty DataFrame for var={} datatime={}, nothing written".format( var , cdtg ) ) + return None + + # Final directory: dbpath / subdir / cdtg + outdir = Path(subdir) / cdtg + outdir.mkdir(parents=True, exist_ok=True) + if fid is not None: + filename = "_".join( (var,cdtg,fid ) ) +".csv" + else: + filename = "_".join( (var,cdtg)) +".csv" + + # Dir+Filename + filepath = "/".join( (str(outdir) , str(filetype) +"_"+str(filename) ) ) + try: + df.to_csv(filepath, index = False , sep=",") + if vrb in [3]: + print( "The dataframe has been written into {} for var: {} , datetime: {}".format(os.path.basename(filepath) , var , cdtg )) + except Exception: + print( "ERROR writing dataframe : var: {} and datetime: {}".format( var , cdtg )) + return None + return str(filepath) + + + + + + def ReadFrame(self, filepaths , verbose =None ): + """ + Read files. If var is provided, filters only matching files. + Returns a single concatenated DataFrame. + """ + vrb = verbose + frames = [] + if isinstance( filepaths, str ): + df = pd.read_csv( filepaths , header='infer', engine='python', sep=",") + return df + + elif isinstance (filepaths, list ): + for fp in filepaths: + if not os.path.isfile (fp): + print("File not found : " , fp ) + frames.append( pd.DataFrame() ) + else: + df = pd.read_csv( fp , header='infer', engine='python', sep=",") + if vrb in [3]: + print( "Dataframe loaded with : {} rows)".format( df.shape[0] )) + frames.append(df) + + if not frames: + return None + + # Concatenate all chunks + all_df = pd.concat(frames, axis=0, ignore_index=True) + print( "Total: {} rows merged".format( all_df.shape[0] )) + return all_df + + diff --git a/modules/obstype_info.py b/modules/obstype_info.py new file mode 100755 index 0000000..8d900de --- /dev/null +++ b/modules/obstype_info.py @@ -0,0 +1,294 @@ +# -*-coding:utf-8 -*- +import sys , os + + +class ObsType: + """ + Class: Contains the variable names and their attributes + NOTE : ONLY THE CONV DATA ARE PROCESSED IN THIS VERSION ! + + Methods: + ConvDict + RenameVarno + SelectConv + """ + + def _init__(self): + # ObsType lists + # Init the observation lists : Conventional and Satem + self.conv_obs=[ + 'gpssol' , + 'synop' , + 'dribu' , + 'airep' , + 'airepl' , + 'radar' , + 'temp' , + 'templ' ] + + self.sat_obs = [ + 'amsua' , + 'amsub' , + 'atms' , + 'iasi' , + 'mwhs' , + 'msh' , + 'seviri' ] + return None + + def ConvDict(self) : + # CONVENTIONAL + self.obs_conv=[ + + + # MAYBE WILL BE CHANGED TO A LIST AND LOOP ! + # SURFACE + {"obs_name" : "gpssol", + "obstype" : 1 , + "codetype" : 110 , + "varno" : 128 , + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None } , + { "obs_name" : "synop_z", + "obstype" : 1 , + "codetype" : [11, 14, 170, 182] , + "varno" : 1, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "synop_v", + "obstype" : 1 , + "codetype" : [11, 14, 170, 182] , + "varno" : 42, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "synop_u", + "obstype" : 1 , + "codetype" : [11, 14, 170, 182] , + "varno" : 41, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "synop_h", + "obstype" : 1 , + "codetype" : [11, 14, 170, 182] , + "varno" : 58, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "synop_t", + "obstype" : 1 , + "codetype" : [11, 14, 170, 182] , + "varno" : 39, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "dribu_z", + "obstype" : 4 , + "codetype" : None, + "varno" : [1, 39, 41, 42], + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "dribu_t", + "obstype" : 4 , + "codetype" : None, + "varno" : 39, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "dribu_u", + "obstype" : 4 , + "codetype" : None, + "varno" : 41, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "dribu_v", + "obstype" : 4 , + "codetype" : None, + "varno" : 42, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + + { "obs_name" : "ascat", + "obstype" : 9 , + "codetype" : None, + "varno" : None, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None } , + # RADAR + { "obs_name" : "radar_rh", + "obstype" : 13, + "codetype" : None, + "varno" : 29, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "radar_dow", + "obstype" : 13, + "codetype" : None, + "varno" : 195 , + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + + + # UPPER AIR + { "obs_name" : "airep_t", + "obstype" : [2] , + "codetype" : None, + "varno" : [2] , + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "airep_u", + "obstype" : 2 , + "codetype" : None, + "varno" : 3, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "airep_v", + "obstype" : 2 , + "codetype" : None, + "varno" : 4 , + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "airepl_t" , + "obstype" : 2 , + "codetype" : None, + "varno" : 2, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [25000, 35000] }, + { "obs_name" : "airepl_u" , + "obstype" : 2 , + "codetype" : None, + "varno" : 3, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [25000, 35000] }, + { "obs_name" : "airepl_v" , + "obstype" : 2 , + "codetype" : None, + "varno" : 4, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [25000, 35000] }, + { "obs_name" : "temp_t", + "obstype" : 5, + "codetype" : None, + "varno" : 2, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "temp_u", + "obstype" : 5, + "codetype" : None, + "varno" : 3, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "temp_v", + "obstype" : 5, + "codetype" : None, + "varno" : 4, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "temp_q", + "obstype" : 5, + "codetype" : None, + "varno" : 7, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : None }, + { "obs_name" : "templ_t" , + "obstype" : 5 , + "codetype" : None, + "varno" : 2, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [40000, 60000]} , + { "obs_name" : "templ_u" , + "obstype" : 5 , + "codetype" : None, + "varno" : 3, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [40000, 60000]} , + { "obs_name" : "templ_v" , + "obstype" : 5 , + "codetype" : None, + "varno" : 4, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [40000, 60000]} , + { "obs_name" : "templ_q" , + "obstype" : 5 , + "codetype" : None, + "varno" : 7, + "vertco_reference_1": None, + "sensor" : None, + "level_range" : [40000, 60000] } ] + + + # Varno dict + self.varno_dict ={ 128 : "ztd" , # gpssol Zenith total delay + 1 : "z" , # geopotential + 42 : "u" , # 10m wind speed u + 41 : "v" , # 10m wind speed v + 58 : "h" , # 2 relative humidity + 39 : "t" , # 2 meter temperature + 2 : "t" , # T temperature (upper air) + 3 : "u" , # U component of wind (upper air) + 4 : "v" , # V component of wind (upper air) + 7 : "q" , # specific humidity + 29 : "rh" , # upper rh (radar) + 195 : "dw" # Dopp radial wind + } + + # Add Units + self.unit= { + 128 : "m" , # gpssol Zenith total delay + 1 : "mgp" , # geopotential + 42 : "m/s" , # 10m wind speed u + 41 : "m/s" , # 10m wind speed v + 58 : "%" , # 2 relative humidity + 39 : "K" , # 2 meter temperature + 2 : "K" , # T temperature (upper air) + 3 : "m/s" , # U component of wind (upper air) + 4 : "m/s" , # V component of wind (upper air) + 7 : "g/kg" , # specific humidity + 29 : "%" , # upper rh (radar) + 195 : "m/s" # Dopp radial wind + } + return self.obs_conv , self.varno_dict , self.unit + + + + def RenameVarno (self, string_var ): + dict_={} + obslist , dict_ = self.ConvDict() + obsname=[ dobs["obs_name"] for dobs in obslist ] + name = string_var.split("_")[0] + vr =int(string_var.split("_")[1] ) + vname =name + "_"+ dict_[vr] + return vname + + + def SelectConv(self, list_ ): + varobs=[] + if not isinstance ( list_ ): + print("Parameter list must be a python list") + sys.exit () + self.list = list_ + obs_list, var_dict , var_unit =self.ConvDict() + + diff --git a/modules/read_odb.py b/modules/read_odb.py new file mode 100755 index 0000000..753eb10 --- /dev/null +++ b/modules/read_odb.py @@ -0,0 +1,256 @@ +# -*- coding: utf-8 -*- +import os , sys +from pathlib import Path +import shutil +from datetime import date, timedelta , datetime +import pandas as pd +import tarfile + + +# odb4py +from odb4py.utils import SqlParser , OdbObject +from odb4py.core import odb_open , odb_dca ,odb_dict + + +# TuneBR modules +from .config_env import TuneEnv + +class Odb: + """ + Class : read_odb. Get the paths from TuneEnv object and extrcat odb rows + + Method: ReadOdbRows. read obs_error , obs-guess, obs-analysis departures + """ + + def __init__(self, Paths ): + self.paths =Paths + self.dbtype ="CCMA" # Can be hard coded since it's NOT possible to use ECMA in TuneBR + # Static query can be used here ! + self.sql_query="SELECT obstype,varno,an_depar,fg_depar,final_obs_error FROM hdr, body, errstat" + + self.ext_list=[ item [0] for item in shutil.get_archive_formats() ] + self.rows_path =Path (Paths["ROWS_PATH"] ) + + + def CheckTarball (self, tarpath ): + try: + with tarfile.open(str(tarpath) , 'r') as t: + try: + member = t.getmember(self.dbtype ) + return member.isdir() + except KeyError: + return False + except: + FileNotFoundError + print("**WARNING : tar file {} not found. ".format(tarpath)) + pass + + + + def CopyOdb (self , dt ): + basedir =Path (self.paths["BASEDIR"] ) + odbpath =Path (self.paths["ODBPATH"] ) + odb_template=Path (self.paths["ODB_TEMPLATE"] ) + tmpdir =Path (self.paths["WORKDIR"] ) + dtedir = tmpdir / dt + + # rreate tmpdir if it doesn't exist + tmpdir.mkdir(parents=True, exist_ok=True) + + # Delete and recreate the date dir + shutil.rmtree(dtedir, ignore_errors=True) + dtedir.mkdir(parents=True) + + # Decompress ODB archive + # (We assume that the CCMA is the first dir inside the tar file , + # if not, an Exception is raised ! ) + FileTemplate=odb_template + filename=str(FileTemplate).replace( "YYYYMMDDHH", dt) + tarfile = odbpath / filename + status = self.CheckTarball( tarfile ) + + if status == False: + raise Exception ("The ODB archive must contain a CCMA directory" ) + + # Which extension .? It is often ".tar" But never know ! + arc_ext =os.path.splitext(tarfile)[1][1:] + if arc_ext not in self.ext_list: + raise Exception ("Supported formats for archive unpacking :", ext_list) # Supports: zip, tar, gztar, bztar, xztar or zstdtar. + + if os.path.isfile(tarfile ): + print( "Decompress the ODB archive ...", tarfile ) + # Unpack in tmp dir + shutil.unpack_archive(tarfile, dtedir , arc_ext) + + + + + + + def CreateDca ( self, dates , ncpu ): + print("Prepare DCA files for all ODBs ...\n") + NCPU = int(ncpu ) + basedir =self.paths["BASEDIR"] + tmpdir= "/".join((basedir, "tmp" )) + # It is possible to have a string , be sure it s a list + if isinstance (dates , str ): + dates = [ dates ] + + # Loop over dates list + for dt in dates: + if len( dt ) != 10: + print("Malformatted start date. Must be YYYYMMDDHH\n") + sys.exit(1) + + #print( "Prepare ODB : Date ... {} \n".format( dt ) ) + self.CopyOdb ( dt ) + dtedir= "/".join((tmpdir , dt )) + # Get attributes and create DCA + dbpath="/".join( (dtedir ,"CCMA" )) + if os.path.isdir ( dbpath ): + # It can hapen to have CCMA directory but a missing CCMA.dd , .sch or .desc + # Better use try + conn = odb_open ( dbpath ) + try: + db = OdbObject (dbpath) + db_attrs= db.get_attrib () + tab_list= db_attrs["tables"] + ic = conn.odb_dca(database =dbpath , + dbtype ="CCMA" , + ncpu = NCPU , + extra_args ="-z -u -q ", # Means Update if existe,remove empty files and run in quite mode + tables =tab_list ) # Only the tables found in the CCMA + conn.odb_close() + + except: + FileNotFoundError + print( f"Missing meta-data files in ODB : {dbpath}" ) + else: + + # Missing odbs + print ("ODB file : " , dbpath , " is missing !\n") + pass + + + + + def OdbExtract (self, dates, rows_path, progress=False , verbose =False ): + if isinstance (dates , str ): + dates = [ dates ] + + # Get Basedir + basedir=self.paths["BASEDIR"] + os.chdir(basedir) + + # Loop over dates list + for dt in dates: + if len( dt) != 10: + print("Malformatted start date. Must be YYYYMMDDHH\n") + sys.exit(1) + + print( "Rows extraction, ODB date : " , dt ) + tmpdir= "/".join( (basedir, "tmp" ) ) + dbpath= "/".join( (tmpdir , dt , "CCMA") ) + + # ODB ENV VARIBALES + # Update the variables. They are visible only inside this scope 'for each datetime' + odb_env={ "ODB_SRCPATH_CCMA" :dbpath , + "ODB_DATAPATH_CCMA":dbpath , + "TO_ODB_ECMWF":"0" , + "ODB_STATIC_LINKING":"1" , + "ODB_CMA":"CCMA" , + "ODB_IO_METHOD":"4" , + "ODB_CTX_DEBUG":"0" , + "VERSION":"1" , + "DEGRE":"1" , + "DIRECT":"0" , + "F_RECLUNIT":"BYTE" } + # Export variables + for k , v in odb_env.items(): os.environ[k] =v + + # Save extracted rows + os.makedirs( self.rows_path , exist_ok=True) + + # prepre a csv file name + filename="_".join( ("odb_rows", dt )) + outfile ="/".join( (rows_path , filename )) + + # Fetch ODB rows + if os.path.isdir ( dbpath ): + # We don't need to extract if we rerun for the same period + # Skip the datetime if the odb rows are already there ! + if not os.path.isfile (outfile): + try: + # odb_dict method throws an exception if no rows returned + conn =odb_open( dbpath ) + data_dict =conn.odb_dict (database=dbpath, + sql_query =self.sql_query , + nfunc =0 , + fmt_float =10 , + queryfile =None , + poolmask =None , + pbar =progress , + verbose =verbose ) + + # Save as csv file ( flush ) + pd.DataFrame(data_dict).to_csv( outfile , + header=False , + index= False , + sep="," , + decimal='.' ) + conn.odb_close() + except: + RuntimeError + print("No data returned for ODB {} \n".format( dt ) ) + pass + else: + + print(f"Rows are already extracted for the ODB {dbpath} \n" ) + continue + else: + print(f"**WARNING : CCMA directory not found for the ODB {dbpath} \n") + + + def ReadOdbRows ( basedir, rows_path, cdtg, target ): + # If file exists else skip datetime (in case of a rerun ) + + infile="/".join(( rows_path , "odb_rows_"+cdtg )) + if os.path.isfile ( infile ): + # Load dataset from file + df = pd.read_csv( infile , sep="," , names=["obstype","varno","an_depar","fg_depar","final_obs_error"] ) + # Get column Series + obstype=df["obstype"] + varno =df["varno" ] + + # Temperature + t_an =df.query( "varno == 2" )["an_depar"] + t_fg =df.query( "varno == 2" )["fg_depar"] + t_err=df.query( "varno == 2" )["final_obs_error"] + + # Specific Q + q_an =df.query( "varno == 7" )["an_depar"] + q_fg =df.query( "varno == 7" )["fg_depar"] + q_err=df.query( "varno == 7" )["final_obs_error"] + + # Brightness Tb + tb_an =df.query( "varno == 119" )["an_depar"] + tb_fg =df.query( "varno == 119" )["fg_depar"] + tb_err=df.query( "varno == 119" )["final_obs_error"] + + # Wind speed U,V + uv_an =df.query( "varno == 3 or varno==4" )["an_depar"] + uv_fg =df.query( "varno == 3 or varno==4" )["fg_depar"] + uv_err=df.query( "varno == 3 or varno==4" )["final_obs_error"] + + # Obs errors + if target=="predef": # Returns predefined errors + return (t_err, tb_err , q_err, uv_err) + + # Returns diagnosed fg departures & analysis departures + elif target=="fg_diag": # fg_dep + return (t_fg , tb_fg , q_fg , uv_fg) + + elif target=="an_diag": # an_dep + return (t_an , tb_an , q_an , uv_an) + else: + return None diff --git a/modules/setting.py b/modules/setting.py new file mode 100755 index 0000000..3c6c755 --- /dev/null +++ b/modules/setting.py @@ -0,0 +1,148 @@ +# -*- coding:utf-8 -*- +import os , sys +from datetime import datetime ,timedelta + +# Obstool modules +from .obstype_info import ObsType + + +class Setting: + """ + class : contains methods to set datetime period and varlist for obstool + Methods : set_period + set_obs_list + """ + def __init__ (self ): + self.obt = ObsType () + + def set_period(self, bdate , edate , cycle_inc=3 ): + # Create datetime liste + if not isinstance (bdate, str) or not isinstance ( edate , str): + btype, etype = type(bdate ) , type( edate) + print("Start and end dates period argument must be strings. Got {} {} ".format( btype, etype ) ) + sys.exit(0) + + if len( bdate) != 10: print("Malformatted start date") ; sys.exit(1) + if len( edate) != 10: print("Malformatted end date") ; sys.exit(2) + period=[] + bdate =datetime.strptime( bdate , "%Y%m%d%H") + edate =datetime.strptime( edate , "%Y%m%d%H") + delta =timedelta(hours=int(cycle_inc)) + while bdate <= edate: + strdate=bdate.strftime("%Y%m%d%H") + period.append( strdate ) + bdate += delta + return period + + + def set_obs_list(self, obs_list ): + obsdict , obsvarno , units = self.obt.ConvDict() + if len(obs_list ) == 0 : + print("Can not process empty observation liste. At least one obstype is required" ) + sys.exit(1) + + if not isinstance ( obs_list , list): + print("Observation types variable expects liste." ) + sys.exit(2) + + + obs_name =[ item["obs_name"] for item in obsdict ] + for obs in obs_list: + if obs.strip() not in obs_name: + print("Can t process the given obstype {}. Not predefined in the list of conventional obstypes".format( obs )) + sys.exit(3) + + # Obs attributes + obs_names=[] + codetype=[] + obstype =[] + varno =[] + lrange =[] + vertco =[] + sensor =[] + for item in obsdict: + obsname = item["obs_name"] + if obsname in obs_list: + obs_names.append(item["obs_name"]) + varno.append (item["varno" ]) + codetype.append (item["codetype"]) + obstype.append (item["obstype" ]) + lrange.append (item["level_range"]) + vertco.append (item["vertco_reference_1"]) + sensor.append (item["sensor" ] ) + return obs_names , codetype , obstype , varno , lrange , vertco ,sensor + + + + +# TODO +# class Satem: +# def __init__(): +# pass +# DATA CETEGORIES SATEM +# What we need ? +# SATEM +# NAME : sensor +# AMSUA : 3 +# AMSUB : 4 +# MHS :15 +# IASI :16 +# ATMS :19 +# MWHS :73 +# SEVIRI :29 +# return None + + + +# CONV OBSTYPE LIST +# obs varno +# gpssol : 128 +# synop : 1,29,58,41,42,7,39 +# dribu : +# airep : 2,3,4 +# airepl : 2,3,4 +# radar : 29,195 +# temp : 2,3,4,39,58,41,42 +# templ : 2,3,4,39,58,41,42 + +class Conv: + """ + class : instantiate Conv obs category and corresponding SQL query + """ + def __init__(self ): + # ---> CAUTION !: NEVER CHANGE THE ORDER OF COLUMNS + # TO ADD NEW ONEs , ONE CAN DO ONLY AN APPEND ! + self.cols =[ "", + "obstype" , # 1 + "codetype" , # 2 + "statid" , # 3 + "varno" , # 4 + "degrees(lat)" , # 5 + "degrees(lon)" , # 6 + "vertco_reference_1" , # 7 + "date" , # 8 + "time@hdr" , # 9 + "an_depar" , # 10 + "fg_depar" , # 11 + "obsvalue" , # 12 + ] + + # Considered tables & additional sql statement + ODB_CONSIDER_TABLES="/hdr/desc/body/radar/radar_body" + self.tables = ["hdr","desc","body",] #"errstat", "radar_body"] + self.tbl_env = "/".join( self.tables ) + self.other_sql = " (an_depar is not NULL) AND (fg_depar is not NULL) AND (datum_event1.fg2big@body = 0)" + + # ObsType list + self.conv_obs=[ 'gpssol' , + 'synop' , + 'dribu' , + 'airep' , + 'airepl' , + 'radar' , + 'temp' , + 'templ' ] + obstype = ObsType() + self.obs, self.varno , self.unit = obstype.ConvDict() + return None + diff --git a/modules/sigma_bo.py b/modules/sigma_bo.py new file mode 100755 index 0000000..7de72f1 --- /dev/null +++ b/modules/sigma_bo.py @@ -0,0 +1,359 @@ +# -*- coding: utf-8 -*- +import sys +import os +import pandas as pd +from math import sqrt +from statistics import stdev , mean + +from .read_odb import Odb +from .config_env import TuneEnv + + +# BASED ON THE PROGRAM ratio.F90 +# IN THE ORIGINAL VERSION B. Strajnar 2010/10/26 + +class Predef: + """ + Class : Returns predefined SIGMAO for t , bt, q and ke + """ + + def __init__(self , paths , dates , lverb , lwrite): + self.basedir=paths['BASEDIR']; + self.rows_path= paths["ROWS_PATH"] + self.lverb =lverb ; + self.lwrite =lwrite; + # INIT LISTS + self.psigma=[] ; + self.dates =dates ; + self.rabso =2147483647.0 ; # Missing data Indicator MDI + + + + def Write2File (self , varname , pname , value ): + # Output path ( BASEDIR/out) + os.makedirs(self.basedir+"/out" , exist_ok=True) + file_=self.basedir+"/out/"+varname+"_"+pname + + outfile=open(file_ , "w") + if isinstance ( value , list): + for i,j in enumerate (value): + outfile.write( str(i+1)+" "+str(j)+"\n") + outfile.close() + else: + outfile.write(str(value) ) + + + + + + # Returned params from ReadOdbRows + def GetSigmaP(self , param): + means=[] + # Get according to their index + if param=="t" : idx=0 + elif param=="bt": idx=1 + elif param=="q" : idx=2 + elif param=="ke": idx=3 + else: + print("Unkown parameter -->",param , ", possible short names : t , bt , q and ke" ) + sys.exit () + for dt in self.dates: + prows= Odb.ReadOdbRows (self.basedir, self.rows_path , dt , "predef") + if isinstance( prows , tuple ): + if isinstance( prows[idx], pd.Series): + pred_rows = prows[idx].dropna().to_list() # Clean NaN and convert to list + self.psigma= self.psigma +pred_rows + + # Get mean for each datetime + if len( self.psigma) != 0: + nobs_p = len( self.psigma ) + dte = dt[0:8] + hh = dt[8:10] + dh = dte+" "+hh+":00:00" + means.append( dh+" "+str( mean(self.psigma))+" "+str(nobs_p) ) + else: + nobs_p = 0 + dte= dt[0:8] + hh = dt[8:10] + dh =dte+" "+hh+":00:00" + means.append( dh+" "+"None"+" "+str(nobs_p)) + + + if self.lwrite == True: + self.Write2File ("so_pred_means_vs_date" ,param , means ) + + if len(self.psigma) != 0: + so_pred=sum( self.psigma)/len(self.psigma) + if self.lwrite == True: self.Write2File ("so_pred_mean" ,param ,so_pred ) + else: + so_pred=None + if self.lwrite == True: self.Write2File ("so_pred_mean" ,param ,so_pred ) + return so_pred + + + + + + +class Diag: + """ + Class : + # Returns SIGMAO AND SIGMAB diagnostics using the "obs differences" method + # METHOD : + # HBH^T = d(b,a) * d(b,o)^T + # R = d(o,a) * d(b,o)^T + """ + + def __init__(self ,paths , dates, lverb , lwrite): + self.basedir=paths["BASEDIR"] + self.rows_path = paths["ROWS_PATH"] + self.lverb =lverb + self.lwrite =lwrite + + # INIT LISTS + self.sigb =[] ; + self.sigo =[] ; + self.sigb_out=[] ; + self.sigo_out=[] ; + self.dates =dates ; + self.rabso = 2147483647.0 ; # Missing data + + def Write2File (self , varname , pname , value ): + # OUT PATH ( BASEDIR/out) + os.makedirs(self.basedir+"/out" , exist_ok=True) + file_=self.basedir+"/out/"+varname+"_"+pname + outfile=open(file_ , "w") + if isinstance ( value , list): + for i,j in enumerate (value): + outfile.write( str(i+1)+" "+str(j)+"\n") + outfile.close() + else: + outfile.write(str(value) ) + + + def ComputeSigmab(self,fg_dep , an_dep): + s=0 + if len(fg_dep )!=0 and len(an_dep) !=0: + for i in range(len(fg_dep )): + diff=fg_dep[i] - an_dep[i] + s =s+ diff * fg_dep[i] + sigmab=sqrt(s/len(fg_dep)) + else: + sigmab=None + return sigmab + + + def ComputeSigmao(self,fg_dep , an_dep): + s=0 + if len(fg_dep) != 0 and len(an_dep) != 0: + for i in range(len(fg_dep)): + s=s + (fg_dep[i]*an_dep[i] ) + sigmao = sqrt( s/len(fg_dep)) + else: + sigmao = None + return sigmao + + + # 0 = t, 1=bt, 2=q, 3=ke + def GetSigmaD(self, param): + if param=="t" : idx=0 ; + if param=="bt": idx=1 ; + if param=="q" : idx=2 ; + if param=="ke": idx=3 ; + ncase=[] + + for dt in self.dates: + # Get the pd.Series according to the idx ==> param + fg_rows =Odb.ReadOdbRows (self.basedir ,self.rows_path , dt, "fg_diag") + an_rows =Odb.ReadOdbRows (self.basedir ,self.rows_path , dt, "an_diag") + + # Check if it's pandas Series inside a tuple (as expected ) + if isinstance ( fg_rows , tuple ) and isinstance ( an_rows , tuple ): + if isinstance (fg_rows[idx], pd.Series ) and isinstance (an_rows[idx], pd.Series ) : + fg_rows = fg_rows [idx] + an_rows = an_rows [idx] + + # Concat an and fg + df = pd.concat([fg_rows , an_rows], axis=1) + + # Clean the df ( align the available values without NaN ) + df = df.dropna() + fg = df.iloc[:, 0].to_list() + an = df.iloc[:, 1].to_list() + + so = self.ComputeSigmao( fg, an ) # Observations + sb = self.ComputeSigmab( fg, an ) # Background + if so != None and sb != None: + self.sigo.append(so ) + self.sigb.append(sb ) + ncase.append(len(fg) ) # len(fg)=len(an) + + if self.lwrite==True: + if so !=None and sb !=None : # Available values + dte= dt[0:8] + hh = dt[8:10] + dh =dte+" "+hh+":00:00" + self.sigo_out.append(dh+" "+str(self.ComputeSigmao( fg, an ))+" "+str(len(fg) ) ) + self.sigb_out.append(dh+" "+str(self.ComputeSigmab( fg, an ))+" "+str(len(fg) ) ) + else: + dte= dt[0:8]; hh = dt[8:10] ; dh =dte+" "+hh+":00:00" + self.sigo_out.append(dh+" "+"None"+" "+str(len(fg) ) ) + self.sigb_out.append(dh+" "+"None"+" "+str(len(fg) ) ) + + if self.lwrite==True: + self.Write2File("so_diag_means_vs_date", param , self.sigo_out) + self.Write2File("sb_diag_means_vs_date", param , self.sigb_out) + + # Mean values + if len( self.sigo) !=0 and len(self.sigb ) !=0: + if self.lwrite==True: + self.Write2File("so_diag_mean", param , mean(self.sigo)) + self.Write2File("sb_diag_mean", param , mean(self.sigb)) + return mean(self.sigb) , mean(self.sigo), sum(ncase) + else: + if self.lwrite==True: + self.Write2File("so_diag_mean", param , self.rabso) + self.Write2File("sb_diag_mean", param , self.rabso) + return None , None , 0 + + + + +class Ratios: + def __init__(self,paths, Nobs, rednmc ,so_pred, so_diag , sb_pred, sb_diag, lwrite ): + self.basedir= paths["BASEDIR"] + self.ptot = sum(Nobs) + self.pt = Nobs[0] + self.pbt = Nobs[1] + self.pq = Nobs[2] + self.pke = Nobs[3] + self.so_tp = so_pred[0] ; self.so_td =so_diag[0] + self.so_btp = so_pred[1] ; self.so_btd=so_diag[1] + self.so_qp = so_pred[2] ; self.so_qd =so_diag[2] + self.so_kep = so_pred[3] ; self.so_ked=so_diag[3] + + self.sb_tp = sb_pred[0] ; self.sb_td =sb_diag[0] + self.sb_qp = sb_pred[1] ; self.sb_qd =sb_diag[2] + self.sb_kep = sb_pred[2] ; self.sb_ked=sb_diag[3] + + self.rednmc=rednmc + + self.rabso= 2147483647.0 + self.lwrite=lwrite + + + def Write2File (self , varname , pname , value ): + # OUT PATH ( BASEDIR/out) + os.system( "mkdir -p "+ self.basedir+"/out" ) + file_=self.basedir+"/out/"+varname+"_"+pname + + outfile=open(file_ , "w" ) + if isinstance ( value , list): + for i,j in enumerate (value): + outfile.write( str(i+1)+" "+str(j)+"\n") + outfile.close() + else: + outfile.write(str(value) ) + + + + + + def ComputeRatios(self, pred , diag , rednmc, target ): + if pred != None and diag != None : + if target == "sigmao": + ratio =diag/pred + return ratio + elif target == "sigmab": + ratio =diag/(pred * float(rednmc)) + return ratio + else: + ratio =None + return ratio + + + def RatioSo(self): + """ + The averaged ratio is computed by weighted sum of + the ratio of each observation subset + roav=sqrt(roq**2*float(pq)/float(ptot)+rot**2*float(pt)/float(ptot) + robt**2*float(pbt)/float(ptot)+roke**2*float(pke)/float(ptot)) + + """ + + target="sigmao" + rot = self.ComputeRatios( self.so_tp , self.so_td ,self.rednmc , target ) + robt= self.ComputeRatios( self.so_btp , self.so_btd ,self.rednmc , target ) + roq = self.ComputeRatios( self.so_qp , self.so_qd ,self.rednmc , target ) + roke= self.ComputeRatios( self.so_kep , self.so_ked ,self.rednmc , target ) + + if rot == None : rot =0. + if roq == None : roq =0. + if robt == None : robt=0. + if roke == None : roke=0. + + if self.ptot != 0: + rrot = rot **2* float(self.pt ) /float(self.ptot) + rrobt = robt**2* float(self.pbt) /float(self.ptot) + rroq = roq **2* float(self.pq ) /float(self.ptot) + rroke = roke**2* float(self.pke) /float(self.ptot) + + roav = sqrt( rroq +rrot + rrobt + rroke ) + return (rot , robt , roq , roke ,roav) + else: + print("Number of total observations is equal to zero") + return None + + if self.lwrite==True: + lines= "ro_t : "+str("%.4f" % rot) +" | Nobs_t : "+str(int(self.pt)) +"\n" \ + "ro_bt : "+str("%.4f" % robt) +" | Nobs_bt : "+str(int(self.pbt)) +"\n" \ + "ro_q : "+str("%.4f" % roq) +" | Nobs_q : "+str(int(self.pq) ) +"\n" \ + "ro_ke : "+str("%.4f" % roke) +" | Nobs_ke : "+str(int(self.pke)) +"\n" \ + "\n" \ + "ro_avg : "+str("%.4f" % roav) +" | Nobs_tot : "+str(int(self.ptot)) +"\n" + self.Write2File( "ratios" ,"so" , lines ) + #return rot , robt , roq , roke ,roav + + + + + def RatioSb(self): + """ + sb_*d: diagnosed sigmab + sb_*p: predefined sigmab + rb* : tuning ratio for sigmab for each parameter + rbav : weighted average ratio for sigmab + + rb=sbd/(sbp*rednmc) + rbav=sqrt(rbq**2*pq/ptot+rbt**2*pt/ptot+rbke**2*pke/ptot) + """ + + target="sigmab" + rbt = self.ComputeRatios( self.sb_tp , self.sb_td ,self.rednmc , target ) + rbq = self.ComputeRatios( self.sb_qp , self.sb_qd ,self.rednmc , target ) + rbke = self.ComputeRatios( self.sb_kep , self.sb_ked ,self.rednmc , target ) + + # AVERAGE + if rbt == None : rbt =0. + if rbq == None : rbq =0. + if rbke == None : rbke=0. + + if self.ptot != 0: + rrbt =rbt**2 *(self.pt / self.ptot ) + rrbq =rbq**2 *(self.pq / self.ptot ) + rrbke=rbke**2*(self.pke / self.ptot ) + rbav = sqrt( rrbt + rrbq + rrbke ) + + if self.lwrite == True: + lines= "rb_t : "+str("%.4f" % rbt) +" | Nobs_t : "+str(int(self.pt)) +"\n" \ + "rb_q : "+str("%.4f" % rbq) +" | Nobs_q : "+str(int(self.pq) ) +"\n" \ + "rb_ke : "+str("%.4f" % rbke) +" | Nobs_ke : "+str(int(self.pke)) +"\n" \ + "\n" \ + "rb_avg : "+str("%.4f" % rbav) +" | Nobs_tot : "+str(int(self.ptot)) +"\n" + self.Write2File( "ratios" ,"sb" , lines ) + + return (rbt , rbq , rbke , rbav) + else: + print( "Number of total innovations is equal to zero" ) + return None + diff --git a/modules/utils.py b/modules/utils.py new file mode 100755 index 0000000..84128e9 --- /dev/null +++ b/modules/utils.py @@ -0,0 +1,494 @@ +# -*-coding :utf-8 -*- +import os , sys , gc +import pandas as pd +import numpy as np +from pathlib import Path +from collections import defaultdict +import multiprocessing as mp +from multiprocessing import Pool , cpu_count, shared_memory +from itertools import chain + + +# odb4py mdules +from odb4py.utils import SqlParser , OdbObject +from odb4py.core import odb_dca, odb_open , odb_close ,odb_dict, odb_gcdist + + + +# Obstool, Desroziers & Jarvinen , Tools & modules +from .build_sql import SqlHandler +from .obstype_info import ObsType +from .setting import Setting , Conv +from .handle_df import * +from .io_base import DataIO + + + + +class DCAFiles: + """ + Class : Checks and creates the DCA (Direct Access Column )files + if they are not already in the ODB + If the DCA files are already there , this class is not called + Methods : + CheckDca + + """ + def __init__(self): + return None + + def CheckDca( self, dbpath , sub_base=None , verbose = False ): + # Prepare DCA files if not in ODB + try: # os.path.isdir( dbpath ): + db = OdbObject ( dbpath ) + dbname = db.get_attrib()["name"] + + if not os.path.isdir ("/".join( [dbpath , "dca"] ) ): + if verbose in [2,3 ]: + print( "No DCA files in {} 'directory'".format(dbpath ) ) + #env.OdbVars["CCMA.IOASSIGN"]="/".join( (dbname, "CCMA.IOASSIGN" ) ) + #env.OdbVars["ECMA.IOASSIGN"]="/".join( (dbname, "ECMA.IOASSIGN" ) ) + status = odb_dca ( database=dbpath , dbtype=dbname , ncpu=8 ) + if status < 0 : + print("Failed to create DCA files ... \n Another attempt will be done with odb_dict !" ) + else : + if verbose ==True : + print("DCA files already in database: '{}'".format( dbname ) ) + else: + pass + except: + FileNotFoundError + print("**WARNING : ODB path {} not found".format(dbpath)) + pass + + + + + +class OdbReader: + """ + Class : Prepare the query according to the user setting (Obs list , period etc) + The SQL query is sent directly to the ECMA or CCMA ODB . + + Returns data are as python dictionary + Methods : get_odb_rows + + """ + def __init__(self , dbpath ,type_ ): + # Path to the directory containing the ODB(s) + self.odb_path = dbpath + self.odb_type = type_ + self.vrb = 1 + if not os.path.isdir (self.odb_path): + print("ODB(s) directory '{}' not found.".format( self.odb_path )) + sys.exit(0) + + # Setting + self.st = Setting () + + # DCA Files + self.dca_f = DCAFiles() + + # CONV + self.conv = Conv() + self.conv_obs = self.conv.conv_obs + self.cols = self.conv.cols + self.tables = self.conv.tables + self.other_sql = self.conv.other_sql + + # IO + self.io = DataIO () + # SQL + self.sql =SqlHandler() + + # DF + self.rd =Rows2Df () + self.cnt =ConcatDf() + + # List of odb rows + self.list_rows =[] + + # return a list of dicts + self.dlist = defaultdict(list) + + def get_odb_rows (self, period , + obs_list , + max_dist , + bin_dist , + #file_io , later ! + file_path , + fextract , + cycle_inc =3, + pbar =False , + verbosity =0, + chunk_size = None + ): + + vrb=verbosity + if vrb not in [0,1,2,3]: + print("Min and max verbosity levels: 0 -> 3. Got : ", vrb ) + print("Fallback to default value: verbosity=", self.vrb ) + + + # Default dict (Collect dataframes with variable as a key ) + df_vars = defaultdict(list) + + # Diags period + period_interval= [ min( period ) , max(period) ] + + # ODB Paths + paths = [ "/".join( (self.odb_path , prd , self.odb_type )) for prd in period ] + + # Obs attributes + obs_names , codetype , obstype , varno , lrange , vertco ,sensor =self.st.set_obs_list(obs_list ) + + list_dict=[] + for i, cma_path in enumerate(paths): + + # Open ODB + conn= odb_open ( cma_path ) + + # UPDATE IOASSIGN , ODB_SRCPATH & ODB_DATAPATH + os.environ["IOASSIGN"]=cma_path+"/IOASSIGN" + os.environ["ODB_SRCPATH_CCMA"] =cma_path + os.environ["ODB_DATAPATH_CCMA"]=cma_path + os.environ["ODB_IDXPATH_CCMA" ]=cma_path + + # Check DCA directory (if not there they will be created ) + dca_f=DCAFiles() + dca_f.CheckDca ( cma_path ) + + if vrb == [ 2 , 3]: + print("ODB PATHS set to :") + print("IOASSIGN :",cma_path+"IOASSIGN" ) + print("ODB_SRCPATH_CCMA :",cma_path ) + print("ODB_DATAPATH_CCMA :",cma_path ) + print("ODB_IDXPATH_CCMA :",cma_path ) + + for jo, obs in enumerate(obs_names) : + query=self.sql.BuildQuery( columns =self.cols , + tables =self.tables , + obs =obs , + obstype =obstype [jo] , + obsvano =varno [jo] , + codetype =codetype [jo] , + lrange =lrange [jo] , + vertco_type ="height" , + vertco =vertco [jo] , + sensor =sensor [jo] , + remaining_sql =self.other_sql ) + nfunc , sql_query = self.sql.CheckQuery( query) + cdtg = period [i] + if vrb in [0, 1, 2, 3,4]: + print( "Process observation type: {} ODB date : {} ".format( obs , cdtg )) + query_file=None ; + poolmask = None ; + pool =None ; + float_fmt= 15 ; # 10 digits float values + verbose = False ; + pbar = False + + # Progress bar & Verbosity inside pyodb + # Progress bar is useful for huge ODBs + if vrb in [3,4]: + verbose= True + if vrb in [2,3,4]: + pbar = True + + + # Write odb rows once , if there is a rerun of the same period + # the odb extraction is skipped + filename = "_".join( ( "df_rows" , obs ,cdtg)) +".csv" + fpath = "/".join( (self.odb_path,cdtg , filename) ) + + if fextract is True and os.path.isfile(fpath): + os.remove ( fpath ) + elif os.path.isfile(fpath) and fextract is False and vrb in [1,2,3,4]: + print("ODB rows already in file for : obstype: {} , date: {}".format( obs, cdtg ) ) + + if os.path.isfile(fpath): + rows = self.io.ReadFrame( fpath ) + # Process rows from file + df_dist = self.rd.DfDist (rows, obs, cdtg, max_dist ) + + # Subset df + spl = SplitDf (df_dist , obs , cdtg ) + ndist = df_dist.dropna(subset=["dist"]).copy() + df_stat = spl.SubsetDf( bin_dist , max_dist) + self.dlist[obs].append( df_stat ) + + + else: + # Get (or Get again if fextracted =True) and process the rows from ODB + try: + rows= conn.odb_dict (cma_path , + sql_query , + nfunc , + float_fmt , + query_file, + poolmask , + pbar , + verbose ) + conn.odb_close() + if not rows: + if vrb in [0,1,2,3,4]: + print( "**WARNING : Data not available for obstype={} and varno={}".format( obs , varno[jo]) ) + empty_df= pd.DataFrame([]) + self.dlist[obs].append( empty_df ) + + except: + RuntimeError + if vrb in [0,1,2,3,4]: + print( "**WARNING : Data not available for obstype={} and varno={}".format( obs , varno[jo] )) + empty_df= pd.DataFrame([]) + self.dlist[obs].append( empty_df ) + pass + else: + # Process rows directly if the file not there + df_dist = self.rd.DfDist (rows, obs, cdtg, max_dist ) + + # Subset df + spl = SplitDf (df_dist , obs , cdtg ) + ndist = df_dist.dropna(subset=["dist"]).copy() + df_stat= spl.SubsetDf( bin_dist , max_dist) + self.dlist[obs].append( df_stat ) + + # Save data. + io_fpath = "/".join( ( self.odb_path, cdtg )) + rows_file = "_".join( ( "df_rows" , obs ,cdtg)) +".csv" + stat_file = "_".join( ( "df_stat" , obs ,cdtg)) +".csv" + + rows_fpath = "/".join( ( self.odb_path,cdtg , rows_file) ) + stat_fpath = "/".join( ( self.odb_path,cdtg , stat_file) ) + + # Write if not done ! + if not os.path.isfile(rows_fpath ): + df_rows = pd.DataFrame( rows ) + self.io.FlushFrame( df_rows , self.odb_path, cdtg , obs, fid=None ,verbose =vrb ,filetype="df_rows") + + if not os.path.isfile (stat_file): + self.io.FlushFrame( df_stat , self.odb_path, cdtg , obs, fid=None ,verbose =vrb ,filetype="df_stat") + return self.dlist + + + + + + + +# The methods below are ouside classes +# Since they may be used by multiple processes +def CreateSharedMfile (name, size): + # Delete and recreate the shared memory file if exists + # under /dev/shm + try: + # If exists close and delete + old = shared_memory.SharedMemory(name=name) # unlink and close + old.close() + old.unlink() + except FileNotFoundError: + pass + return shared_memory.SharedMemory(name=name, create=True, size=size) + + + + +def GcdistChunkShared (i0, i1, N): + """ + Method : Worker shared memory . + Only receives i0, i1 and global N. + Accesses shared memory numpy arrays directly. + returns: The matrix chunked blocks and begin/end index + """ + + # + CreateSharedMfile + shm_lon = shared_memory.SharedMemory(name="shm_lons") + shm_lat = shared_memory.SharedMemory(name="shm_lats") + + # Recreate numpy arrays as buffers + lon = np.ndarray((N,), dtype=np.float64, buffer=shm_lon.buf) + lat = np.ndarray((N,), dtype=np.float64, buffer=shm_lat.buf) + + # Get latlon by chunk indices + sub_lon = lon[i0:i1] + sub_lat = lat[i0:i1] + + # The method for computing the great circle distances is written in C + # The C code is more or less the same as the one used in R "sp" package + # --In order to get the same values as in the original obstool , the method + # in sp package was adapted to odb4py python package (The method is called odb_gcdist) + block = odb_gcdist (sub_lon, sub_lat, lon, lat) + + # Close memory buffers handle + shm_lon.close() + shm_lat.close() + return (i0, i1, block) + + + + +class DistMatrix: + """@Class : DistMatrix: computes the interdistances between all + the latlon pairs ( for N coordinates > 3000 , the process starts to slow down) + + Approach to speed up : Divide the coordinates arrays into chuncks and + compute each block in a sparated process + + @Returns: A reconstructed matrix of all distances between latlon pairs. + + Methods : + GcdistParallel + + """ + + def __init__(self, lons, lats): + self.lons = np.asarray(lons, dtype=np.float64) # Init with 64 bits + self.lats = np.asarray(lats, dtype=np.float64) + self.N = len(self.lons) + + def GcdistParallel(self, chunk_size=200, workers=None): + N = self.N + lon = self.lons + lat = self.lats + + # If the number of CPUs is set otherwise get from the machine + if workers is None: + workers = mp.cpu_count() + + # Create shared memory file + # More Safe + shm_lons =CreateSharedMfile ( size=lon.nbytes , name="shm_lons" ) + shm_lats =CreateSharedMfile ( size=lon.nbytes , name="shm_lats" ) + + # Copy data once + np.ndarray(lon.shape, dtype=np.float64, buffer=shm_lons.buf)[:] = lon + np.ndarray(lat.shape, dtype=np.float64, buffer=shm_lats.buf)[:] = lat + + # Chunk ranges + chunk_ranges = [(i0, min(i0 + chunk_size, N)) for i0 in range(0, N, chunk_size)] + + # Init the final matrix + distmat = np.zeros((N, N), dtype=np.float64) + try: + with mp.Pool(processes=workers) as pool: + # Use starmap method and give chunk ranges as args + results = pool.starmap( GcdistChunkShared , [(i0, i1, N) for (i0, i1) in chunk_ranges]) + + # Fill the matrix by block + for (i0, i1, block) in results: + distmat[i0:i1] = block + finally: + # Cleanup and unlink memory + shm_lons.close() + shm_lons.unlink() + shm_lats.close() + shm_lats.unlink() + + return distmat + + + +class Rows2Df : + + """ + Class: Build the dataframes with O-G and O-A departures. + To speed up the computation of the distances between the + latlon pairs , the matrix is computed by chunks + + Returns a DF of O-G , O-A , indices and distances + Methods : + DfDist + DfPrep + """ + + def __init__(self): + # I/O + self.io = DataIO () + return None + + + + def DfDist (self , rows , + obs , + cdtg , + max_dist , + verbosity =0 , + write_file=None, + file_path =None, + filetype =None + ): + #pd.set_option('display.max_rows', 20 ) + vrb=verbosity + if vrb not in [0,1,2,3,4]: + print("**WARNING : Min and max verbosity levels: 0 -> 4 Got : ", vrb ) + print("Fallback to default value: verbosity=", vrb ) + + if rows is None: + print("Rows from ODB not available for var : {}".format(var ) ) + sys.exit() + else: + df_rows = pd.DataFrame(rows ) + lats = np.array(rows["degrees(lat)" ]) + lons = np.array(rows["degrees(lon)" ]) + an_d = np.array(rows["an_depar@body"]) + fg_d = np.array(rows["fg_depar@body"]) + + # ARGS order : lons1, lats1, lon2, lat2 : Compute distances between latlon pairs + d = DistMatrix (lons , lats ) + + # The distances ar rounded to 0 decimals + # Else , there will be a number of observations sets in each bins. + # It leads to different statistics between R and Python whil perfomrming the cut , subset etc . + dist = d.GcdistParallel().round(0) + + N = len(lats) + dist_1d = dist.ravel() + + # Indices as in as.table(matdist) : row/col in R version + d1, d2 = np.indices((N, N)) + d1 = d1.ravel().astype(np.int32, copy=False) + d2 = d2.ravel().astype(np.int32, copy=False) + + # Invert d1, d2 to match the ones in R + ndist = pd.DataFrame({ + "n1" : d2, + "n2" : d1, + "dist": dist_1d }) + + # Subset using the max distance + df_dist = ndist[ndist["dist"] <= max_dist].copy() + + # Add OA/FG departures to corresponding indices + df_dist["OA1"] = an_d[df_dist["n1"].values] + df_dist["OA2"] = an_d[df_dist["n2"].values] + df_dist["FG1"] = fg_d[df_dist["n1"].values] + df_dist["FG2"] = fg_d[df_dist["n2"].values] + + # Set datatypes + df_dist = df_dist.astype({ + "n1": "int32", + "n2": "int32", + "dist": "int32", + "OA1": "float64", + "OA2": "float64", + "FG1": "float64", + "FG2": "float64" + }) + # For safety + df_dist = df_dist[df_dist["dist"] <= max_dist].copy() + + # Round to 6 decimal precision + return df_dist.round(6) + + + + def DfPrep (self , + frame_liste, + dir_=None , + period=None , + var_list=None ): + # Concat df by vars for the whole period + cnt= ConcatDf () + cdf= cnt.ConcatFromListe ( frame_liste) + return cdf diff --git a/scripts/CMakeLists.txt b/scripts/CMakeLists.txt index 7ce6714..91a656a 100644 --- a/scripts/CMakeLists.txt +++ b/scripts/CMakeLists.txt @@ -12,6 +12,8 @@ install_scripts (jbdiagnose.sh plotjbbal.py plotjbdiag.py) install_scripts (cv_header_list.sh) install_scripts (varbcdiagnose.sh plotvarbccoeff.py) install_scripts (rundiacov.sh diag_diacov.sh plotdiacov.py) -install_scripts(clddet_diagnosis.sh clddet_plotter.py) +install_scripts (clddet_diagnosis.sh clddet_plotter.py) install_scripts (datool_dfs.py) +install_scripts (tune_br.py) +install_scripts (obstool.py) #install_scripts (diacov.sh festat.sh) diff --git a/scripts/obstool.py b/scripts/obstool.py new file mode 100755 index 0000000..43d4774 --- /dev/null +++ b/scripts/obstool.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +import os, sys +import pandas as pd +from datetime import datetime +from pathlib import Path +import argparse + +# odb4py +from odb4py.utils import SqlParser , OdbObject + + +# Insert parent directory to get "modules" directory +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + +# Obstool modules +from modules.setting import Setting +from modules.utils import OdbReader , Rows2Df +from modules.conv_stats import DHLStat + + + +# Parse arguments from the command line +def Usage(): + help_= "Description :\n\ + Diagnostic tool 'obstool'.\n\ + Estimates the spatiale covariances , correlations and \n\ + standard deviations from the obs-guess and the obs-analysis\n\ + departures.\n\ + The output :\n\ + Estimated optimal thinning distance for a given observation\n\ + type variable.\n\ + \n\ + The implemented approaches are ones proposed by :\n\ + Hollingsworth/Lönnberg (1986) and Desroziers (2005)\n" + return help_ + +def ParseArgs(): + # The var list must be given as a liste sperated with ":" + var=lambda s:s.split(":") + ft =lambda x: x.lower() == "true" + parser = argparse.ArgumentParser( description=" ") + + obslist=["synop_t" ,"synop_z" ,"synop_u" ,"synop_v","synop_h", + "airep_t" ,"airep_u" ,"airep_v" , + "airepl_t","airepl_u","airepl_v", + "temp_t" ,"temp_u" ,"temp_v" ,"temp_q" , + "templ_u" ,"templ_u" ,"templ_v" ,"templ_q", + "dribu_t" ,"dribu_u" ,"dribu_v" ,"dribu_z", + "radar_rh","radar_dow" ] + + + # Args significations + arg1 ="Path to ODB directory. (default=.)" + arg2 ="ODB format/type. (default=CCMA)" + arg3 ="Begin date. Format YYYYMMDDHH. ( default=1970010100)" + arg4 ="End date. Format YYYYMMDDHH. ( default=1970010103)" + arg5 ="Observation category. (default=conv)" + arg6 ="List of observation types separated by ':' ( e.g: airep_t:synop_u)" + arg7 ="Cycle increment in hours. (default= 3)" + arg8 ="Maximum distance for diagnostics in [Km] (default=100 Km)" + arg9 ="Bin distance for diagnostics in [Km] (default=10 Km)" + arg10="Verbosity level. 0 to 4" + arg11="--force_extract ( without value ) :" \ + "The ODB rows are archived when obstool is running. " \ + "If the same period has been chosen the ODB rows " \ + "won't be extracted except if --force_extract is used in command line" + + # Required + parser.add_argument("--odb_path" ,required=True ,type=str ,default="." ,help=arg1) + parser.add_argument("--odb_type" ,required=True ,type=str ,default="CCMA" ,choices=["CCMA"], help=arg2 ) + parser.add_argument("--bdate" ,required=True ,type=str ,default="1970010100",help=arg3) + parser.add_argument("--edate" ,required=True ,type=str ,default="1970010103",help=arg4) + parser.add_argument("--obs_category" ,required=True ,type=str ,default="conv" ,choices=["conv","satem"], help=arg5) + parser.add_argument("--var_list" ,required=True ,type=var ,default="synop_t" ,choices=None ,help=arg6) + parser.add_argument("--cycle_inc" ,required=True ,type=int ,default=3 ,choices=None ,help=arg7) + + # Optional + parser.add_argument("--max_dist" ,required=False ,type=float ,default=100.0 ,choices=None ,help=arg8) + parser.add_argument("--bin_dist" ,required=False ,type=float ,default=10.0 ,choices=None ,help=arg9) + parser.add_argument("--verbose" ,required=False ,type=int ,default=1 ,choices=None ,help=arg10) + parser.add_argument("--force_extract",action="store_true", help=arg11 ) + # + if len(sys.argv) == 1: + Usage() + parser.print_help() + sys.exit(1) + return parser.parse_args() + + +# Start +start_time = datetime.now() + +# Get args +args= ParseArgs () + +# Odb +odbpath= args.odb_path +odbtype= args.odb_type + +# Period and cycle +bdate = args.bdate +edate = args.edate +cycle_inc=args.cycle_inc + +# Set variable list +var_list= [] +for v in args.var_list: var_list.append( v.lstrip () ) + +# Extract ODB rows again +fextract=args.force_extract + +# verbose +verb=args.verbose + +# Max and binning distances +max_dist= args.max_dist +bin_dist= args.bin_dist + + + + +# Instantiate +st = Setting () +rd = Rows2Df () +rr = OdbReader(odbpath , odbtype ) + +# Set period +period=st.set_period( bdate, edate ) + +# Set and check the var list +st.set_obs_list( var_list ) + +# Collection of pre-selected frames with distances <= max_dist +frame_liste = rr.get_odb_rows (period , + var_list , + max_dist , + bin_dist , + odbpath , + fextract , + cycle_inc , + pbar =True , + verbosity = verb ) + +# Concat Df for the final stats +cdf = rd.DfPrep( frame_liste ) + +# Time duration of ODB extraction +ext_time = datetime.now() + +# Diff +extraction_time = ext_time - start_time +print(f"Extraction time : {extraction_time}") + +# Get the stats by var +for var in var_list: + dhl= DHLStat ( cdf[var] , max_dist = max_dist ,bin_dist = bin_dist ) + cov= dhl.getCov( var , inplace=False ) + sig= dhl.getSig( var , inplace=False ) + cor= dhl.getCor( var , inplace=False ) + diag=dhl.getStatFrame(var) + # Save stats in a csv file + diag.to_csv ( "py_"+var+".csv", index=False , header=True ) + + +# Finish +end_time = datetime.now() + +# Time duration of stats computation +duration = end_time - ext_time +print(f"Statistics calculation duration : {duration}") + +# Total +tot_duration = end_time - start_time +print(f"Total duration runtime : {tot_duration}") +quit() diff --git a/scripts/tune_br.py b/scripts/tune_br.py new file mode 100755 index 0000000..fbce5f4 --- /dev/null +++ b/scripts/tune_br.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os , sys +import configparser +import argparse +from datetime import datetime ,timedelta +from statistics import mean +from pathlib import Path + +# TuneBR modules +# Insert parent directory to get "modules" directory +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + +from modules.sigma_bo import Predef, Diag , Ratios +from modules.read_odb import Odb ,TuneEnv +from modules.gsacov import GSA + + + + +# Get the config file +parser = argparse.ArgumentParser() +parser.add_argument( + "config_file", + type=Path, + help="Path to tuneBR config file ( .ini format )" +) + +args = parser.parse_args() + +if not args.config_file.is_file(): + parser.error(f"Config file '{args.config_file}' not found") + + +ini_file=args.config_file + +# Start +StartTime = datetime.now() + + +# Parse config file +config=configparser.ConfigParser( +interpolation=configparser.ExtendedInterpolation() ) + +# All items in upper case +config.optionxform = str +config.read(ini_file) + +# Init TuneBR env +env = TuneEnv ( config ) +PathDict, ModelDict = env.__Dicts__() +bdate =env.BeginDate +edate =env.EndDate +cycle_inc=env.cycle_inc + + +# Init GSA fortran routine +statfile =env.stabal +nlev =env.nflev +nsmax =env.nsmax +deltax =env.deltax +lverb =env.llverb +lwrite =env.lwrite +outfile =env.outfile +prog_bar =env.prog_bar +ndca_cpu =env.dca_cpu +rows_path=env.rows_path + + +# Get background standard deviations (Profiles + means ) +g=GSA(PathDict , statfile ,nsmax ,nlev , deltax , lverb ,lwrite ) +tsig_ver , sb_pred_t = g.GetSigmaB (2) # Temperature KPAR=2 +qsig_ver , sb_pred_q = g.GetSigmaB (3) # Specific q KPAR=3 +vsig_ver , sb_pred_v = g.GetSigmaB (4) # Vorticity KPAR=4 +dsig_ver , sb_pred_d = g.GetSigmaB (5) # Divergence KPAR=5 +kesig_ver, sb_pred_ke = g.GetSigmaB (999) # UV Components not in stabal file --> set arbitrary unique number 999 +print("Extraction of SIGMA_B values, done !\n") + + +# Create datetime list +cdtg=[] +bdate =datetime.strptime( bdate , "%Y%m%d%H") +edate =datetime.strptime( edate , "%Y%m%d%H") +delta =timedelta(hours=int(cycle_inc)) +while bdate <= edate: + strdate=bdate.strftime("%Y%m%d%H") + cdtg.append( strdate ) + bdate += delta + + +# ODB extraction +db=Odb ( PathDict ) + +# Get ODB rows +db.CreateDca ( cdtg , ndca_cpu ) +db.OdbExtract ( cdtg ,rows_path , prog_bar ) + +# Predefined SIGMA_O +print("Compute predefined SIGMA_O ...!\n") +so_pred_t =Predef ( PathDict , cdtg , lverb , lwrite).GetSigmaP ("t" ) +so_pred_bt=Predef ( PathDict , cdtg , lverb , lwrite).GetSigmaP ("bt") +so_pred_q =Predef ( PathDict , cdtg , lverb , lwrite).GetSigmaP ("q" ) +so_pred_ke=Predef ( PathDict , cdtg , lverb , lwrite).GetSigmaP ("ke") + + +# Compute SIGMA_O AND SIGMA_B diagnostics +print("Compute SIGMA_O, SIGMA_B diagnostics ...!\n") +sb_diag_t , so_diag_t , pt =Diag(PathDict ,cdtg, lverb , lwrite).GetSigmaD("t" ) +sb_diag_bt , so_diag_bt, pbt=Diag(PathDict ,cdtg, lverb , lwrite).GetSigmaD("bt") +sb_diag_q , so_diag_q , pq =Diag(PathDict ,cdtg, lverb , lwrite).GetSigmaD("q" ) +sb_diag_ke , so_diag_ke, pke=Diag(PathDict ,cdtg, lverb , lwrite).GetSigmaD("ke") + + +# Use the same notation as in the RC-LACE version +sb_pred=[sb_pred_t,sb_pred_q,sb_pred_ke] # PREDEFINED Sb (sb_bt PREDEFINED DOESN'T EXIST FOR BRIGHTNESS T) +so_pred=[so_pred_t,so_pred_bt,so_pred_q,so_pred_ke] # // So +sb_diag=[sb_diag_t,sb_diag_bt,sb_diag_q,sb_diag_ke] # DIAGNOSED Sb +so_diag=[so_diag_t,so_diag_bt,so_diag_q,so_diag_ke] # // So + +# Total wind obs is : nobs ke/2 +nobs =[ pt , pbt , pq , pke/2. ] + +# OBS Mean +nobs_mean =int(sum(nobs)/len(nobs)) + +# INIT RATIO OBJECT WITH CORRESPONDING PREDEF AND DIAG LISTS +rednmc=env.rednmc +r=Ratios(PathDict, nobs, rednmc , so_pred , so_diag , sb_pred , sb_diag ,lwrite) + +# GET RATIOS +sigo = r.RatioSo() +sigb = r.RatioSb() + +rot , robt , roq , roke,roav = sigo[0],sigo[1],sigo[2],sigo[3],sigo[4] +rbt , rbq , rbke , rbav = sigb[0],sigb[1],sigb[2],sigb[3] + +# Output +text = ( + 60*"-" + "\n" + + "Var | cases | Ratio_o | Ratio_b".center(50, " ") + "| \n" + + 60*"-" + "\n" + + f"t |{str(pt).center(15)}{str(round(rot,5)).center(15)}{str(round(rbt,5)).center(15)} \n" + + f"bt |{str(pbt).center(15)}{str(round(robt,5)).center(15)}{'None'.center(15)} \n" + + f"q |{str(pq).center(15)}{str(round(roq,5)).center(15)}{str(round(rbq,5)).center(15)} \n" + + f"ke |{str(pke).center(15)}{str(round(roke,5)).center(15)}{str(round(rbke,5)).center(15)} \n" + + 60*"-" + "\n" + + f"Mean |{str(nobs_mean).center(15)}{str(round(roav,5)).center(15)}{str(round(rbav,5)).center(15)}\n" + + 60*"-" +) + +# Simple print on screen +print( text) + +# Write in file given in config + date1_date2 +if lwrite == True : + period = str(cdtg[0]) +"_"+ str(cdtg[-1] ) + outfile= "_".join( ( outfile, period)) + dir_ = os.getcwd() + print( f"SIGMA B/O written in file: {dir_}/{outfile} " ) + with open( outfile , "w", encoding="utf-8") as f: + f.write(text) + +# +if lverb == True: + print("\n"+"The different I/O files are written in "+os.getenv("PWD")+"/out") + print( " ") + EndTime = datetime.now() + Duration=EndTime - StartTime + print( " EXECUTION TIME : \n" , Duration ) +# END +quit() + diff --git a/share/CMakeLists.txt b/share/CMakeLists.txt index 5c8bb18..4ca9fbb 100644 --- a/share/CMakeLists.txt +++ b/share/CMakeLists.txt @@ -1,2 +1,2 @@ -install( DIRECTORY levdef DESTINATION share PATTERN "CMakeLists.txt" EXCLUDE ) - +install( DIRECTORY levdef DESTINATION share PATTERN "CMakeLists.txt" EXCLUDE ) +install( DIRECTORY tunebr_conf DESTINATION share PATTERN "CMakeLists.txt" EXCLUDE ) diff --git a/share/tunebr_conf/config.ini b/share/tunebr_conf/config.ini new file mode 100755 index 0000000..8a07bd0 --- /dev/null +++ b/share/tunebr_conf/config.ini @@ -0,0 +1,48 @@ +# Start and end dates +[DATES] +DATESTART=2024011003 +DATEEND =2024011009 + +[PATHS] +# Basedir +BASEDIR=/scratch/cvah/tunebr + +# ODB path +ODBPATH=${BASEDIR}/odb + +# Workdir +WORKDIR=${BASEDIR}/tmp + +# Temporarly extracted ODB files +ROWS_PATH=${BASEDIR}/odb_rows + +# B matrix file +STATFILE =/hpcperm/cvah/tuneBR_1.1.0/stab/stabcvt + +# The filename where the CCMA is assumed archived in +ODB_TEMPLATE=CCMA_minim_YYYYMMDDHH.tar + + +# Some model setting +[MODEL] +CYCLE_INC=3 +DELTAX =1.3 +NSMAX =287 +NFLEV =87 +REDNMC =1.0 + +[OPTIONS] +# Number of CPUs for DCA files creation +DCA_CPU= 16 +# Progress bar when the ODB rows are fetched ! +PROGBAR=True + +# Sigma B/O will be written in file +OUTFILE=sigmabo_ratios + +# Enable verbosity +LLVERB=False + +# Enable file writing +# The output files of the diffrent steps will be written in $BASEDIR/out +LWRITE=True diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..1fe4a69 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,70 @@ +# CmakeFile to compile readgsa.F90 +# Creates a virtualenv , where f2py compiles the routine +# When the venv is created , the numpy module is installed. +# ==> Eveything is done inside the build directory ! +cmake_minimum_required (VERSION 3.18) + +project (readgsa LANGUAGES C Fortran) + + +# Find system python +find_package(Python REQUIRED COMPONENTS Interpreter) + +message(STATUS "System python: ${Python_EXECUTABLE}") + +# Create venv inside build +set(VENV_DIR ${CMAKE_BINARY_DIR}/venv) +if(NOT EXISTS ${VENV_DIR}/bin/python) + message(STATUS "Creating virtualenv ...") + execute_process( + COMMAND ${Python_EXECUTABLE} -m venv ${VENV_DIR} + RESULT_VARIABLE venv_status + ) + if(NOT venv_status EQUAL 0) + message(FATAL_ERROR "Failed to create virtualenv") + endif() + execute_process( COMMAND ${VENV_DIR}/bin/pip install --upgrade pip ) + execute_process( COMMAND ${VENV_DIR}/bin/pip install numpy ) +endif() + + +# Switch python to venv +set(Python_EXECUTABLE ${VENV_DIR}/bin/python CACHE FILEPATH "" FORCE) + +# Look for numpy inside venv +find_package(Python REQUIRED COMPONENTS Interpreter Development NumPy) +message(STATUS "Python from venv: ${Python_EXECUTABLE}") + +# Locate f2py includes +execute_process( COMMAND ${Python_EXECUTABLE} -c "import numpy.f2py;print(numpy.f2py.get_include())" OUTPUT_VARIABLE F2PY_INCLUDE_DIR OUTPUT_STRIP_TRAILING_WHITESPACE) + +# Common variables +set(f2py_module_name "readgsa") +set(fsrc_path "${CMAKE_SOURCE_DIR}/src/readgsa.F90") +set(fortran_src_file "${fsrc_path}" ) +set(f2py_module_c "${f2py_module_name}module.c") + +add_custom_command( OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/${f2py_module_c}" COMMAND ${Python_EXECUTABLE} -m "numpy.f2py" "${fortran_src_file}" + -m ${f2py_module_name} + --lower + DEPENDS readgsa.F90 # Fortran source +) + +# Set up target +Python_add_library( readgsa MODULE WITH_SOABI + "${CMAKE_CURRENT_BINARY_DIR}/${f2py_module_c}" # Generated .c module + "${F2PY_INCLUDE_DIR}/fortranobject.c" # From NumPy + "${fortran_src_file}" # Fortran source +) + +# Include +target_include_directories(${f2py_module_name} PRIVATE ${F2PY_INCLUDE_DIR} ${Python_NumPy_INCLUDE_DIRS}) + +# Libs +target_link_libraries(${f2py_module_name} PRIVATE Python::Module Python::NumPy) + +# Set the suffix to " .so " without machine_system architecture/realease info ..etc +set_target_properties(readgsa PROPERTIES PREFIX "" SUFFIX ".so") + +# Add it to the python modules +install(TARGETS readgsa LIBRARY DESTINATION ${CMAKE_SOURCE_DIR}/modules ) diff --git a/src/readgsa.F90 b/src/readgsa.F90 new file mode 100644 index 0000000..26a00eb --- /dev/null +++ b/src/readgsa.F90 @@ -0,0 +1,122 @@ + SUBROUTINE READCOV(pcov,kdim,nsmax,kpar,cfile ,llverb,kret) + IMPLICIT NONE +! read a set of nonsep matrices for the requested parameter +! from a GSA file +! F. Bouttier Oct 96 +! input : cfile = GSAfilename +! kpar = GSA parameter selector (IPAR coding) +! output : pcov = covariance matrices +! kret = ok if 0, not found if 1 + +! ROUTINE BASED ON +! THE MAIN PROGRAM stdav.F90 : G. Boloni based on "fediacov" 2010/09/09 +! : B. Strajnar 2010/10/26 +! USE AS PYTHON EXTENSION : I. Dehmous 2023/10/03 + + + INTEGER, PARAMETER :: jplun=67 + INTEGER, PARAMETER :: JPRB = SELECTED_REAL_KIND(13,300) + INTEGER, PARAMETER :: JPIM = SELECTED_INT_KIND(9) + INTEGER(KIND=JPIM),INTENT(IN ) :: kdim,nsmax ,kpar + INTEGER(KIND=JPIM),INTENT(OUT) :: kret + INTEGER(KIND=JPIM) :: idim1,idim2,ipar1,ipar2,icor1,icor2, icheck + REAL(KIND=JPRB) ,INTENT(OUT) :: pcov(kdim,kdim,0:nsmax) + CHARACTER(LEN=*) :: cfile + CHARACTER(LEN=8) :: clid + CHARACTER(LEN=70) :: clcom + LOGICAL ,INTENT(IN) :: llverb + LOGICAL :: llfound + INTEGER(KIND=JPIM) :: iorig,idate,itime,inbset + INTEGER(KIND=JPIM) :: jset , jj , jn , jk + INTEGER(KIND=JPIM) :: inbmat,iweight,itypmat,isetdist,ilendef + + +! 1.OPEN FILE AND READ GSA HEADER +! THE FILE IS WRITTEN IN BIG ENDIAN MODE +! This option is added also in compiler options anyway + + OPEN(jplun,FILE=cfile,FORM='unformatted',STATUS='old',& + & ACCESS='sequential' ,convert='BIG_ENDIAN') + + IF (llverb) THEN + WRITE(*,*) 'Searching for parameter :',kpar, & + & 'in GSA file: ',cfile,'...' + ENDIF + + READ(jplun) clid + READ(jplun) clcom +! This read statement below is not always valid!!! +! In some versions of festat it is not written to the stabal file!!! + READ(jplun) iorig,idate,itime,inbset + IF (llverb) THEN + WRITE(*,*) 'ID =',clid + WRITE(*,*) 'COMMENT=',clcom + WRITE(*,*) 'ORIGINE,DATE,TIME,INBSET',iorig,idate,itime,inbset + ENDIF + +! 2. SEARCH THE SET + llfound=.false. + DO jset=1,inbset + READ(jplun) inbmat,iweight,itypmat,isetdist,ilendef + !IF (llverb) THEN + WRITE(*,*) 'Scanning set No...',jset + !ENDIF + IF (llverb) THEN + WRITE(*,*) ' ',inbmat,iweight,itypmat,isetdist,ilendef + ENDIF + READ(jplun) idim1,idim2,ipar1,ipar2,icor1,icor2 + IF (llverb) THEN + WRITE(*,*) ' ',idim1,idim2,ipar1,ipar2,icor1,icor2 + ENDIF + IF ((ipar1==kpar).AND.(ipar2==kpar)) llfound=.true. + +! skip coordinate records + READ(jplun) + READ(jplun) + IF (llfound) THEN + IF (llverb ) THEN + WRITE(*,*) 'Found param :', kpar ,' In the set No :', jset + WRITE(*,*) 'Reading data ...' + ENDIF + +! THE OLD ABORT ROUTINE CAUSES PYTHON CRASH + IF (inbmat/=nsmax+1)call abor1('Error - bad matrix number') + IF (itypmat/=1 )call abor1('Error - not vert cov matrices') + IF (isetdist/=2 )call abor1('Error - not distrib by wave n') + IF (ilendef/=1 )call abor1('Error - bad def list length') + IF (idim1/=kdim )call abor1('Error - bad 1st dimension') + IF (idim2/=kdim )call abor1('Error - bad 2nd dimension') + IF (icor1/=icor2 )call abor1('Error - coords differ') + + DO jn=0,nsmax + READ(jplun) + READ(jplun) ((pcov(jj,jk,jn),jk=1,kdim),jj=1,kdim),icheck + IF (icheck/=3141592) call abor1('Error - wrong check word') + ENDDO + EXIT + ELSE +! SKIP THE SET + DO jn=1,inbmat + IF(ilendef>0) read(jplun) + READ(jplun) + ENDDO + ENDIF + ENDDO + CLOSE(jplun) + IF (llfound) THEN + kret=0 + ELSE + WRITE(*,*) 'WARNING : Could not find matrix set' + pcov(:,:,:)=0. + kret=1 + ENDIF + RETURN + END SUBROUTINE READCOV + +!-------------------------abor1 ------------------------- + SUBROUTINE abor1(cdmess) + character(len=*) :: cdmess + WRITE(*,*) cdmess + !call exit(1) + END +! Finish diff --git a/tests/tools/CMakeLists.txt b/tests/tools/CMakeLists.txt index 212890c..400ddf0 100644 --- a/tests/tools/CMakeLists.txt +++ b/tests/tools/CMakeLists.txt @@ -8,6 +8,8 @@ list( APPEND tools_tests_scripts test_plotdiacov_help.sh test_varbcdiagnose_help.sh test_plotvarbcpred_help.sh + test_tune_br_help.sh + test_obstool_help.sh ) diff --git a/tests/tools/test_obstool_help.sh b/tests/tools/test_obstool_help.sh new file mode 100755 index 0000000..69fa20a --- /dev/null +++ b/tests/tools/test_obstool_help.sh @@ -0,0 +1,17 @@ +#!/usr/bin/env bash + +set -uex + +# A unique working directory +wd=$(pwd) +test_wd=$(pwd)/test_obstool_help + +mkdir -p ${test_wd} +cd ${test_wd} + +obstool.py -h > /dev/null && echo "obstool help/usage works" + + +# Clean up +cd ${wd} +rm -rf ${test_wd} diff --git a/tests/tools/test_tune_br_help.sh b/tests/tools/test_tune_br_help.sh new file mode 100755 index 0000000..7de9bb2 --- /dev/null +++ b/tests/tools/test_tune_br_help.sh @@ -0,0 +1,17 @@ +#!/usr/bin/env bash + +set -uex + +# A unique working directory +wd=$(pwd) +test_wd=$(pwd)/test_tune_br_help + +mkdir -p ${test_wd} +cd ${test_wd} + +tune_br.py -h > /dev/null && echo "tune_br help/usage works" + + +# Clean up +cd ${wd} +rm -rf ${test_wd}