Changeset 343 in lmdz_wrf

Timestamp:

Mar 9, 2015, 5:32:39 PM (10 years ago)

Author:

lfita

Message:

Adding 'nan' checking on single station observed values

File:

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

trunk/tools/validation_sim.py

@@ -11 +11 @@
 from optparse import OptionParser
 from netCDF4 import Dataset as NetCDFFile
+from scipy import stats as sts
+import numpy.ma as ma
 
 main = 'validarion_sim.py'
@@ -294 +296 @@
 
     return valpos
+
+def datetimeStr_datetime(StringDT):
+    """ Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
+    >>> datetimeStr_datetime('1976-02-17_00:00:00')
+    1976-02-17 00:00:00
+    """
+    import datetime as dt
+
+    fname = 'datetimeStr_datetime'
+
+    dateD = np.zeros((3), dtype=int)
+    timeT = np.zeros((3), dtype=int)
+
+    dateD[0] = int(StringDT[0:4])
+    dateD[1] = int(StringDT[5:7])
+    dateD[2] = int(StringDT[8:10])
+
+    trefT = StringDT.find(':')
+    if not trefT == -1:
+#        print '  ' + fname + ': refdate with time!'
+        timeT[0] = int(StringDT[11:13])
+        timeT[1] = int(StringDT[14:16])
+        timeT[2] = int(StringDT[17:19])
+
+    if int(dateD[0]) == 0:
+        print warnmsg
+        print '    ' + fname + ': 0 reference year!! changing to 1'
+        dateD[0] = 1 
+ 
+    newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])
+
+    return newdatetime
+
+def datetimeStr_conversion(StringDT,typeSi,typeSo):
+    """ Function to transform a string date to an another date object
+    StringDT= string with the date and time
+    typeSi= type of datetime string input
+    typeSo= type of datetime string output
+      [typeSi/o]
+        'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
+        'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
+        'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
+        'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
+        'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
+        'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
+        'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H], 
+          [H], ':', [M], [M], ':', [S], [S]
+    >>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
+    [1976    2   17    8   32    5]
+    >>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
+    1980/03/05 18-00-00
+    """
+    import datetime as dt
+
+    fname = 'datetimeStr_conversion'
+
+    if StringDT[0:1] == 'h':
+        print fname + '_____________________________________________________________'
+        print datetimeStr_conversion.__doc__
+        quit()
+
+    if typeSi == 'cfTime':
+        timeval = np.float(StringDT.split(',')[0])
+        tunits = StringDT.split(',')[1].split(' ')[0]
+        Srefdate = StringDT.split(',')[1].split(' ')[2]
+
+# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
+##
+        yrref=Srefdate[0:4]
+        monref=Srefdate[5:7]
+        dayref=Srefdate[8:10]
+
+        trefT = Srefdate.find(':')
+        if not trefT == -1:
+#            print '  ' + fname + ': refdate with time!'
+            horref=Srefdate[11:13]
+            minref=Srefdate[14:16]
+            secref=Srefdate[17:19]
+            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
+              '_' + horref + ':' + minref + ':' + secref)
+        else:
+            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
+              + '_00:00:00')
+
+        if tunits == 'weeks':
+            newdate = refdate + dt.timedelta(weeks=float(timeval))
+        elif tunits == 'days':
+            newdate = refdate + dt.timedelta(days=float(timeval))
+        elif tunits == 'hours':
+            newdate = refdate + dt.timedelta(hours=float(timeval))
+        elif tunits == 'minutes':
+            newdate = refdate + dt.timedelta(minutes=float(timeval))
+        elif tunits == 'seconds':
+            newdate = refdate + dt.timedelta(seconds=float(timeval))
+        elif tunits == 'milliseconds':
+            newdate = refdate + dt.timedelta(milliseconds=float(timeval))
+        else:
+              print errormsg
+              print '    timeref_datetime: time units "' + tunits + '" not ready!!!!'
+              quit(-1)
+
+        yr = newdate.year
+        mo = newdate.month
+        da = newdate.day
+        ho = newdate.hour
+        mi = newdate.minute
+        se = newdate.second
+    elif typeSi == 'matYmdHMS':
+        yr = StringDT[0]
+        mo = StringDT[1]
+        da = StringDT[2]
+        ho = StringDT[3]
+        mi = StringDT[4]
+        se = StringDT[5]
+    elif typeSi == 'YmdHMS':
+        yr = int(StringDT[0:4])
+        mo = int(StringDT[4:6])
+        da = int(StringDT[6:8])
+        ho = int(StringDT[8:10])
+        mi = int(StringDT[10:12])
+        se = int(StringDT[12:14])
+    elif typeSi == 'Y-m-d_H:M:S':
+        dateDT = StringDT.split('_')
+        dateD = dateDT[0].split('-')
+        timeT = dateDT[1].split(':')
+        yr = int(dateD[0])
+        mo = int(dateD[1])
+        da = int(dateD[2])
+        ho = int(timeT[0])
+        mi = int(timeT[1])
+        se = int(timeT[2])
+    elif typeSi == 'Y-m-d H:M:S':
+        dateDT = StringDT.split(' ')
+        dateD = dateDT[0].split('-')
+        timeT = dateDT[1].split(':')
+        yr = int(dateD[0])
+        mo = int(dateD[1])
+        da = int(dateD[2])
+        ho = int(timeT[0])
+        mi = int(timeT[1])
+        se = int(timeT[2])
+    elif typeSi == 'Y/m/d H-M-S':
+        dateDT = StringDT.split(' ')
+        dateD = dateDT[0].split('/')
+        timeT = dateDT[1].split('-')
+        yr = int(dateD[0])
+        mo = int(dateD[1])
+        da = int(dateD[2])
+        ho = int(timeT[0])
+        mi = int(timeT[1])
+        se = int(timeT[2])
+    elif typeSi == 'WRFdatetime':
+        yr = int(StringDT[0])*1000 + int(StringDT[1])*100 + int(StringDT[2])*10 +    \
+          int(StringDT[3])
+        mo = int(StringDT[5])*10 + int(StringDT[6])
+        da = int(StringDT[8])*10 + int(StringDT[9])
+        ho = int(StringDT[11])*10 + int(StringDT[12])
+        mi = int(StringDT[14])*10 + int(StringDT[15])
+        se = int(StringDT[17])*10 + int(StringDT[18])
+    else:
+        print errormsg
+        print '  ' + fname + ': type of String input date "' + typeSi +              \
+          '" not ready !!!!'
+        quit(-1)
+
+    if typeSo == 'matYmdHMS':
+        dateYmdHMS = np.zeros((6), dtype=int)
+        dateYmdHMS[0] =  yr
+        dateYmdHMS[1] =  mo
+        dateYmdHMS[2] =  da
+        dateYmdHMS[3] =  ho
+        dateYmdHMS[4] =  mi
+        dateYmdHMS[5] =  se
+    elif typeSo == 'YmdHMS':
+        dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) +        \
+          str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
+    elif typeSo == 'Y-m-d_H:M:S':
+        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
+          str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
+          str(se).zfill(2)
+    elif typeSo == 'Y-m-d H:M:S':
+        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
+          str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
+          str(se).zfill(2)
+    elif typeSo == 'Y/m/d H-M-S':
+        dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' +               \
+          str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
+          str(se).zfill(2) 
+    elif typeSo == 'WRFdatetime':
+        dateYmdHMS = []
+        yM = yr/1000
+        yC = (yr-yM*1000)/100
+        yD = (yr-yM*1000-yC*100)/10
+        yU = yr-yM*1000-yC*100-yD*10
+
+        mD = mo/10
+        mU = mo-mD*10
+        
+        dD = da/10
+        dU = da-dD*10
+
+        hD = ho/10
+        hU = ho-hD*10
+
+        miD = mi/10
+        miU = mi-miD*10
+
+        sD = se/10
+        sU = se-sD*10
+
+        dateYmdHMS.append(str(yM))
+        dateYmdHMS.append(str(yC))
+        dateYmdHMS.append(str(yD))
+        dateYmdHMS.append(str(yU))
+        dateYmdHMS.append('-')
+        dateYmdHMS.append(str(mD))
+        dateYmdHMS.append(str(mU))
+        dateYmdHMS.append('-')
+        dateYmdHMS.append(str(dD))
+        dateYmdHMS.append(str(dU))
+        dateYmdHMS.append('_')
+        dateYmdHMS.append(str(hD))
+        dateYmdHMS.append(str(hU))
+        dateYmdHMS.append(':')
+        dateYmdHMS.append(str(miD))
+        dateYmdHMS.append(str(miU))
+        dateYmdHMS.append(':')
+        dateYmdHMS.append(str(sD))
+        dateYmdHMS.append(str(sU))
+    else:
+        print errormsg
+        print '  ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
+        quit(-1)
+
+    return dateYmdHMS
 
 def coincident_CFtimes(tvalB, tunitA, tunitB):
