diff --git a/src/bitsea/validation/online/RunSatValidation.py b/src/bitsea/validation/online/RunSatValidation.py index 4bdeb381..45593f5b 100644 --- a/src/bitsea/validation/online/RunSatValidation.py +++ b/src/bitsea/validation/online/RunSatValidation.py @@ -1,5 +1,5 @@ import argparse - +# For MER def argument(): parser = argparse.ArgumentParser(description = 'run sat validation on N procs') parser.add_argument( '--inputdir', '-i', diff --git a/src/bitsea/validation/online/RunSatValidation_MER.py b/src/bitsea/validation/online/RunSatValidation_MER.py new file mode 100644 index 00000000..b3053e39 --- /dev/null +++ b/src/bitsea/validation/online/RunSatValidation_MER.py @@ -0,0 +1,115 @@ +import argparse +# For MER +def argument(): + parser = argparse.ArgumentParser(description = 'run sat validation on N procs') + parser.add_argument( '--inputdir', '-i', + type = str, + required = True, + help = 'Archive of P_l from model') + + parser.add_argument( '--outputdir', '-o', + type = str, + required = True, + help = 'Output dir of sat validation') + + + parser.add_argument( '--satdir', '-s', + type = str, + required = True, + help = 'Daily satellite dir') + + parser.add_argument( '--climdir', '-c', + type = str, + required = True, + help = 'Clim dir') + + parser.add_argument( '--maskfile', '-m', + type = str, + required = True, + help = 'Maskfile') + parser.add_argument( '--var', '-v', + type = str, + required = True, + choices = ['P_l','kd490']) + + parser.add_argument( '--rundate', '-r', + type = str, + required = True, + help = 'OPA_RUNDATE') + + return parser.parse_args() + +args = argument() + +from datetime import datetime,timedelta + +from bitsea.commons.mask import Mask +from bitsea.commons.submask import SubMask +from bitsea.basins import V2 as OGS +from bitsea.commons.utils import addsep +import numpy as np +from SatValidation import SatValidation + +try: + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nranks =comm.size +except: + rank = 0 + nranks = 1 + + + +SAT_DAILY_DIR = addsep(args.satdir) +OUTDIR = addsep(args.outputdir) +CLIM_DIR = addsep(args.climdir) +ARCHIVEDIR = addsep(args.inputdir) +OPA_RUNDATE = args.rundate + + +TheMask=Mask(args.maskfile) + +#basins +#nSUB = len(OGS.P.basin_list) +nSUB = len(OGS.adr.basin_list) +jpk,jpj,jpi =TheMask.shape +mask20_2D = TheMask.mask_at_level(20.0) +dtype = [(sub.name, bool) for sub in OGS.P] +SUB = np.zeros((jpj,jpi),dtype=dtype) + +for sub in OGS.Pred: + SUB[sub.name] = SubMask(sub,maskobject=TheMask).mask_at_level(0) + if 'atl' in sub.name: continue + SUB['med'] = SUB['med'] | SUB[sub.name] + +COASTNESS_LIST=['coast','open_sea','everywhere'] + +dtype = [(coast, bool) for coast in COASTNESS_LIST] +COASTNESS = np.ones((jpj,jpi),dtype=dtype) +COASTNESS['coast'] = ~mask20_2D +COASTNESS['open_sea'] = mask20_2D + +SUFFIX={'P_l' : '_cmems_obs-oc_med_bgc-plankton_nrt_l3-multi-1km_P1D.nc', + 'votemper': '_cmems_obs-oc_med_bgc-transp_nrt_l3-multi-1km_P1D.nc' #? + } # for votemper set the correct SAT_file + +def climfilename(CLIM_DIR, month,var): + if var=='P_l': + return CLIM_DIR + '/ave.yyyy' + month + '15-00:00:00.nc' + else: + return None + +delay=-2 +opa_rundate=datetime.strptime(OPA_RUNDATE,'%Y%m%d') +FC_RUNDATE = opa_rundate + timedelta(days=delay) # FC_RUNDATE is the forecast at time 0, of 2 days ago +daydate = FC_RUNDATE +day = daydate.strftime('%Y%m%d') +month = daydate.strftime('%m') +avefile = "ave.%s-12:00:00.%s.nc" %(day, args.var) +LOCAL_FC = "%s%s/%s" %(ARCHIVEDIR, FC_RUNDATE.strftime('%Y%m%d'),avefile) +SAT_FILE = SAT_DAILY_DIR + day + SUFFIX[args.var] +CLIM_FILE= climfilename(CLIM_DIR, month, args.var) +f_name = OUTDIR + 'Validation_f' + str(fc_day+1) + '_' + FC_RUNDATE.strftime('%Y%m%d') + '_on_daily_Sat.' + day + '.nc' +SatValidation(args.var, LOCAL_FC,SAT_FILE,CLIM_FILE,TheMask,f_name,SUB,COASTNESS_LIST,COASTNESS,nSUB) + diff --git a/src/bitsea/validation/online/SatValidation_MER.py b/src/bitsea/validation/online/SatValidation_MER.py new file mode 100644 index 00000000..eb730256 --- /dev/null +++ b/src/bitsea/validation/online/SatValidation_MER.py @@ -0,0 +1,142 @@ +import numpy as np +import bitsea.Sat.SatManager as Sat +import bitsea.matchup.matchup as matchup +from bitsea.commons.dataextractor import DataExtractor +from bitsea.basins import V2 as OGS #MODIFICARE CON I SOTTOBACINI del NAd (al momento รจ adr1) +import scipy.io.netcdf as NC +import os + +def weighted_mean(Conc, Weight): + + Weight_sum = Weight.sum() + Mass = (Conc * Weight).sum() + Weighted_Mean = Mass/Weight_sum + return Weighted_Mean + +def weighted_var(Conc, Weight): + ConcMean = weighted_mean(Conc,Weight) + Weight_sum = Weight.sum() + Mass = ((Conc-ConcMean)**2 * Weight).sum() + Weighted_Var = Mass/Weight_sum + return Weighted_Var + +def SatValidation(var, modfile,satfile,climafile,TheMask,outfile,SUB,COASTNESS_LIST,COASTNESS,nSUB): + assert var in ['P_l','kd490'] + varname_sat={'P_l':'CHL','kd490':'KD490'} + if not os.path.exists(satfile): + print("Not existing file: " + satfile) + return + + print(outfile) + nCOAST = len(COASTNESS_LIST) + De = DataExtractor(TheMask,filename=modfile,varname=var,dimvar=2) + Model = De.values[:,:] + + if climafile is not None: + Declim = DataExtractor(TheMask,filename=climafile,varname=var,dimvar=2) + Clima = Declim.values[:,:] +# try: +# Sat24 = Sat.readfromfile(satfile,varname_sat[var]) # weekly +# except: + Sat24 = Sat.convertinV4format(Sat.readfromfile(satfile, varname_sat[var])) # daily + + cloudsLand = np.isnan(Sat24) + Sat24[cloudsLand] = -999.0 + cloudsLand = Sat24==-999.0 + modelLand = Model > 1.0e+19 + nodata = cloudsLand | modelLand + + + + VALID_POINTS = np.zeros((nSUB,nCOAST),np.float32) + + MODEL_MEAN = np.zeros((nSUB,nCOAST),np.float32) + SAT___MEAN = np.zeros((nSUB,nCOAST),np.float32) + MODEL_VARIANCE = np.zeros((nSUB,nCOAST),np.float32) + SAT___VARIANCE = np.zeros((nSUB,nCOAST),np.float32) + BGC_CLASS4_CHL_BIAS_SURF_BASIN = np.zeros((nSUB,nCOAST),np.float32) + BGC_CLASS4_CHL_RMS_SURF_BASIN = np.zeros((nSUB,nCOAST),np.float32) + BGC_CLASS4_CHL_CORR_SURF_BASIN = np.zeros((nSUB,nCOAST),np.float32) + + + ANOMALY_CORR = np.zeros((nSUB,nCOAST),np.float32) + + for icoast, coast in enumerate(COASTNESS_LIST): + +# for isub, sub in enumerate(OGS.P): + for isub, sub in enumerate(OGS.adr): + selection = SUB[sub.name] & (~nodata) & COASTNESS[coast] + M = matchup.matchup(Model[selection], Sat24[selection]) + VALID_POINTS[isub,icoast] = M.number() + if M.number() > 0 : + BGC_CLASS4_CHL_RMS_SURF_BASIN[isub,icoast] = M.RMSE() + BGC_CLASS4_CHL_BIAS_SURF_BASIN[isub,icoast] = M.bias() + BGC_CLASS4_CHL_CORR_SURF_BASIN[isub,icoast] = M.correlation() + + weight = TheMask.area[selection] + MODEL_MEAN[isub,icoast] = weighted_mean( M.Model,weight) + SAT___MEAN[isub,icoast] = weighted_mean( M.Ref, weight) + MODEL_VARIANCE[isub,icoast] = weighted_var( M.Model,weight) + SAT___VARIANCE[isub,icoast] = weighted_var( M.Ref, weight) + + + + if climafile is not None: + Mclima = matchup.matchup(Model[selection]-Clima[selection], \ + Sat24[selection]-Clima[selection]) + ANOMALY_CORR[isub,icoast] = Mclima.correlation() + + else: + BGC_CLASS4_CHL_RMS_SURF_BASIN[isub,icoast] = np.nan + BGC_CLASS4_CHL_BIAS_SURF_BASIN[isub,icoast]= np.nan + BGC_CLASS4_CHL_CORR_SURF_BASIN[isub,icoast]= np.nan + MODEL_MEAN[isub,icoast] = np.nan + SAT___MEAN[isub,icoast] = np.nan + MODEL_VARIANCE[isub,icoast] = np.nan + SAT___VARIANCE[isub,icoast] = np.nan + + ANOMALY_CORR[isub,icoast] = np.nan + + + ncOUT = NC.netcdf_file(outfile,'w') + ncOUT.createDimension('nsub', nSUB) + ncOUT.createDimension('ncoast', nCOAST) + s='' + for sub in OGS.P: s =s+sub.name + "," + setattr(ncOUT,'sublist',s) + s='' + for coast in COASTNESS_LIST: s =s+coast + "," + setattr(ncOUT,'coastlist',s) + + + ncvar= ncOUT.createVariable('number', 'i', ('nsub','ncoast')) + ncvar[:]= VALID_POINTS + + ncvar = ncOUT.createVariable('BGC_CLASS4_CHL_RMS_SURF_BASIN', 'f', ('nsub','ncoast')) + ncvar[:] = BGC_CLASS4_CHL_RMS_SURF_BASIN + + ncvar = ncOUT.createVariable('BGC_CLASS4_CHL_BIAS_SURF_BASIN', 'f', ('nsub','ncoast')) + ncvar[:] = BGC_CLASS4_CHL_BIAS_SURF_BASIN + + ncvar = ncOUT.createVariable('BGC_CLASS4_CHL_CORR_SURF_BASIN', 'f', ('nsub','ncoast')) + ncvar[:] = BGC_CLASS4_CHL_CORR_SURF_BASIN + + ncvar = ncOUT.createVariable('MODEL_MEAN', 'f', ('nsub','ncoast')) + ncvar[:] = MODEL_MEAN + + ncvar=ncOUT.createVariable('SAT___MEAN', 'f', ('nsub','ncoast')) + ncvar[:] = SAT___MEAN + + ncvar = ncOUT.createVariable('MODEL_VARIANCE', 'f', ('nsub','ncoast')) + ncvar[:] = MODEL_VARIANCE + + ncvar=ncOUT.createVariable('SAT___VARIANCE', 'f', ('nsub','ncoast')) + ncvar[:] = SAT___VARIANCE + + + if climafile is not None: + ncvar=ncOUT.createVariable('ANOMALY_CORRELATION', 'f', ('nsub','ncoast')) + ncvar[:] = ANOMALY_CORR + + ncOUT.close() + diff --git a/src/bitsea/validation/online/Validation_C3_MER.sh b/src/bitsea/validation/online/Validation_C3_MER.sh new file mode 100644 index 00000000..e6b50949 --- /dev/null +++ b/src/bitsea/validation/online/Validation_C3_MER.sh @@ -0,0 +1,42 @@ +# SST: + +cd $OPA_BITSEA/validation/online +# SAT DIR for Temperature: +SAT_DAILY_DIR=${ONLINE_REPO}/SAT/SST/NRT/DAILY/CHECKED_24 +# Model dir: ave_files +inputdir=${ONLINE_VALIDATION_DIR}/ARCHIVE_FC_FOR_SAT/ +# Output dir: name of the day to validate +outputdir=${ONLINE_VALIDATION_DIR}/SST + +# OPA_RUNDATE is the actual date in which I am running the model +opa_prex_or_die "python RunSatValidation_MER.py -i $inputdir -o $outputdir -s $SAT_DAILY_DIR -c null -m $MASKFILE -r $OPA_RUNDATE -v votemper" +INPDIR_SAT_V_DAILY=${OPA_INPDIR_ROOT}/analysis/VALIDATION/SAT/SST/NRT/ +opa_cp "$outputdir/Validation_*nc $INPDIR_SAT_V_DAILY" + +################## +#CHL: + +SAT_DAILY_DIR=${ONLINE_REPO}/SAT/CHL/NRT/DAILY/CHECKED_24 + +inputdir=${ONLINE_VALIDATION_DIR}/ARCHIVE_FC_FOR_SAT/ +outputdir=${ONLINE_VALIDATION_DIR}/CHL +# VERIFICARE il CLIMDIR: temporary we set null +#climdir=${ONLINE_VALIDATION_DIR}/CHL_CLIM/2D/ + +#opa_prex_or_die "python RunSatValidation.py -i $inputdir -o $outputdir -s $SAT_DAILY_DIR -c $climdir -m $MASKFILE -r $OPA_RUNDATE -v P_l" +opa_prex_or_die "python RunSatValidation_MER.py -i $inputdir -o $outputdir -s $SAT_DAILY_DIR -c null -m MASKFILE -r $OPA_RUNDATE -v P_l" + + +INPDIR_SAT_V_DAILY=${OPA_INPDIR_ROOT}/analysis/VALIDATION/SAT/CHL/NRT/ +opa_cp "$outputdir/Validation_*nc $INPDIR_SAT_V_DAILY" +opa_prex "python sat_plotter.py -i ${INPDIR_SAT_V_DAILY} -o ${outputdir} -d $OPA_RUNDATE -v P_l" +opa_cp " ${outputdir}/*png $ONLINE_VALIDATION_DIR/MEDEAF_PUB/Sat_validation" + + +# MAP CREATION, for SST and CHLa: +# for the SATmap, use as basis the script: +# https://github.com/inogs/bit.sea/blob/master/src/bitsea/validation/deliverables/sat_model_RMSD_and_plot.py +# for Model, use as basis "build_layer_map.py" +#mit_prex_or_die "mpirun python $MEDEAF_DIR/build_layer_map.py -i $NETCDF_DIR -o $OUTPUTDIR -d $MIT_RUNDATE -m $MASKFILE -p $PLOTLISTFILE " + +