mvr on o2o n3n

src/bitsea/instruments/check.py

@@ -20,8 +20,11 @@ def __init__(self,OUTDIR:Path, verboselevel=1, threshold_nitrate=1):
         if verboselevel = 1, a NetCDF will be printed out
         else, no files will be dumped
         '''
-        self.outdir = OUTDIR
-        OUTDIR.mkdir(parents=True, exist_ok=True)
+        if OUTDIR=="":
+            self.outdir=""
+        else:    
+            self.outdir = OUTDIR
+            OUTDIR.mkdir(parents=True, exist_ok=True)
         self.verboselevel=verboselevel
         self.threshold_nitrate = threshold_nitrate
 

src/bitsea/postproc/aveScan.py

@@ -415,7 +415,7 @@ def dumpPointProfiles(var):
     ncOUT_Pprofiles.createDimension("z"         ,jpk)            
     global VAR
     #VAR = VAR.copy()
-    VAR[~tmask]=np.NaN    
+    VAR[~tmask]=np.nan    
     P=PointProfiles(var)
     ncvar    = ncOUT_Pprofiles.createVariable(var  ,'f',('Ncruise','z'))    
     ncvar[:] = P

src/bitsea/validation/online/MVR_gen.py

@@ -50,6 +50,7 @@ def argument():
 from bitsea.commons import genUserDateList as DL
 from bitsea.commons import timerequestors as requestors
 import calendar
+import sys
 
 INDIRSAT = addsep(args.inputsat)
 INDIRFLOAT = addsep(args.inputfloat)
@@ -73,6 +74,8 @@ def argument():
 TLfloat = TimeList.fromfilenames(TI,INDIRFLOAT,"BioFloat_Weekly*nc", \
                     prefix="BioFloat_Weekly_validation_", \
                     dateformat='%Y%m%d')
+
+#### Sat
 TLsat = {}
 
 LISTforecast = ['f%d' %ii for ii in range(1,4)]
@@ -93,15 +96,17 @@ def argument():
                         prefix="Validation_" + ff + "_YYYYMMDD_on_daily_Sat.", \
                         dateformat='%Y%m%d')
 
+#### end Sat
+
 METRICS_NAMES = [
     'number of data values',
     'mean of product',
     'mean of reference',
     'mean squared error',
     'variance of product',
     'variance of reference',
-    'correlation',
     'anomaly correlation',
+    'correlation',
     ]
 
 
@@ -145,8 +150,8 @@ def argument():
     DICTmetricnames[tt]['mean squared error'] = 'BGC_CLASS4_CHL_RMS_SURF_BASIN' + stringlog
     DICTmetricnames[tt]['variance of product']   = 'MODEL_VARIANCE' + stringlog
     DICTmetricnames[tt]['variance of reference'] = 'SAT___VARIANCE' + stringlog
-    DICTmetricnames[tt]['correlation'] = 'BGC_CLASS4_CHL_CORR_SURF_BASIN' + stringlog
     DICTmetricnames[tt]['anomaly correlation'] = 'ANOMALY_CORRELATION' + stringlog
+    DICTmetricnames[tt]['correlation'] = 'BGC_CLASS4_CHL_CORR_SURF_BASIN' + stringlog
 
 
 file1 = TLsat['f1'].filelist[0]
@@ -176,7 +181,7 @@ def argument():
 
     for idate,datef in enumerate(AllDates):
         req = requestors.Daily_req(datef.year,datef.month,datef.day)
-        TIMELIST[idate] = datef.toordinal() - datetime.datetime(1970,1,1).toordinal()
+        TIMELIST[idate] = datef.toordinal() - datetime.datetime(1950,1,1).toordinal()
         for iforecast,ff in enumerate(LISTforecast):
             ii = TLsat[ff].select_one(req)
             if ii==None: continue
@@ -200,11 +205,11 @@ def argument():
 
 LISTmetrics_FLOAT = [
     'nobs',
-    'refmean',
     'modelmean',
-    'refvar',
+    'refmean',
+    'mse',
     'modelvar',
-    'rmse',
+    'refvar',
     'bias',
     'BGC_CLASS4_CHL_CORR_SURF_BASIN_LOG',
     'anomaly_corr',
@@ -219,11 +224,11 @@ def argument():
 DICTmetricnames['number of data values'] = 'nobs'
 DICTmetricnames['mean of product']   = 'modelmean'
 DICTmetricnames['mean of reference'] = 'refmean'
-DICTmetricnames['mean squared error'] = 'rmse'
+DICTmetricnames['mean squared error'] = 'mse'
 DICTmetricnames['variance of product']   = 'modelvar'
 DICTmetricnames['variance of reference'] = 'refvar'
-DICTmetricnames['correlation']         = 'corr'
 DICTmetricnames['anomaly correlation'] = 'anomaly_corr'
+DICTmetricnames['correlation']         = 'corr'
 
 
 DICTsubbasins = {}
@@ -253,108 +258,197 @@ def argument():
 DEPTHSlist = []
 for layername in M_depths:
     top,bottom_m = layername.split('-')
-    if float(top)>150 : continue
+    #if float(top)>150 : continue
     bottom = float(bottom_m[:-1])
     DEPTHSlist.append(bottom)
 Ndepths = len(DEPTHSlist)
 
 Nmetrics = len(METRICS_NAMES)
 Nsub = len(AREA_NAMES)
 
-MetricsFLOAT = np.zeros((Ndates,Ndepths,Nmetrics,Nsub))
-MetricsFLOAT[:] = np.nan
-TIMELIST = np.zeros((Ndates,),np.float32)
-
+varlist=['P_l','N3n','O2o']
+Nvar=len(varlist)
 
+MetricsFLOAT = np.zeros((Ndates,Ndepths,Nmetrics,Nsub,Nvar))
+MetricsFLOAT[:] = 1.e+20
+TIMELIST = np.zeros((Ndates,),np.float32)
 
 for idate,datef in enumerate(AllDates):
-    TIMELIST[idate] = datef.toordinal() - datetime.datetime(1970,1,1).toordinal()
+    TIMELIST[idate] = datef.toordinal() - datetime.datetime(1950,1,1).toordinal()
     ii,diffseconds = TLfloat.find(datef,returndiff=True)
     if diffseconds/24/3600>3.5: continue
     filef = TLfloat.filelist[ii]
-
     M = NC.Dataset(filef,"r")
     for mm,metric in enumerate(METRICS_NAMES):
         for isub,subname in enumerate(AREA_NAMES):
             if subname=='Aegean Sea': continue
             varname = DICTmetricnames[0][metric].decode()
-            MetricsFLOAT[idate,:,mm,isub] = M.variables[varname][0,indexAREAS[subname],:Ndepths].data.copy()
+            for ivar, VAR in enumerate (varlist):
+                MetricsFLOAT[idate,:,mm,isub,ivar] = M.variables[varname][ivar,indexAREAS[subname],:Ndepths].data.copy()
 
 masknan = np.isnan(MetricsFLOAT)
 MetricsFLOAT[masknan] = 1.e+20
 
+
 outfile = 'product_quality_stats_' + productname + '_' + STARTTIME + '_' + END__TIME + '.nc'
 outfilepath = OUTDIR + outfile
+
+print(outfilepath)
+
+#import sys
+#sys.exit()
+
 S=NC.Dataset(outfilepath,"w")
 
 S.createDimension("time"         ,None)
-S.createDimension("string_length",max_strlenght)
-S.createDimension("areas"        ,len(AREA_NAMES))
-S.createDimension("metrics"      ,len(METRICS_NAMES))
-S.createDimension("depths"       ,Ndepths)
-S.createDimension("surface"      ,1)
-S.createDimension("forecasts"    ,len(leadtimes))
-
+S.createDimension("area"        ,len(AREA_NAMES))
+S.createDimension("metric"      ,len(METRICS_NAMES))
+S.createDimension("depth"       ,1)
+S.createDimension("layer"       ,Ndepths)
+S.createDimension("forecast"    ,len(leadtimes))
+
+subbasins_renamed = [
+    "Alboran Sea",
+    "South West Med western part",
+    "South West Med eastern part",
+    "North West Med",
+    "Tyrrhenian Sea northern part",
+    "Tyrrhenian Sea southern part",
+    "Adriatic Sea northern part",
+    "Adriatic Sea southern part",
+    "Aegean Sea",
+    "Ionian Sea western part",
+    "Ionian Sea south-eastern part",
+    "Ionian Sea north-eastern part",
+    "Levantine Sea western part",
+    "Levantine Sea central-northern part",
+    "Levantine Sea central-southern part",
+    "Levantine Sea eastern part",
+    "Full domain"
+]
 
-ncvar = S.createVariable("area_names",'c',('areas','string_length'))
 x=np.array(ALLstring,dtype=str)
-ncvar[:] = x.view('U1').reshape(x.size,-1)[:Nsub,:]
-setattr(S.variables['area_names'],"long_name"  ,"area names")
-setattr(S.variables['area_names'],"description","region over which statistics are aggregated")
-
-ncvar = S.createVariable("metric_names",'c',('metrics','string_length'))
-ncvar[:] = x.view('U1').reshape(x.size,-1)[Nsub:,:]
-setattr(S.variables['metric_names'],"long_name","metric names")
 
-ncvar = S.createVariable("forecasts",'f',('forecasts',))
+from netCDF4 import stringtochar
+# writing AREA as Variable
+vlen_str = S.createVLType(str, "vlen_str")
+ncvar = S.createVariable("area", vlen_str, ("area",))
+for i, name in enumerate(subbasins_renamed):
+    ncvar[i] = name
+
+#ncvar = S.createVariable("area", str, ("area",)) # 
+#x=np.array(subbasins_renamed ,dtype=str)
+#ncvar[:] = subbasins_renamed  #x[:Nsub]
+setattr(S.variables['area'],"long_name"  ,"area")
+setattr(S.variables['area'], "description", "Geographical areas that are included in the Region of reference")
+
+# writing metric as Variable
+ncvar = S.createVariable("metric", str, ('metric'))
+ncvar[:] = x[Nsub:]
+setattr(S.variables['metric'], "long_name", "List of CLASS4 metrics")
+
+# writing forecast as Variable
+ncvar = S.createVariable("forecast",'f4',('forecast',) , fill_value=-999.0)
 ncvar[:] = leadtimes
-setattr(S.variables['forecasts'],"long_name","forecast lead time")
-setattr(S.variables['forecasts'],"units"    ,"hours")
+setattr(S.variables['forecast'],"description","forecast lead time")
+setattr(S.variables['forecast'],"units"    ,"hours")
+
+# writing depth as Variable
+ncvar = S.createVariable("depth",'f',('depth',))
+ncvar[:]=np.array(0)
+setattr(S.variables['depth'],"long_name"  ,"depth")
+setattr(S.variables['depth'],"units"      ,"m")
+setattr(S.variables['depth'], "description", "Surface level (0 m) over which statistics are calculated")
+
 
-ncvar = S.createVariable("depths",'f',('depths',))
+# writing layer as Variable 
+ncvar = S.createVariable("layer",'f',('layer',))
 ncvar[:]=np.array(DEPTHSlist)
-setattr(S.variables['depths'],"long_name"  ,"depths")
-setattr(S.variables['depths'],"positive"   ,"down")
-setattr(S.variables['depths'],"units"      ,"m")
-setattr(S.variables['depths'],"description","depth of the base of the vertical layer over which statistics are aggregated")
+setattr(S.variables['layer'],"long_name"  ,"layer")
+setattr(S.variables['layer'],"positive"   ,"down")
+setattr(S.variables['layer'],"units"      ,"m")
+setattr(S.variables['layer'],"description","depth of the base of the vertical layer over which statistics are aggregated")
 
-ncvar = S.createVariable("time",'f',('time',))
+# writing time as Variable
+ncvar = S.createVariable('time', 'f4', ('time',), fill_value=-999.0)
 ncvar[:]=TIMELIST
 setattr(S.variables['time'],"long_name"    ,"validity time")
 setattr(S.variables['time'],"standard_name","time")
-setattr(S.variables['time'],"units"        ,"days since 1970-01-01")
+setattr(S.variables['time'],"units"        ,"days since 1950-01-01 12:00:00")
 setattr(S.variables['time'],"axis"         ,"T")
 
-for tt in ['LOG','NOLOG']:
+# writing stats_chlorophyll-a_sat as Variable ONLY NOLOG
+for tt in ['NOLOG']:  #for tt in ['LOG','NOLOG']:
     if tt=='LOG':
         stringlog = '_' + tt
     else:
         stringlog = ''
-    statsname = 'stats_surface_chlorophyll' + stringlog.lower()
+    statsname = 'stats_chlorophyll-a_sat-l3' + stringlog.lower()
     if tt=='LOG':
-        parametername = 'log of Surface Chlorophyll'
+        parametername = 'log of Surface chlorophyll'
     else:
         parametername = "Surface Chlorophyll"
-    ncvar=S.createVariable(statsname,'f',('time','forecasts','surface','metrics','areas'),fill_value=1.e+20)
+    ncvar = S.createVariable(statsname, 'f4', ('time', 'forecast', 'depth', 'metric', 'area'), fill_value=-999.0)
     ncvar[:,:,0,:,:] = MetricsSAT[tt]
     setattr(S.variables[statsname], "parameter",parametername)
-    setattr(S.variables[statsname], "reference","Satellite observations from OC-TAC")
-    setattr(S.variables[statsname], "units"    ,"log(mg Chl*m^-3)")
+    setattr(S.variables[statsname], "reference","OCEANCOLOUR_MED_BGC_L3_NRT_009_141")
+    setattr(S.variables[statsname], "metric_name","CHL-SURF-D-CLASS4-SAT-PQD_EAN-XYZ-MED")
+    setattr(S.variables[statsname], "units"    ,"mg Chl*m^-3") #log removed
+
 
-statsname = 'stats_profile_chlorophyll'
+###### saving FLOAT
+# writing Chlorophyll as BGC-Argo  Variable
+statsname = 'stats_chlorophyll-a_ins-pf'
 parametername = "Chlorophyll"
-ncvar=S.createVariable(statsname,'f',('time','forecasts','depths','metrics','areas'),fill_value=1.e+20)
-ncvar[:,:,:,:,:] = np.nan
-ncvar[:,0,:,:,:] = MetricsFLOAT
+ncvar=S.createVariable(statsname,'f4',('time', 'forecast', 'layer', 'metric', 'area'), fill_value=-999.0)
+
+ncvar[:,:,:,:,:]   = np.nan
+ncvar[:,0,:,:,:]   = MetricsFLOAT[:,:,:,:,0]
+ncvar[:,:,-3:,:,:] = np.nan
 
 setattr(S.variables[statsname], "parameter",parametername)
-setattr(S.variables[statsname], "reference","BGC-Argo floats")
+setattr(S.variables[statsname], "reference","BGC-Argo")
+setattr(S.variables[statsname], "metric_name","CHL-X_Y-D-CLASS4-PROF-PQD_EAN-XY-MED")
 setattr(S.variables[statsname], "units"    ,"mg Chl*m^-3")
 
+
+# writing nitrate as BGC-Argo  Variable
+statsname = 'stats_nitrate_ins-pf'
+parametername = "Nitrate"
+ncvar = S.createVariable(statsname, 'f4', ('time', 'forecast', 'layer', 'metric', 'area'), fill_value=-999.0)
+
+ncvar[:,:,:,:,:] = np.nan
+ncvar[:,0,:,:,:] = MetricsFLOAT[:,:,:,:,1]
+
+setattr(S.variables[statsname], "parameter",parametername)
+setattr(S.variables[statsname], "reference","BGC-Argo")
+setattr(S.variables[statsname], "metric_name","NO3-X_Y-D-CLASS4-PROF-PQD_EAN-XY-MED")
+setattr(S.variables[statsname], "units"    ,"mmol NO3*m^-3")
+
+# writing Oxygen as BGC-Argo  Variable
+statsname = 'stats_oxygen_ins-pf'
+parametername = "Oxygen"
+ncvar = S.createVariable(statsname, 'f4', ('time', 'forecast', 'layer', 'metric', 'area'), fill_value=-999.0)
+ncvar[:,:,:,:,:] = np.nan
+ncvar[:,0,:,:,:] = MetricsFLOAT[:,:,:,:,2]
+
+
+setattr(S.variables[statsname], "parameter",parametername)
+setattr(S.variables[statsname], "reference","BGC-Argo")
+setattr(S.variables[statsname], "metric_name","O2-X_Y-D-CLASS4-PROF-PQD_EAN-XY-MED")
+setattr(S.variables[statsname], "units"    ,"mmol O2*m^-3")
+
+
+# writing metadata
 setattr(S,"contact","[email protected]")
 setattr(S,"product",productname)
 setattr(S,"start_date",STARTTIME)
 setattr(S,"end_date"  ,END__TIME)
-setattr(S,"filename"  ,outfile)
+setattr(S, "institution"  , "OGS")  
+from datetime import date
+today = date.today().isoformat()
+setattr(S, "bulletin_date", today)
+setattr(S, "frequency"    , "monthly")
+#setattr(S,"filename"  ,outfile)
 S.close()
 

src/bitsea/validation/online/RunSatValidation.py

@@ -95,7 +95,7 @@ def argument():
 
 def climfilename(CLIM_DIR, month,var):
     if var=='P_l':
-        return CLIM_DIR + '/ave.yyyy' + month + '15-00:00:00.nc'
+        return CLIM_DIR + f"/ave.yyyy{month}15-00:00:00.{var}.nc"
     else:
         return None