@@ -576 +813 @@
         'WRFp': pressure from WRF variables
         'WRFrh': relative humidty fom WRF variables
+        'WRFT': CF-time from WRF variables
         'WRFt': temperature from WRF variables
         'WRFtd': dew-point temperature from WRF variables
@@ -645 +883 @@
                 shape = ncobj.variables['P'].shape
 
+            elif varn == 'WRFT':
+# To compute CF-times from WRF kind
+#
+                import datetime as dt
+
+                times = ncobj.variables['Times']
+                dimt = times.shape[0]
+                varNOcheckv = np.zeros((dimt), dtype=np.float64)
+                self.unitsval = 'seconds since 1949-12-01 00:00:00'
+                refdate = datetimeStr_datetime('1949-12-01_00:00:00')
+
+                dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
+                shape = tuple([dimt])
+
+                for it in range(dimt):
+                    datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime',    \
+                      'YmdHMS')
+                    dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
+                    difft = dateval - refdate
+                    varNOcheckv[it] = difft.days*3600.*24. + difft.seconds +        \
+                          np.float(int(difft.microseconds/10.e6))
+
             elif varn == 'WRFt':
 #        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
@@ -740 +1000 @@
   '[constant] to variables values; divc,[constant]: divide by [constant] to ' +      \
   'variables values'
-varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFt', 'WRFtd', 'WRFws',        \
+varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFT', 'WRFt', 'WRFtd', 'WRFws',\
   'WRFwd', 'WRFz']
 varNOcheckinf = "'WRFdens': air density from WRF variables; 'WRFght': geopotential"+ \
   " height from WRF variables; 'WRFp': pressure from WRF variables; 'WRFrh': " +     \
-  "relative humidty fom WRF variables; 'WRFt': temperature from WRF variables; " +   \
-  "'WRFtd': dew-point temperature from WRF variables; 'WRFws': wind speed from " +   \
-  "WRF variables; 'WRFwd': wind speed direction from WRF variables; 'WRFz': " +     \
-  "height from WRF variables"
+  "relative humidty fom WRF variables; 'WRFT': CF-time from WRF variables'WRFt'; " + \
+  " temperature from WRF variables; 'WRFtd': dew-point temperature from WRF " +      \
+  "variables; 'WRFws': wind speed from WRF variables; 'WRFwd': wind speed " +        \
+  "direction from WRF variables; 'WRFz': height from WRF variables"
 
 dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
@@ -925 +1185 @@
     if not osim.dimensions.has_key(dims[dn][0]):
         print errormsg
-        print '  ' + main + ": simulation file does not have dimension '" +          \
-          dims[dn][0] + "' !!"
+        print '  ' + main + ": simulation file '" + opts.fsim + "' does not have " + \
+          "dimension '" + dims[dn][0] + "' !!"
+        print '    it has: ',osim.dimensions
         quit(-1)
+
     if not osim.variables.has_key(vardims[dn][0]) and not                            \
       searchInlist(varNOcheck,vardims[dn][0]):
         print errormsg
-        print '  ' + main + ": simulation file does not have varibale dimension '" + \
-          vardims[dn][0] + "' !!"
+        print '  ' + main + ": simulation file '" + opts.fsim + "' does not have " + \
+          "varibale dimension '" + vardims[dn][0] + "' !!"
+        print '    it has variables:',osim.variables
         quit(-1)
     if searchInlist(varNOcheck,vardims[dn][0]):
