# L. Fita, LMD-Jussieu. February 2015 ## e.g. sfcEneAvigon # validation_sim.py -d X@west_east@None,Y@south_north@None,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,T@time@time -k single-station -l 4.878773,43.915876,12. -o /home/lluis/DATA/obs/HyMeX/IOP15/sfcEnergyBalance_Avignon/OBSnetcdf.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v HFX@H,LH@LE,GRDFLX@G ## e.g. AIREP # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@alti,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/AIREP/2012/10/AIREP_121018.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@t,WRFtd@td,WRFws@u,WRFwd@dd ## e.g. ATRCore # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/ATRCore/V3/ATR_1Hz-HYMEXBDD-SOP1-v3_20121018_as120051.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@air_temperature@subc@273.15,WRFp@air_pressure,WRFrh@relative_humidity,WRFrh@relative_humidity_Rosemount,WRFwd@wind_from_direction,WRFws@wind_speed ## e.g. BAMED # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/BAMED/BAMED_SOP1_B12_TOT5.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@tas_north,WRFp@pressure,WRFrh@hus,U@uas,V@vas import numpy as np import os import re 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' errormsg = 'ERROR -- errror -- ERROR -- error' warnmsg = 'WARNING -- warning -- WARNING -- warning' # version version=1.0 # Filling values for floats, integer and string fillValueF = 1.e20 fillValueI = -99999 fillValueS = '---' StringLength = 50 # Number of grid points to take as 'environment' around the observed point Ngrid = 1 def searchInlist(listname, nameFind): """ Function to search a value within a list listname = list nameFind = value to find >>> searInlist(['1', '2', '3', '5'], '5') True """ for x in listname: if x == nameFind: return True return False def set_attribute(ncvar, attrname, attrvalue): """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists ncvar = object netcdf variable attrname = name of the attribute attrvalue = value of the attribute """ import numpy as np from netCDF4 import Dataset as NetCDFFile attvar = ncvar.ncattrs() if searchInlist(attvar, attrname): attr = ncvar.delncattr(attrname) attr = ncvar.setncattr(attrname, attrvalue) return ncvar def basicvardef(varobj, vstname, vlname, vunits): """ Function to give the basic attributes to a variable varobj= netCDF variable object vstname= standard name of the variable vlname= long name of the variable vunits= units of the variable """ attr = varobj.setncattr('standard_name', vstname) attr = varobj.setncattr('long_name', vlname) attr = varobj.setncattr('units', vunits) return def writing_str_nc(varo, values, Lchar): """ Function to write string values in a netCDF variable as a chain of 1char values varo= netCDF variable object values = list of values to introduce Lchar = length of the string in the netCDF file """ Nvals = len(values) for iv in range(Nvals): stringv=values[iv] charvals = np.chararray(Lchar) Lstr = len(stringv) charvals[Lstr:Lchar] = '' for ich in range(Lstr): charvals[ich] = stringv[ich:ich+1] varo[iv,:] = charvals return def index_3mat(matA,matB,matC,val): """ Function to provide the coordinates of a given value inside three matrix simultaneously index_mat(matA,matB,matC,val) matA= matrix with one set of values matB= matrix with the other set of values matB= matrix with the third set of values val= triplet of values to search >>> index_mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),np.arange(200,227).reshape(3,3,3),[22,122,222]) [2 1 1] """ fname = 'index_3mat' matAshape = matA.shape matBshape = matB.shape matCshape = matC.shape for idv in range(len(matAshape)): if matAshape[idv] != matBshape[idv]: print errormsg print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \ 'and B:',matBshape[idv],'does not coincide!!' quit(-1) if matAshape[idv] != matCshape[idv]: print errormsg print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \ 'and C:',matCshape[idv],'does not coincide!!' quit(-1) minA = np.min(matA) maxA = np.max(matA) minB = np.min(matB) maxB = np.max(matB) minC = np.min(matC) maxC = np.max(matC) if val[0] < minA or val[0] > maxA: print warnmsg print ' ' + fname + ': first value:',val[0],'outside matA range',minA,',', \ maxA,'!!' if val[1] < minB or val[1] > maxB: print warnmsg print ' ' + fname + ': second value:',val[1],'outside matB range',minB,',', \ maxB,'!!' if val[2] < minC or val[2] > maxC: print warnmsg print ' ' + fname + ': second value:',val[2],'outside matC range',minC,',', \ maxC,'!!' dist = np.zeros(tuple(matAshape), dtype=np.float) dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2 + \ (matC - np.float(val[2]))**2) mindist = np.min(dist) matlist = list(dist.flatten()) ifound = matlist.index(mindist) Ndims = len(matAshape) valpos = np.zeros((Ndims), dtype=int) baseprevdims = np.zeros((Ndims), dtype=int) for dimid in range(Ndims): baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims]) if dimid == 0: alreadyplaced = 0 else: alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid]) valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid]) return valpos def index_2mat(matA,matB,val): """ Function to provide the coordinates of a given value inside two matrix simultaneously index_mat(matA,matB,val) matA= matrix with one set of values matB= matrix with the pother set of values val= couple of values to search >>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111]) [2 1 1] """ fname = 'index_2mat' matAshape = matA.shape matBshape = matB.shape for idv in range(len(matAshape)): if matAshape[idv] != matBshape[idv]: print errormsg print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \ 'and B:',matBshape[idv],'does not coincide!!' quit(-1) minA = np.min(matA) maxA = np.max(matA) minB = np.min(matB) maxB = np.max(matB) Ndims = len(matAshape) # valpos = np.ones((Ndims), dtype=int)*-1. valpos = np.zeros((Ndims), dtype=int) if val[0] < minA or val[0] > maxA: print warnmsg print ' ' + fname + ': first value:',val[0],'outside matA range',minA,',', \ maxA,'!!' return valpos if val[1] < minB or val[1] > maxB: print warnmsg print ' ' + fname + ': second value:',val[1],'outside matB range',minB,',', \ maxB,'!!' return valpos dist = np.zeros(tuple(matAshape), dtype=np.float) dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2) mindist = np.min(dist) if mindist != mindist: print ' ' + fname + ': wrong minimal distance',mindist,'!!' return valpos else: matlist = list(dist.flatten()) ifound = matlist.index(mindist) baseprevdims = np.zeros((Ndims), dtype=int) for dimid in range(Ndims): baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims]) if dimid == 0: alreadyplaced = 0 else: alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid]) valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid]) return valpos def index_mat(matA,val): """ Function to provide the coordinates of a given value inside a matrix index_mat(matA,val) matA= matrix with one set of values val= couple of values to search >>> index_mat(np.arange(27),22.3) 22 """ fname = 'index_mat' matAshape = matA.shape minA = np.min(matA) maxA = np.max(matA) Ndims = len(matAshape) # valpos = np.ones((Ndims), dtype=int)*-1. valpos = np.zeros((Ndims), dtype=int) if val < minA or val > maxA: print warnmsg print ' ' + fname + ': first value:',val,'outside matA range',minA,',', \ maxA,'!!' return valpos dist = np.zeros(tuple(matAshape), dtype=np.float) dist = (matA - np.float(val))**2 mindist = np.min(dist) if mindist != mindist: print ' ' + fname + ': wrong minimal distance',mindist,'!!' return valpos matlist = list(dist.flatten()) valpos = matlist.index(mindist) return valpos def index_mat_exact(mat,val): """ Function to provide the coordinates of a given exact value inside a matrix index_mat(mat,val) mat= matrix with values val= value to search >>> index_mat(np.arange(27).reshape(3,3,3),22) [2 1 1] """ fname = 'index_mat' matshape = mat.shape matlist = list(mat.flatten()) ifound = matlist.index(val) Ndims = len(matshape) valpos = np.zeros((Ndims), dtype=int) baseprevdims = np.zeros((Ndims), dtype=int) for dimid in range(Ndims): baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims]) if dimid == 0: alreadyplaced = 0 else: alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid]) valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid]) 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): """ Function to make coincident times for two different sets of CFtimes tvalB= time values B tunitA= time units times A to which we want to make coincidence tunitB= time units times B >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00', 'hours since 1949-12-01 00:00:00') [ 0. 3600. 7200. 10800. 14400. 18000. 21600. 25200. 28800. 32400.] >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00', 'hours since 1979-12-01 00:00:00') [ 9.46684800e+08 9.46688400e+08 9.46692000e+08 9.46695600e+08 9.46699200e+08 9.46702800e+08 9.46706400e+08 9.46710000e+08 9.46713600e+08 9.46717200e+08] """ import datetime as dt fname = 'coincident_CFtimes' trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3] trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3] tuA = tunitA.split(' ')[0] tuB = tunitB.split(' ')[0] if tuA != tuB: if tuA == 'microseconds': if tuB == 'microseconds': tB = tvalB*1. elif tuB == 'seconds': tB = tvalB*10.e6 elif tuB == 'minutes': tB = tvalB*60.*10.e6 elif tuB == 'hours': tB = tvalB*3600.*10.e6 elif tuB == 'days': tB = tvalB*3600.*24.*10.e6 else: print errormsg print ' ' + fname + ": combination of time untis: '" + tuA + \ "' & '" + tuB + "' not ready !!" quit(-1) elif tuA == 'seconds': if tuB == 'microseconds': tB = tvalB/10.e6 elif tuB == 'seconds': tB = tvalB*1. elif tuB == 'minutes': tB = tvalB*60. elif tuB == 'hours': tB = tvalB*3600. elif tuB == 'days': tB = tvalB*3600.*24. else: print errormsg print ' ' + fname + ": combination of time untis: '" + tuA + \ "' & '" + tuB + "' not ready !!" quit(-1) elif tuA == 'minutes': if tuB == 'microseconds': tB = tvalB/(60.*10.e6) elif tuB == 'seconds': tB = tvalB/60. elif tuB == 'minutes': tB = tvalB*1. elif tuB == 'hours': tB = tvalB*60. elif tuB == 'days': tB = tvalB*60.*24. else: print errormsg print ' ' + fname + ": combination of time untis: '" + tuA + \ "' & '" + tuB + "' not ready !!" quit(-1) elif tuA == 'hours': if tuB == 'microseconds': tB = tvalB/(3600.*10.e6) elif tuB == 'seconds': tB = tvalB/3600. elif tuB == 'minutes': tB = tvalB/60. elif tuB == 'hours': tB = tvalB*1. elif tuB == 'days': tB = tvalB*24. else: print errormsg print ' ' + fname + ": combination of time untis: '" + tuA + \ "' & '" + tuB + "' not ready !!" quit(-1) elif tuA == 'days': if tuB == 'microseconds': tB = tvalB/(24.*3600.*10.e6) elif tuB == 'seconds': tB = tvalB/(24.*3600.) elif tuB == 'minutes': tB = tvalB/(24.*60.) elif tuB == 'hours': tB = tvalB/24. elif tuB == 'days': tB = tvalB*1. else: print errormsg print ' ' + fname + ": combination of time untis: '" + tuA + \ "' & '" + tuB + "' not ready !!" quit(-1) else: print errormsg print ' ' + fname + ": time untis: '" + tuA + "' not ready !!" quit(-1) else: tB = tvalB*1. if trefA != trefB: trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S') trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S') difft = trefTB - trefTA diffv = difft.days*24.*3600.*10.e6 + difft.seconds*10.e6 + difft.microseconds print ' ' + fname + ': different reference refA:',trefTA,'refB',trefTB print ' difference:',difft,':',diffv,'microseconds' if tuA == 'microseconds': tB = tB + diffv elif tuA == 'seconds': tB = tB + diffv/10.e6 elif tuA == 'minutes': tB = tB + diffv/(60.*10.e6) elif tuA == 'hours': tB = tB + diffv/(3600.*10.e6) elif tuA == 'dayss': tB = tB + diffv/(24.*3600.*10.e6) else: print errormsg print ' ' + fname + ": time untis: '" + tuA + "' not ready !!" quit(-1) return tB def slice_variable(varobj, dimslice): """ Function to return a slice of a given variable according to values to its dimensions slice_variable(varobj, dimslice) varobj= object wit the variable dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension [value]: * [integer]: which value of the dimension * -1: all along the dimension * -9: last value of the dimension * [beg]@[end] slice from [beg] to [end] """ fname = 'slice_variable' if varobj == 'h': print fname + '_____________________________________________________________' print slice_variable.__doc__ quit() vardims = varobj.dimensions Ndimvar = len(vardims) Ndimcut = len(dimslice.split('|')) dimsl = dimslice.split('|') varvalsdim = [] dimnslice = [] for idd in range(Ndimvar): for idc in range(Ndimcut): dimcutn = dimsl[idc].split(':')[0] dimcutv = dimsl[idc].split(':')[1] if vardims[idd] == dimcutn: posfrac = dimcutv.find('@') if posfrac != -1: inifrac = int(dimcutv.split('@')[0]) endfrac = int(dimcutv.split('@')[1]) varvalsdim.append(slice(inifrac,endfrac)) dimnslice.append(vardims[idd]) else: if int(dimcutv) == -1: varvalsdim.append(slice(0,varobj.shape[idd])) dimnslice.append(vardims[idd]) elif int(dimcutv) == -9: varvalsdim.append(int(varobj.shape[idd])-1) else: varvalsdim.append(int(dimcutv)) break varvalues = varobj[tuple(varvalsdim)] return varvalues, dimnslice def func_compute_varNOcheck(ncobj, varn): """ Function to compute variables which are not originary in the file ncobj= netCDF object file varn = variable to compute: '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 'WRFz': height from WRF variables """ fname = 'compute_varNOcheck' if varn == 'WRFdens': # print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \ # 'DNW)/g ...' grav = 9.81 # Just we need in in absolute values: Size of the central grid cell ## dxval = ncobj.getncattr('DX') ## dyval = ncobj.getncattr('DY') ## mapfac = ncobj.variables['MAPFAC_M'][:] ## area = dxval*dyval*mapfac dimensions = ncobj.variables['MU'].dimensions mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:]) dnw = ncobj.variables['DNW'][:] varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \ dtype=np.float) levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float) for it in range(mu.shape[0]): for iz in range(dnw.shape[1]): levval.fill(np.abs(dnw[it,iz])) varNOcheck[it,iz,:,:] = levval varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav elif varn == 'WRFght': # print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...' varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:] dimensions = ncobj.variables['PH'].dimensions elif varn == 'WRFp': # print ' ' + fname + ': Retrieving pressure value from WRF as P + PB' varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:] dimensions = ncobj.variables['P'].dimensions elif varn == 'WRFrh': # print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\ # ' equation (T,P) ...' p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) qv = ncobj.variables['QVAPOR'][:] data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65)) data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1) varNOcheckv = qv/data2 dimensions = ncobj.variables['P'].dimensions elif varn == 'WRFt': # print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...' p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) dimensions = ncobj.variables['T'].dimensions elif varn == 'WRFz': grav = 9.81 # print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...' varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav dimensions = ncobj.variables['PH'].dimensions else: print erromsg print ' ' + fname + ": variable '" + varn + "' nor ready !!" quit(-1) return varNOcheck class compute_varNOcheck(object): """ Class to compute variables which are not originary in the file ncobj= netCDF object file varn = variable to compute: 'WRFdens': air density from WRF variables 'WRFght': geopotential height from WRF variables '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 'WRFws': wind speed from WRF variables 'WRFwd': wind direction from WRF variables 'WRFz': height from WRF variables """ fname = 'compute_varNOcheck' def __init__(self, ncobj, varn): if ncobj is None: self = None self.dimensions = None self.shape = None self.__values = None else: if varn == 'WRFdens': # print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \ # 'DNW)/g ...' grav = 9.81 # Just we need in in absolute values: Size of the central grid cell ## dxval = ncobj.getncattr('DX') ## dyval = ncobj.getncattr('DY') ## mapfac = ncobj.variables['MAPFAC_M'][:] ## area = dxval*dyval*mapfac dimensions = ncobj.variables['MU'].dimensions shape = ncobj.variables['MU'].shape mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:]) dnw = ncobj.variables['DNW'][:] varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \ dtype=np.float) levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float) for it in range(mu.shape[0]): for iz in range(dnw.shape[1]): levval.fill(np.abs(dnw[it,iz])) varNOcheck[it,iz,:,:] = levval varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav elif varn == 'WRFght': # print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...' varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:] dimensions = ncobj.variables['PH'].dimensions shape = ncobj.variables['PH'].shape elif varn == 'WRFp': # print ' ' + fname + ': Retrieving pressure value from WRF as P + PB' varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:] dimensions = ncobj.variables['P'].dimensions shape = ncobj.variables['P'].shape elif varn == 'WRFrh': # print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\ # ' equation (T,P) ...' p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) qv = ncobj.variables['QVAPOR'][:] data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65)) data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1) varNOcheckv = qv/data2 dimensions = ncobj.variables['P'].dimensions 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) ...' p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) dimensions = ncobj.variables['T'].dimensions shape = ncobj.variables['P'].shape elif varn == 'WRFtd': # print ' ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...' # tacking from: http://en.wikipedia.org/wiki/Dew_point p0=100000. p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] temp = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) qv = ncobj.variables['QVAPOR'][:] tk = temp - 273.15 data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65)) data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1) rh = qv/data2 pa = rh * data1/100. varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121)) dimensions = ncobj.variables['T'].dimensions shape = ncobj.variables['P'].shape elif varn == 'WRFws': # print ' ' + main + ': computing wind speed from WRF as SQRT(U**2 + V**2) ...' uwind = ncobj.variables['U'][:] vwind = ncobj.variables['V'][:] dx = uwind.shape[3] dy = vwind.shape[2] # de-staggering ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx]) va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:]) varNOcheckv = np.sqrt(ua*ua + va*va) dimensions = tuple(['Time','bottom_top','south_north','west_east']) shape = ua.shape elif varn == 'WRFwd': # print ' ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...' uwind = ncobj.variables['U'][:] vwind = ncobj.variables['V'][:] dx = uwind.shape[3] dy = vwind.shape[2] # de-staggering ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx]) va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:]) theta = np.arctan2(ua, va) dimensions = tuple(['Time','bottom_top','south_north','west_east']) shape = ua.shape varNOcheckv = 360.*(1. + theta/(2.*np.pi)) elif varn == 'WRFz': grav = 9.81 # print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...' varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav dimensions = ncobj.variables['PH'].dimensions shape = ncobj.variables['PH'].shape else: print errormsg print ' ' + fname + ": variable '" + varn + "' nor ready !!" quit(-1) self.dimensions = dimensions self.shape = shape self.__values = varNOcheckv def __getitem__(self,elem): return self.__values[elem] ####### ###### ##### #### ### ## # strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \ " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')" kindobs=['multi-points', 'single-station', 'trajectory'] strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \ "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\ "'trajectory': following a trajectory" simopers = ['sumc','subc','mulc','divc'] opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \ 'substract [constant] to variables values; mulc,[constant]: multipy by ' + \ '[constant] to variables values; divc,[constant]: divide by [constant] to ' + \ 'variables values' 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': 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 " + \ "each source ([DIM]='X','Y','Z','T'; None, no value)" vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " + \ "names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " + \ "no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")" varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " + \ "validate and if necessary operation and value operations: " + opersinf + \ "(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', 'minimum', 'maximum', 'mean', 'mean2', \ 'standard deviation'] parser = OptionParser() parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES") parser.add_option("-D", "--vardimensions", dest="vardims", help=vardimshelp, metavar="VALUES") parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, help=strkObs, metavar="FILE") parser.add_option("-l", "--stationLocation", dest="stloc", help="longitude, latitude and height of the station (only for 'single-station')", metavar="FILE") parser.add_option("-o", "--observation", dest="fobs", help="observations file to validate", metavar="FILE") parser.add_option("-s", "--simulation", dest="fsim", help="simulation file to validate", metavar="FILE") parser.add_option("-t", "--trajectoryfile", dest="trajf", help="file with grid points of the trajectory in the simulation grid ('simtrj')", metavar="FILE") parser.add_option("-v", "--variables", dest="vars", help=varshelp, metavar="VALUES") (opts, args) = parser.parse_args() ####### ####### ## MAIN ####### ofile='validation_sim.nc' if opts.dims is None: print errormsg print ' ' + main + ': No list of dimensions are provided!!' print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\ ' is needed' quit(-1) else: simdims = {} obsdims = {} print main +': couple of dimensions _______' dims = {} ds = opts.dims.split(',') for d in ds: dsecs = d.split('@') if len(dsecs) != 3: print errormsg print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!' print ' [DIM]@[dimnsim]@[dimnobs]' quit(-1) dims[dsecs[0]] = [dsecs[1], dsecs[2]] simdims[dsecs[0]] = dsecs[1] obsdims[dsecs[0]] = dsecs[2] print ' ',dsecs[0],':',dsecs[1],',',dsecs[2] if opts.vardims is None: print errormsg print ' ' + main + ': No list of variables with dimension values are provided!!' print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \ '[vardimtsim]@[vardimtobs] is needed' quit(-1) else: print main +': couple of variable dimensions _______' vardims = {} ds = opts.vardims.split(',') for d in ds: dsecs = d.split('@') if len(dsecs) != 3: print errormsg print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!' print ' [DIM]@[vardimnsim]@[vardimnobs]' quit(-1) vardims[dsecs[0]] = [dsecs[1], dsecs[2]] print ' ',dsecs[0],':',dsecs[1],',',dsecs[2] if opts.obskind is None: print errormsg print ' ' + main + ': No kind of observations provided !!' quit(-1) else: obskind = opts.obskind if obskind == 'single-station': if opts.stloc is None: print errormsg print ' ' + main + ': No station location provided !!' quit(-1) else: stationdesc = [np.float(opts.stloc.split(',')[0]), \ np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])] if opts.fobs is None: print errormsg print ' ' + main + ': No observations file is provided!!' quit(-1) else: if not os.path.isfile(opts.fobs): print errormsg print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!" quit(-1) if opts.fsim is None: print errormsg print ' ' + main + ': No simulation file is provided!!' quit(-1) else: if not os.path.isfile(opts.fsim): print errormsg print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!" quit(-1) if opts.vars is None: print errormsg print ' ' + main + ': No list of couples of variables is provided!!' print ' a ',' list of values [varsim]@[varobs],... is needed' quit(-1) else: valvars = [] vs = opts.vars.split(',') for v in vs: vsecs = v.split('@') if len(vsecs) < 2: print errormsg print ' ' + main + ': wrong number of values in:',vsecs, \ ' at least 2 are needed !!' print ' [varsim]@[varobs]@[[oper][val]]' quit(-1) if len(vsecs) > 2: if not searchInlist(simopers,vsecs[2]): print errormsg print main + ": operation on simulation values '" + vsecs[2] + \ "' not ready !!" quit(-1) valvars.append(vsecs) # Openning observations trajectory ## oobs = NetCDFFile(opts.fobs, 'r') valdimobs = {} for dn in dims: print dn,':',dims[dn] if dims[dn][1] != 'None': if not oobs.dimensions.has_key(dims[dn][1]): print errormsg print ' ' + main + ": observations file does not have dimension '" + \ dims[dn][1] + "' !!" quit(-1) if vardims[dn][1] != 'None': if not oobs.variables.has_key(vardims[dn][1]): print errormsg print ' ' + main + ": observations file does not have varibale " + \ "dimension '" + vardims[dn][1] + "' !!" quit(-1) valdimobs[dn] = oobs.variables[vardims[dn][1]][:] else: if dn == 'X': valdimobs[dn] = stationdesc[0] elif dn == 'Y': valdimobs[dn] = stationdesc[1] elif dn == 'Z': valdimobs[dn] = stationdesc[2] osim = NetCDFFile(opts.fsim, 'r') valdimsim = {} for dn in dims: if not osim.dimensions.has_key(dims[dn][0]): print errormsg 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 '" + opts.fsim + "' does not have " + \ "varibale dimension '" + vardims[dn][0] + "' !!" print ' it has variables:',osim.variables quit(-1) if searchInlist(varNOcheck,vardims[dn][0]): valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0]) else: valdimsim[dn] = osim.variables[vardims[dn][0]][:] # General characteristics dimtobs = valdimobs['T'].shape[0] dimtsim = valdimsim['T'].shape[0] print main +': observational time-steps:',dimtobs,'simulation:',dimtsim notfound = np.zeros((dimtobs), dtype=int) if obskind == 'multi-points': trajpos = np.zeros((2,dimt),dtype=int) for it in range(dimtobs): trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \ [valdimobs['X'][it],valdimobss['Y'][it]]) elif obskind == 'single-station': stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \ valdimobs['X']]) stationpos = np.zeros((2), dtype=int) iid = 0 for idn in osim.variables[vardims['X'][0]].dimensions: if idn == dims['X'][0]: stationpos[1] = stsimpos[iid] elif idn == dims['Y'][0]: stationpos[0] = stsimpos[iid] iid = iid + 1 print main + ': station point in simulation:', stationpos print ' station position:',valdimobs['X'],',',valdimobs['Y'] print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \ valdimsim['Y'][tuple(stsimpos)] elif obskind == 'trajectory': if opts.trajf is not None: if not os.path.isfile(opts.fsim): print errormsg print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!" quit(-1) else: otrjf = NetCDFFile(opts.fsim, 'r') trajpos[0,:] = otrjf.variables['obssimtrj'][0] trajpos[1,:] = otrjf.variables['obssimtrj'][1] otrjf.close() else: if dims.has_key('Z'): trajpos = np.zeros((3,dimtobs),dtype=int) for it in range(dimtobs): if np.mod(it*100./dimtobs,10.) == 0.: print ' trajectory done: ',it*100./dimtobs,'%' stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \ [valdimobs['Y'][it],valdimobs['X'][it]]) stationpos = np.zeros((2), dtype=int) iid = 0 for idn in osim.variables[vardims['X'][0]].dimensions: if idn == dims['X'][0]: stationpos[1] = stsimpos[iid] elif idn == dims['Y'][0]: stationpos[0] = stsimpos[iid] iid = iid + 1 if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1 trajpos[0,it] = stationpos[0] trajpos[1,it] = stationpos[1] # In the simulation 'Z' varies with time ... non-hydrostatic model! ;) # trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \ # stationpos[1]], valdimobs['Z'][it]) else: trajpos = np.zeros((2,dimtobs),dtype=int) for it in range(dimtobs): stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \ [valdimobs['Y'][it],valdimobss['X'][it]]) stationpos = np.zeros((2), dtype=int) iid = 0 for idn in osim.variables[vardims['X'][0]].dimensions: if idn == dims['X'][0]: stationpos[1] = stsimpos[iid] elif idn == dims['Y'][0]: stationpos[0] = stsimpos[iid] iid = iid + 1 if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1 trajpos[0,it] = stationspos[0] trajpos[1,it] = stationspos[1] print main + ': not found',np.sum(notfound),'points of the trajectory' # Getting times tobj = oobs.variables[vardims['T'][1]] obstunits = tobj.getncattr('units') 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 ## coindtvalues0 = [] for it in range(dimtsim): ot = 0 for ito in range(ot,dimtobs-1): if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \ simobstimes[it]: ot = ito tdist = simobstimes[it] - valdimobs['T'][ito] coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \ tdist]) coindtvalues = np.array(coindtvalues0, dtype=np.float) Ncoindt = len(coindtvalues[:,0]) print main + ': found',Ncoindt,'coincident times between simulation and observations' if Ncoindt == 0: print warnmsg print main + ': no coincident times found !!' print ' stopping it' quit(-1) # Validating ## onewnc = NetCDFFile(ofile, 'w') # Dimensions newdim = onewnc.createDimension('time',None) newdim = onewnc.createDimension('bnds',2) newdim = onewnc.createDimension('obstime',None) newdim = onewnc.createDimension('couple',2) newdim = onewnc.createDimension('StrLength',StringLength) 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',6) # Variable dimensions ## newvar = onewnc.createVariable('obstime','f8',('time')) basicvardef(newvar, 'obstime', 'time observations', obstunits ) set_attribute(newvar, 'calendar', 'standard') set_attribute(newvar, 'bounds', 'time_bnds') newvar[:] = coindtvalues[:,3] newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength')) basicvardef(newvar, 'couple', 'couples of values', '-') writing_str_nc(newvar, ['sim','obs'], StringLength) newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength')) 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': if dims.has_key('Z'): newdim = onewnc.createDimension('trj',3) else: newdim = onewnc.createDimension('trj',2) newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj')) basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-') newvar[:] = trajpos.transpose() if dims.has_key('Z'): newdim = onewnc.createDimension('simtrj',4) trjsim = np.zeros((4,Ncoindt), dtype=int) trjsimval = np.zeros((4,Ncoindt), dtype=np.float) else: newdim = onewnc.createDimension('simtrj',3) trjsim = np.zeros((3,Ncoindt), dtype=int) trjsimval = np.zeros((3,Ncoindt), dtype=np.float) Nvars = len(valvars) for ivar in range(Nvars): simobsvalues = [] varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1] print ' ' + varsimobs + '... .. .' if not oobs.variables.has_key(valvars[ivar][1]): print errormsg print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \ "' !!" quit(-1) if not osim.variables.has_key(valvars[ivar][0]): if not searchInlist(varNOcheck, valvars[ivar][0]): print errormsg print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \ "' !!" quit(-1) else: ovsim = compute_varNOcheck(osim, valvars[ivar][0]) else: ovsim = osim.variables[valvars[ivar][0]] for idn in ovsim.dimensions: if not searchInlist(simdims.values(),idn): print errormsg print ' ' + main + ": dimension '" + idn + "' of variable '" + \ valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!" quit(-1) 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), \ dtype = np.float) 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': for it in range(Ncoindt): slicev = dims['X'][0] + ':' + str(trajpos[2,it]) + '|' + \ dims['Y'][0]+ ':' + str(trajpos[1,it]) + '|' + \ dims['T'][0]+ ':' + str(coindtvalues[it][0]) slicevar, dimslice = slice_variable(ovsim, slicev) simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]]) 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]) slicevar, dimslice = slice_variable(ovsim, slicev) simobsSvalues[it,:,:] = slicevar elif obskind == 'single-station': for it in range(Ncoindt): 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) 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 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': if dims.has_key('Z'): for it in range(Ncoindt): ito = int(coindtvalues[it,1]) if notfound[ito] == 0: trajpos[2,ito] = index_mat(valdimsim['Z'][coindtvalues[it,0],:, \ trajpos[1,ito],trajpos[0,ito]], valdimobs['Z'][ito]) slicev = dims['X'][0]+':'+str(trajpos[0,ito]) + '|' + \ dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' + \ dims['Z'][0]+':'+str(trajpos[2,ito]) + '|' + \ dims['T'][0]+':'+str(int(coindtvalues[it,0])) slicevar, dimslice = slice_variable(ovsim, slicev) simobsvalues.append([ slicevar, ovobs[int(ito)]]) minx = np.max([trajpos[0,ito]-Ngrid,0]) maxx = np.min([trajpos[0,ito]+Ngrid+1,ovsim.shape[3]]) miny = np.max([trajpos[1,ito]-Ngrid,0]) maxy = np.min([trajpos[1,ito]+Ngrid+1,ovsim.shape[2]]) minz = np.max([trajpos[2,ito]-Ngrid,0]) maxz = np.min([trajpos[2,ito]+Ngrid+1,ovsim.shape[1]]) slicev = dims['X'][0] + ':' + str(minx) + '@' + str(maxx) + '|' +\ dims['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '|' + \ dims['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '|' + \ dims['T'][0] + ':' + str(int(coindtvalues[it,0])) slicevar, dimslice = slice_variable(ovsim, slicev) sliceS = [] sliceS.append(it) sliceS.append(slice(0,maxz-minz)) sliceS.append(slice(0,maxy-miny)) sliceS.append(slice(0,maxx-minx)) simobsSvalues[tuple(sliceS)] = slicevar if ivar == 0: trjsim[0,it] = trajpos[0,ito] trjsim[1,it] = trajpos[1,ito] trjsim[2,it] = trajpos[2,ito] trjsim[3,it] = coindtvalues[it,0] else: simobsvalues.append([fillValueF, fillValueF]) simobsSvalues[it,:,:,:]= np.ones((Ngrid*2+1,Ngrid*2+1,Ngrid*2+1),\ dtype = np.float)*fillValueF else: for it in range(Ncoindt): if notfound[it] == 0: ito = coindtvalues[it,1] slicev = dims['X'][0]+':'+str(trajpos[2,ito]) + '|' + \ dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' + \ dims['T'][0]+':'+str(coindtvalues[ito,0]) slicevar, dimslice = slice_variable(ovsim, slicev) simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]]) slicev = dims['X'][0] + ':' + str(trajpos[0,it]-Ngrid) + '@' + \ str(trajpos[0,it]+Ngrid) + '|' + dims['Y'][0] + ':' + \ str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + \ '|' + dims['T'][0] + ':' + str(coindtvalues[it,0]) slicevar, dimslice = slice_variable(ovsim, slicev) simobsSvalues[it,:,:] = slicevar else: simobsvalues.append([fillValue, fillValue]) simobsSvalues[it,:,:] = np.ones((Ngrid*2+1,Ngrid*2+1), \ dtype = np.float)*fillValueF print simobsvalues[varsimobs][:][it] arrayvals = np.array(simobsvalues) if len(valvars[ivar]) > 2: const=np.float(valvars[ivar][3]) if valvars[ivar][2] == 'sumc': simobsSvalues = simobsSvalues + const arrayvals[:,0] = arrayvals[:,0] + const elif valvars[ivar][2] == 'subc': simobsSvalues = simobsSvalues - const arrayvals[:,0] = arrayvals[:,0] - const elif valvars[ivar][2] == 'mulc': simobsSvalues = simobsSvalues * const arrayvals[:,0] = arrayvals[:,0] * const elif valvars[ivar][2] == 'divc': simobsSvalues = simobsSvalues / const arrayvals[:,0] = arrayvals[:,0] / const else: print errormsg print ' ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!" 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[:,1], fillValueF) obsmask2 = obsmask*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((9), dtype=np.float) diffvals = np.zeros((Ncoindt), dtype=np.float) diffvals = arrayvals[:,0] - obsmask simobsstats[0] = simstats[0] - obsstats[0] 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): aroundstats[0,it] = np.min(simobsSvalues[it,]) aroundstats[1,it] = np.max(simobsSvalues[it,]) aroundstats[2,it] = np.mean(simobsSvalues[it,]) aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,]) aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \ aroundstats[2,it]) # Statistics around obs values aroundostats = np.zeros((6,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] = obsmask.min() aroundostats[2,it] = obsmask.max() aroundostats[3,it] = obsmask.mean() aroundostats[4,it] = obsmask2.mean() aroundostats[5,it] = np.sqrt(aroundostats[4,it] - aroundostats[3,it]* \ aroundostats[3,it]) # sim Values to netCDF newvar = onewnc.createVariable(valvars[ivar][0], 'f', ('time'), \ fill_value=fillValueF) descvar = 'simulated: ' + valvars[ivar][0] basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units')) newvar[:] = arrayvals[:,0] # obs Values to netCDF newvar = onewnc.createVariable(valvars[ivar][1], 'f', ('time'), \ fill_value=fillValueF) descvar = 'observed: ' + valvars[ivar][1] basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units')) newvar[:] = arrayvals[:,1] # Around values if not onewnc.variables.has_key(valvars[ivar][0] + 'around'): if dims.has_key('Z'): if not onewnc.dimensions.has_key('zaround'): newdim = onewnc.createDimension('zaround',Ngrid*2+1) newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \ ('time','zaround','yaround','xaround'), fill_value=fillValueF) else: newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \ ('time','yaround','xaround'), fill_value=fillValueF) descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0] basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units')) newvar[:] = simobsSvalues # 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 (' + str(Ngrid) + ', ' + str(Ngrid) + \ ') simulated statisitcs: ' + valvars[ivar][0] basicvardef(newvar, valvars[ivar][0] + 'staround', descvar, \ ovobs.getncattr('units')) newvar[:] = aroundstats.transpose() if not searchInlist(onewnc.variables, 'time_bnds'): newvar = onewnc.createVariable('time_bnds','f8',('time','bnds')) basicvardef(newvar, 'time_bnds', 'time', obstunits ) set_attribute(newvar, 'calendar', 'standard') newvar[:] = simobsTtvalues # 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')) set_attribute(newvar, 'cell_methods', 'statistics') newvar[:] = aroundostats.transpose() onewnc.sync() newvar = onewnc.createVariable('simtrj','i',('time','simtrj')) basicvardef(newvar, 'simtrj', 'coordinates [X,Y,Z,T] of the coincident trajectory ' +\ 'in sim', obstunits) newvar[:] = trjsim.transpose() # Global attributes ## set_attribute(onewnc,'author_nc','Lluis Fita') set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \ 'LMD-Jussieu, UPMC, Paris') set_attribute(onewnc,'country_nc','France') set_attribute(onewnc,'script_nc',main) set_attribute(onewnc,'version_script',version) set_attribute(onewnc,'information', \ 'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html') set_attribute(onewnc,'simfile',opts.fsim) set_attribute(onewnc,'obsfile',opts.fobs) onewnc.sync() onewnc.close() print main + ": successfull writting of '" + ofile + "' !!"