Sat validation for mer

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',

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)
+

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()
+

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 "
+
+