Changeset 344 in lmdz_wrf

Timestamp:

Mar 9, 2015, 7:02:22 PM (10 years ago)

Author:

lfita

Message:

Adding statistics for sim, obs, simobs and around temporally from observations

File:

• trunk/tools/validation_sim.py (modified) (10 diffs)

trunk/tools/validation_sim.py

@@ -1018 +1018 @@
   "(WRFdiagnosted variables also available: " + varNOcheckinf + ")"
 statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
+gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae',  \
+  'rmse', 'r_pearsoncorr', 'p_pearsoncorr']
+ostatsn = ['number of points', 'start', 'end', 'minimum', 'maximum', 'mean',         \
+  'mean2', 'standard deviation']
 
 parser = OptionParser()
@@ -1334 +1338 @@
 newdim = onewnc.createDimension('xaround',Ngrid*2+1)
 newdim = onewnc.createDimension('yaround',Ngrid*2+1)
+newdim = onewnc.createDimension('gstats',9)
 newdim = onewnc.createDimension('stats',5)
+newdim = onewnc.createDimension('tstats',8)
 
 # Variable dimensions
@@ -1350 +1356 @@
 basicvardef(newvar, 'statistics', 'statitics from values', '-')
 writing_str_nc(newvar, statsn, StringLength)
+
+newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
+basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
+writing_str_nc(newvar, gstatsn, StringLength)
+
+newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
+basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
+writing_str_nc(newvar, ostatsn, StringLength)
 
 if obskind == 'trajectory':
@@ -1402 +1416 @@
 
     ovobs = oobs.variables[valvars[ivar][1]]
+
+# Simulated values spatially around coincident times
     if dims.has_key('Z'):
         simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1, Ngrid*2+1),         \
@@ -1407 +1423 @@
     else:
         simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1), dtype = np.float)
+
+# Observed values temporally around coincident times
+    simobsTvalues = {}
+    simobsTtvalues = np.zeros((Ncoindt,2), dtype=np.float)
 
     if obskind == 'multi-points':
@@ -1439 +1459 @@
             slicevar, dimslice = slice_variable(ovsim, slicev)
             simobsSvalues[it,:,:] = slicevar
+
+            if it == 0:
+                itoi = 0
+                itof = int(coindtvalues[it,1]) / 2
+            elif it == Ncoindt-1:
+                itoi = int( (ito + int(coindtvalues[it-1,1])) / 2)
+                itof = int(coindtvalues[it,1])
+            else:
+                itod = int( (ito - int(coindtvalues[it-1,1])) / 2 ) 
+                itoi = ito - itod
+                itod = int( (int(coindtvalues[it+1,1]) - ito) / 2 )
+                itof = ito + itod
+
+            slicev = dims['T'][1] + ':' + str(itoi) + '@' + str(itof + 1)
+
+            slicevar, dimslice = slice_variable(ovobs, slicev)
+            simobsTvalues[str(it)] = slicevar
+
+            simobsTtvalues[it,0] = valdimobs['T'][itoi]
+            simobsTtvalues[it,1] = valdimobs['T'][itof]
 
     elif obskind == 'trajectory':
@@ -1532 +1572 @@
 
 # statisics obs
-    obsmask = ma.masked_equal(arrayvals, fillValueF)
+    obsmask = ma.masked_equal(arrayvals[:,1], fillValueF)
     obsmask2 = obsmask*obsmask
-
-    print 'Lluis',obsmask
 
     obsstats = np.zeros((5), dtype=np.float)
@@ -1545 +1583 @@
 
 # Statistics sim-obs
-    simobsstats = np.zeros((5), dtype=np.float)
+    simobsstats = np.zeros((9), dtype=np.float)
     diffvals = np.zeros((Ncoindt), dtype=np.float)
 
-    diffvals = arrayvals[:,0] - arrayvals[:,1]
+    diffvals = arrayvals[:,0] - obsmask
+
     simobsstats[0] = simstats[0] - obsstats[0]