@@ -940 +1203 @@
 
 # General characteristics
-dimtobs = len(valdimobs['T'])
-dimtsim = len(valdimsim['T'])
+dimtobs = valdimobs['T'].shape[0]
+dimtsim = valdimsim['T'].shape[0]
 
 print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
@@ -1026 +1289 @@
 tobj = oobs.variables[vardims['T'][1]]
 obstunits = tobj.getncattr('units')
-tobj = osim.variables[vardims['T'][0]]
-simtunits = tobj.getncattr('units')
-
-simobstimes = coincident_CFtimes(valdimsim['T'], obstunits, simtunits)
+if vardims['T'][0] == 'WRFT':
+    tsim = valdimsim['T'][:]
+    simtunits = 'seconds since 1949-12-01 00:00:00'
+else:
+    tsim = osim.variables[vardims['T'][0]]
+    simtunits = tobj.getncattr('units')
+
+simobstimes = coincident_CFtimes(valdimsim['T'][:], obstunits, simtunits)
 
 # Concident times
@@ -1157 +1424 @@
     elif obskind == 'single-station':
         for it in range(Ncoindt):
-            slicev = dims['X'][0]+':'+str(stationpos[1]) + '|' +                      \
-              dims['Y'][0]+':'+str(stationpos[0]) + '|' +                             \
-              dims['T'][0]+':'+str(coindtvalues[it][0])
+            ito = int(coindtvalues[it,1])
+            slicev = dims['X'][0] + ':' + str(stationpos[1]) + '|' +                 \
+              dims['Y'][0] + ':' + str(stationpos[0]) + '|' +                        \
+              dims['T'][0] + ':' + str(int(coindtvalues[it][0]))
             slicevar, dimslice = slice_variable(ovsim, slicev)
-            slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
-              str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
-              str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '|' +      \
-              dims['T'][0] + ':' + str(coindtvalues[it,0])
+            if ovobs[int(ito)] != ovobs[int(ito)]: 
+                simobsvalues.append([ slicevar, fillValueF])
+            else:
+                simobsvalues.append([ slicevar, ovobs[int(ito)]])
+            slicev = dims['X'][0] + ':' + str(stationpos[1]-Ngrid) + '@' +           \
+              str(stationpos[1]+Ngrid+1) + '|' + dims['Y'][0] + ':' +                \
+              str(stationpos[0]-Ngrid) + '@' + str(stationpos[0]+Ngrid+1) + '|' +    \
+              dims['T'][0] + ':' + str(int(coindtvalues[it,0]))
             slicevar, dimslice = slice_variable(ovsim, slicev)
             simobsSvalues[it,:,:] = slicevar
+
     elif obskind == 'trajectory':
         if dims.has_key('Z'):
@@ -1250 +1523 @@
             quit(-1)
 
+# statisics sim
+    simstats = np.zeros((5), dtype=np.float)
+    simstats[0] = np.min(arrayvals[:,0])
+    simstats[1] = np.max(arrayvals[:,0])
+    simstats[2] = np.mean(arrayvals[:,0])
+    simstats[3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
+    simstats[4] = np.sqrt(simstats[3] - simstats[2]*simstats[2])
+
+# statisics obs
+    obsmask = ma.masked_equal(arrayvals, fillValueF)
+    obsmask2 = obsmask*obsmask
+
+    print 'Lluis',obsmask
+
+    obsstats = np.zeros((5), dtype=np.float)
+    obsstats[0] = obsmask.min()
+    obsstats[1] = obsmask.max()
+    obsstats[2] = obsmask.mean()
+    obsstats[3] = obsmask2.mean()
+    obsstats[4] = np.sqrt(obsstats[3] - obsstats[2]*obsstats[2])
+
+# Statistics sim-obs
+    simobsstats = np.zeros((5), dtype=np.float)
+    diffvals = np.zeros((Ncoindt), dtype=np.float)
+
+    diffvals = arrayvals[:,0] - arrayvals[:,1]
+    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
     aroundstats = np.zeros((5,Ncoindt), dtype=np.float)