-    simobsstats[1] = np.mean(np.abs(diffvals))
-    simobsstats[2] = np.sqrt(simstats[3] + obsstats[3] - 2.*simstats[2]*obsstats[2])
-    simobsstats[3], simobsstats[4] = sts.pearsonr(arrayvals[:,0], arrayvals[:,1])
-
-    print 'Lluis means sim:',simstats[0],'obs:',obsstats[0]
-    print 'Lluis tests rmse:',simobsstats[2], np.sqrt(np.mean(diffvals*diffvals))
-
-# Statistics around values
+    simobsstats[1] = np.mean(arrayvals[:,0]*obsmask)
+    simobsstats[2] = np.min(diffvals)
+    simobsstats[3] = np.max(diffvals)
+    simobsstats[4] = np.mean(diffvals)
+    simobsstats[5] = np.mean(np.abs(diffvals))
+    simobsstats[6] = np.sqrt(np.mean(diffvals*diffvals))
+    simobsstats[7], simobsstats[8] = sts.pearsonr(arrayvals[:,0], arrayvals[:,1])
+
+# Statistics around sim values
     aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
     for it in range(Ncoindt):
@@ -1566 +1606 @@
         aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]*           \
           aroundstats[2,it])
+
+# Statistics around obs values
+    aroundostats = np.zeros((8,Ncoindt), dtype=np.float)
+
+    for it in range(Ncoindt):
+        obsmask = ma.masked_equal(simobsTvalues[str(it)], fillValueF)
+        obsmask2 = obsmask*obsmask
+
+        aroundostats[0,it] = len(obsmask.flatten())
+        aroundostats[1,it] = simobsTtvalues[it,0]
+        aroundostats[2,it] = simobsTtvalues[it,1]
+        aroundostats[3,it] = obsmask.min()
+        aroundostats[4,it] = obsmask.max()
+        aroundostats[5,it] = obsmask.mean()
+        aroundostats[6,it] = obsmask2.mean()
+        aroundostats[7,it] = np.sqrt(aroundostats[6,it] - aroundostats[5,it]*        \
+          aroundostats[5,it])
 
 # Values to netCDF
@@ -1590 +1647 @@
         newvar[:] = simobsSvalues
 
-# Statistics
-    if not searchInlist(onewnc.variables,valvars[ivar][0]):
+# sim Statistics
+    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
+        newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('stats'),   \
+          fill_value=fillValueF)
+        descvar = 'simulated statisitcs: ' + valvars[ivar][0]
+        basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, ovobs.getncattr('units'))
+        newvar[:] = simstats
+
+# obs Statistics
+    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
+        newvar = onewnc.createVariable(valvars[ivar][1] + 'stobs', 'f', ('stats'),   \
+          fill_value=fillValueF)
+        descvar = 'observed statisitcs: ' + valvars[ivar][1]
+        basicvardef(newvar, valvars[ivar][1] + 'stobs', descvar,                     \
+          ovobs.getncattr('units'))
+        newvar[:] = obsstats
+
+# sim-obs Statistics
+    if not searchInlist(onewnc.variables,varsimobs + 'st'):
+        newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('gstats'),            \
+          fill_value=fillValueF)
+        descvar = 'simulated-observed statisitcs: ' + varsimobs
+        basicvardef(newvar, varsimobs + 'st', descvar, ovobs.getncattr('units'))
+        newvar[:] = simobsstats
+
+# around sim Statistics
+    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'staround'):
         newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f',           \
           ('time','stats'), fill_value=fillValueF)
-        descvar = 'around simulated statisitcs: ' + valvars[ivar][0]
-        basicvardef(newvar, varsimobs + 'staround', descvar, ovobs.getncattr('units'))
+        descvar = 'around (' +  str(Ngrid) + ', ' + str(Ngrid) +                     \
+          ') simulated statisitcs: ' + valvars[ivar][0]
+        basicvardef(newvar, valvars[ivar][0] + 'staround', descvar,                  \
+          ovobs.getncattr('units'))
         newvar[:] = aroundstats.transpose()
+
+# around obs Statistics
+    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'staround'):
+        newvar = onewnc.createVariable(valvars[ivar][1] + 'staround', 'f',           \
+          ('time','tstats'), fill_value=fillValueF)
+        descvar = 'around temporal observed statisitcs: ' + valvars[ivar][1]
+        basicvardef(newvar, valvars[ivar][1] + 'staround', descvar,                  \
+          ovobs.getncattr('units'))
+        newvar[:] = aroundostats.transpose()
 
         onewnc.sync()