Changeset 271 in lmdz_wrf


Ignore:
Timestamp:
Feb 25, 2015, 11:21:23 AM (10 years ago)
Author:
lfita
Message:

Adding 3D variables on the 'compute_tevolboxtraj'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var_tools.py

    r257 r271  
    97639763          [time] [xpoint] [ypoint]
    97649764        [Tbeg]: equivalent first time-step of the trajectory within the netCDF file
    9765         [lonname],[latname],[timename]: longitude, latitude and time variables names
     9765        [lonname],[latname],[zname],[timename]: longitude, latitude, z and time variables names
    97669766        [timekind]: kind of time
    97679767          'cf': cf-compilant
    … …  
    97709770        [circler]: radius in grid points of a centerd circle
    97719771      ncfile= netCDF file to use
    9772       varn= variable name
     9772      varn= ',' list of variables' name ('all', for all variables)
    97739773        EMERGENCY version, assuming 3D [time,lat,lon] variable !
    97749774    """
    … …  
    97889788    lonn = values.split(',')[1]
    97899789    latn = values.split(',')[2]
    9790     timn = values.split(',')[3]
     9790    zn = values.split(',')[3]
     9791    timn = values.split(',')[4]
    97919792    timekind = values.split(',')[4]
    97929793    boxs = int(values.split(',')[5])
    … …  
    98279828    timobj = objfile.variables[timn]
    98289829
     9830    if varn.find(',') != -1:
     9831        varns = varn.split(',')
     9832    else:
     9833        if varn == 'all':
     9834            varns = objfile.variables
     9835        else:
     9836            varns = varn
     9837
    98299838    if len(lonobj.shape) == 3:
    98309839        lonv = lonobj[0,:,:]
    … …  
    98449853          ncfile + '" !!!!!'
    98459854        quit(-1)
    9846 
    9847     varobj = objfile.variables[varn]
    98489855
    98499856# Selecting accordingly a trajectory
    … …  
    98609867    trajobj = open(trajfile,'r')
    98619868
     9869# Trajectory values/slices
     9870##
     9871    iline = 0
     9872    it = 0
     9873
     9874    xrangeslice = []
     9875    yrangeslice = []
     9876    xrangeslice2D = []
     9877    yrangeslice2D = []
     9878
     9879    cxrangeslice = []
     9880    cyrangeslice = []
     9881    cxrangeslice2D = []
     9882    cyrangeslice2D = []
     9883
     9884    cslicev = []
     9885    cslice2D = []
     9886    cslicevnoT = []
     9887
    98629888    gtrajvals = np.zeros((Ttraj,3), dtype=int)
    98639889    trajvals = np.zeros((Ttraj,3), dtype=np.float)
    9864     varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
    9865     lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
    9866     latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
    9867     rvarvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
    9868     rlonvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
    9869     rlatvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
    9870 
    9871     statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
    9872     rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
    9873 
    9874     iline = 0
    9875     it = 0
     9890
    98769891    for line in trajobj:
    9877         slicev = []
    9878         slice2D = []
    9879         slicevnoT = []
    98809892
    98819893## Skipping first line
    … …  
    99509962                    xend2D = boxs
    99519963
    9952                 slicev.append(slice(yinit, yend))
    9953                 slicev.append(slice(xinit, xend))
    9954 
    9955                 slicevnoT.append(slice(yinit, yend))
    9956                 slicevnoT.append(slice(xinit, xend))
    9957 
    9958                 slice2D.append(slice(yinit2D, yend2D))
    9959                 slice2D.append(slice(xinit2D, xend2D))
    9960 
    9961                 varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
    9962                 varvals[it,:,:] = varvalst
    9963 
    9964 # box stats values
    9965                 maskedvals = ma.masked_values (varvalst, fillValue)
    9966                 statvarvals[it,0] = varvalst[box2,box2]
    9967                 statvarvals[it,1] = maskedvals.min()
    9968                 statvarvals[it,2] = maskedvals.max()
    9969                 statvarvals[it,3] = maskedvals.mean()
    9970                 maskedvals2 = maskedvals*maskedvals
    9971                 statvarvals[it,4] = maskedvals2.mean()
    9972                 statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                      \
    9973                   statvarvals[it,3]*statvarvals[it,3])
    9974 
    9975                 varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
    9976                 lonvals[it,:,:] = varvalst
    9977                 varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
    9978                 latvals[it,:,:] = varvalst
    9979 
    9980             else:
    9981                 slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1))
    9982                 slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1))
    9983                 slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
    9984                 slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
    9985 
    9986                 slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 + 1))
    9987                 slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 + 1))
    9988 
    9989                 varvalst = varobj[tuple(slicev)]
    9990 # box values
    9991 
    9992                 varvals[it,:,:] = varvalst
    9993 #            print 'values at time t______'
    9994 #            print varvalst
    9995 
    9996 # box stats values
    9997                 statvarvals[it,0] = varvalst[box2,box2]
    9998                 statvarvals[it,1] = np.min(varvalst)
    9999                 statvarvals[it,2] = np.max(varvalst)
    10000                 statvarvals[it,3] = np.mean(varvalst)
    10001                 statvarvals[it,4] = np.mean(varvalst*varvalst)
    10002                 statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                      \
    10003                   statvarvals[it,3]*statvarvals[it,3])
    10004 
    10005                 varvalst = lonv[tuple(slicevnoT)]
    10006                 lonvals[it,:,:] = varvalst
    10007                 varvalst = latv[tuple(slicevnoT)]           
    10008                 latvals[it,:,:] = varvalst
     9964                yrangeslice.append([yinit, yend])
     9965                xrangeslice.append([xinit, xend])
     9966                yrangeslice2D.append([yinit2D, yend2D])
     9967                xrangeslice2D.append([xinit2D, xend2D])
    100099968
    100109969# circle values
    10011             slicev = []
    10012             slice2D = []
    10013             slicevnoT = []
    10014 
    10015             slicev.append(gtrajvals[it,0])
     9970            cslicev.append(gtrajvals[it,0])
    100169971            circdist = radius_dist(dimy, dimx, gtrajvals[it,2], gtrajvals[it,1])
    100179972
    100189973            if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1 \
    100199974              or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
    10020 
    10021                 rvarvalst = np.ones((Nrad*2+1, Nrad*2+1), dtype=np.float)*fillValue
    100229975
    100239976                if gtrajvals[it,2]-Nrad < 0:
    … …  
    1004910002                    xend2D = 2*Nrad+1
    1005010003
    10051                 slicev.append(slice(yinit, yend))
    10052                 slicev.append(slice(xinit, xend))
    10053 
    10054                 slicevnoT.append(slice(yinit, yend))
    10055                 slicevnoT.append(slice(xinit, xend))
    10056 
    10057                 slice2D.append(slice(yinit2D, yend2D))
    10058                 slice2D.append(slice(xinit2D, xend2D))
    10059 
    10060                 rvarvalst[tuple(slice2D)] = varobj[tuple(slicev)]
    10061                 rvarvalst[tuple(slice2D)] = np.where(circdist[tuple(slice2D)] >      \
    10062                   np.float(Nrad), fillValue, rvarvalst[tuple(slice2D)])
    10063 
    10064                 rvarvals[it,:,:] = rvarvalst
    10065 
    10066 # circle stats values
    10067                 maskedvals = ma.masked_values (rvarvalst, fillValue)
    10068                 rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
    10069                 rstatvarvals[it,1] = maskedvals.min()
    10070                 rstatvarvals[it,2] = maskedvals.max()
    10071                 rstatvarvals[it,3] = maskedvals.mean()
    10072                 maskedvals2 = maskedvals*maskedvals
    10073                 rstatvarvals[it,4] = maskedvals2.mean()
    10074                 rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                    \
    10075                   rstatvarvals[it,3]*rstatvarvals[it,3])
    10076 
    10077                 rvarvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
    10078                 rlonvals[it,:,:] = rvarvalst
    10079                 rvarvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
    10080                 rlatvals[it,:,:] = rvarvalst
    10081 
    10082             else:
    10083                 slicev.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1))
    10084                 slicev.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1))
    10085                 slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad+1))
    10086                 slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad+1))
    10087 
    10088                 slice2D.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1))
    10089                 slice2D.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1))
    10090 
    10091                 rvarvalst = varobj[tuple(slicev)]
    10092                 cdist = circdist[tuple(slicevnoT)]
    10093 # circle values
    10094                 rvarvalst = np.where(cdist > np.float(Nrad), fillValue, rvarvalst)
    10095 
    10096                 rvarvals[it,:,:] = rvarvalst
    10097 
    10098 # circle stats values
    10099                 maskedvals = ma.masked_values (rvarvalst, fillValue)
    10100                 rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
    10101                 rstatvarvals[it,1] = maskedvals.min()
    10102                 rstatvarvals[it,2] = maskedvals.max()
    10103                 rstatvarvals[it,3] = maskedvals.mean()
    10104                 maskedvals2 = maskedvals*maskedvals
    10105                 rstatvarvals[it,4] = maskedvals2.mean()
    10106                 rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                    \
    10107                   rstatvarvals[it,3]*rstatvarvals[it,3])
    10108 
    10109                 rvarvalst = lonv[tuple(slicevnoT)]
    10110                 rlonvals[it,:,:] = rvarvalst
    10111                 rvarvalst = latv[tuple(slicevnoT)]           
    10112                 rlatvals[it,:,:] = rvarvalst
    10113 
    10114             it = it + 1
    10115 #            print 'statistics:',rstatvarvals[it,:]
     10004                cyrangeslice.append([yinit, yend])
     10005                cxrangeslice.append([xinit, xend])
     10006                cyrangeslice2D.append([yinit2D, yend2D])
     10007                cxrangeslice2D.append([xinit2D, xend2D])
    1011610008
    1011710009        iline = iline + 1
    … …  
    1015310045    newvar[:] = latvals
    1015410046
     10047    statnames = ['minbox', 'maxbox', 'meanbox', 'mean2box', 'stdevbox']
     10048    vstlname = ['minimum value within', 'maximum value within', 'mean value within', \
     10049      'squared mean value within', 'standard deviation value within']
     10050    cstatnames = ['mincircle','maxcircle','meancircle','mean2circle','stdevcircle']
     10051
     10052# Getting values
     10053##
     10054    ivar = 0
     10055    for vn in varns:
     10056        varobj = objfile.variables[vn]
     10057        Nvardims = len(varobj.shape)
     10058
     10059        slicev = []
     10060        slice2D = []
     10061        slicevnoT = []
     10062        cslicev = []
     10063        cslice2D = []
     10064        cslicevnoT = []
     10065
     10066        if Nvardims == 4:
     10067# Too optimistic, but at this stage...
     10068            dimt = varobj.shape[0]
     10069            dimz = varobj.shape[1]
     10070 
     10071            if not objofile.dimensions.has_key('z'):
     10072                newdim = objofile.createDimension('z', dimt)
     10073
     10074                varzobj = objfile.variables[zn]
     10075                newvar = objofile.createVariable(zn, 'f8', ('z'))
     10076
     10077                if searchInlist(varzobj.ncattrs(),'standard_name'):
     10078                    vsname = varzobj.getncattr('standard_name')
     10079                else:
     10080                    vsname = variables_values(zn)[1]
     10081                if searchInlist(varzobj.ncattrs(),'long_name'):
     10082                    vlname = varzobj.getncattr('long_name')
     10083                else:
     10084                    vlname = variables_values(zn)[4].replace('|',' ')
     10085                if searchInlist(varzobj.ncattrs(),'units'):
     10086                    vunits = varzobj.getncattr('units')
     10087                else:
     10088                    vunits = variables_values(zn)[5].replace('|',' ')
     10089
     10090                newattr = basicvardef(newvar, vsname, vlname, vunits)
     10091                if len(varzobj.shape) == 1:
     10092                    newvar[:] = varzobj[:]
     10093                elif len(varzobj.shape) == 2:
     10094                    newvar[:] = varzobj[0,:]
     10095                elif len(varzobj.shape) == 3:
     10096                    newvar[:] = varzobj[0,0,:]
     10097
     10098            varvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
     10099            lonvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
     10100            latvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
     10101            rvarvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10102            rlonvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10103            rlatvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10104
     10105            statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
     10106            rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
     10107
     10108            for it in range(dimt):
     10109# box values
     10110                slicev.appen(slice(0:dimz))
     10111                slicev.appen(slice(yrangeslice[it]))
     10112                slicev.appen(slice(xrangeslice[it]))
     10113
     10114                slicevnoT.appen(slice(0:dimz))
     10115                slicevnoT.appen(slice(yrangeslice[it]))
     10116                slicevnoT.appen(slice(xrangeslice[it]))
     10117
     10118                slice2D.appen(slice(0:dimz))
     10119                slice2D.appen(slice(yrangeslice2D[it]))
     10120                slice2D.appen(slice(xrangeslice2D[it]))
     10121
     10122                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
     10123                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
     10124
     10125                    rvarvalst = np.ones((dimz,Nrad*2+1,Nrad*2+1),dtype=np.float)*    \
     10126                      fillValue
     10127
     10128                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
     10129                    varvals[it,:,:,:] = varvalst
     10130   
     10131# box stats values
     10132                    maskedvals = ma.masked_values (varvalst, fillValue)
     10133                    statvarvals[it,:,0] = varvalst[box2,box2]
     10134                    statvarvals[it,:,1] = maskedvals.min()
     10135                    statvarvals[it,:,2] = maskedvals.max()
     10136                    statvarvals[it,:,3] = maskedvals.mean()
     10137                    maskedvals2 = maskedvals*maskedvals
     10138                    statvarvals[it,:,4] = maskedvals2.mean()
     10139                    statvarvals[it,:,5] = np.sqrt(statvarvals[it,:,4] -              \
     10140                      statvarvals[it,:,3]*statvarvals[it,:,3])
     10141
     10142                    varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
     10143                    lonvals[it,:,:,:] = varvalst
     10144                    varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
     10145                    latvals[it,:,:,:] = varvalst
     10146
     10147                else:
     10148                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
     10149                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
     10150                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
     10151                      box2+1))
     10152                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
     10153                      box2+1))
     10154
     10155                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
     10156                      1))
     10157                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
     10158                      1))
     10159
     10160                    varvalst = varobj[tuple(slicev)]
     10161# box values
     10162
     10163                    varvals[it,:,:,:] = varvalst
     10164#                    print 'values at time t______'
     10165#                    print varvalst
     10166
     10167# box stats values
     10168                    statvarvals[it,:,0] = varvalst[box2,box2]
     10169                    statvarvals[it,:,1] = np.min(varvalst)
     10170                    statvarvals[it,:,2] = np.max(varvalst)
     10171                    statvarvals[it,:,3] = np.mean(varvalst)
     10172                    statvarvals[it,:,4] = np.mean(varvalst*varvalst)
     10173                    statvarvals[it,:,5] = np.sqrt(statvarvals[it,:,4] -              \
     10174                      statvarvals[it,:,3]*statvarvals[it,:,3])
     10175
     10176                    varvalst = lonv[tuple(slicevnoT)]
     10177                    lonvals[it,:,:] = varvalst
     10178                    varvalst = latv[tuple(slicevnoT)]           
     10179                    latvals[it,:,:] = varvalst
     10180
     10181# Circle values
     10182                cslicev.appen(slice(0:dimz))
     10183                cslicev.appen(slice(cyrangeslice[it]))
     10184                cslicev.appen(slice(cxrangeslice[it]))
     10185
     10186                cslicevnoT.appen(slice(0:dimz))
     10187                cslicevnoT.appen(slice(cyrangeslice[it]))
     10188                cslicevnoT.appen(slice(cxrangeslice[it]))
     10189
     10190                cslice2D.appen(slice(0:dimz))
     10191                cslice2D.appen(slice(cyrangeslice2D[it]))
     10192                cslice2D.appen(slice(cxrangeslice2D[it]))
     10193
     10194                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1   \
     10195                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
     10196
     10197                    rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
     10198                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
     10199                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslice2D)] >\
     10200                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
     10201
     10202                    rvarvals[it,:,:,:] = rvarvalst
     10203
     10204# circle stats values
     10205                    maskedvals = ma.masked_values (rvarvalst, fillValue)
     10206                    rstatvarvals[it,:,0] = rvarvalst[Nrad,Nrad]
     10207                    rstatvarvals[it,:,1] = maskedvals.min()
     10208                    rstatvarvals[it,:,2] = maskedvals.max()
     10209                    rstatvarvals[it,:,3] = maskedvals.mean()
     10210                    maskedvals2 = maskedvals*maskedvals
     10211                    rstatvarvals[it,:,4] = maskedvals2.mean()
     10212                    rstatvarvals[it,:,5] = np.sqrt(rstatvarvals[it,:,4] -            \
     10213                      rstatvarvals[it,:,3]*rstatvarvals[it,:,3])
     10214
     10215                    rvarvalst[tuple(cslice2D)] = lonv[tuple(cslicevnoT)]
     10216                    rlonvals[it,:,:,:] = rvarvalst
     10217                    rvarvalst[tuple(cslice2D)] = latv[tuple(cslicevnoT)]           
     10218                    rlatvals[it,:,:,:] = rvarvalst
     10219
     10220                else:
     10221                    slicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
     10222                    slicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
     10223                    slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
     10224                      Nrad+1))
     10225                    slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
     10226                      Nrad+1))
     10227
     10228                    slice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
     10229                    slice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
     10230
     10231                    rvarvalst = varobj[tuple(cslicev)]
     10232                    cdist = circdist[tuple(cslicevnoT)]
     10233# circle values
     10234                    rvarvalst = np.where(cdist > np.float(Nrad),fillValue,rvarvalst)
     10235
     10236                    rvarvals[it,:,:,:] = rvarvalst
     10237
     10238# circle stats values
     10239                    maskedvals = ma.masked_values (rvarvalst, fillValue)
     10240                    rstatvarvals[it,:,0] = rvarvalst[Nrad,Nrad]
     10241                    rstatvarvals[it,:,1] = maskedvals.min()
     10242                    rstatvarvals[it,:,2] = maskedvals.max()
     10243                    rstatvarvals[it,:,3] = maskedvals.mean()
     10244                    maskedvals2 = maskedvals*maskedvals
     10245                    rstatvarvals[it,:,4] = maskedvals2.mean()
     10246                    rstatvarvals[it,:,5] = np.sqrt(rstatvarvals[it,4] -              \
     10247                      rstatvarvals[it,:,3]*rstatvarvals[it,3])
     10248
     10249                    rvarvalst = lonv[tuple(cslicevnoT)]
     10250                    rlonvals[it,:,:,:] = rvarvalst
     10251                    rvarvalst = latv[tuple(cslicevnoT)]           
     10252                    rlatvals[it,:,:,:] = rvarvalst
     10253
     10254#            print 'statistics:',rstatvarvals[it,:]
     10255
    1015510256# variable box values
    10156     newvar = objofile.createVariable(varn + 'box', 'f4', ('time', 'y', 'x'),         \
    10157       fill_value=fillValue)
    10158     if searchInlist(varobj.ncattrs(),'standard_name'):
    10159         vsname = varobj.getncattr('standard_name')
    10160     else:
    10161         vsname = variables_values(varn)[1]
    10162     if searchInlist(varobj.ncattrs(),'long_name'):
    10163         vlname = varobj.getncattr('long_name')
    10164     else:
    10165         vlname = variables_values(varn)[4].replace('|',' ')
    10166     if searchInlist(varobj.ncattrs(),'units'):
    10167         vunits = varobj.getncattr('units')
    10168     else:
    10169         vunits = variables_values(varn)[5].replace('|',' ')
    10170 
    10171     newattr = basicvardef(newvar, vsname, vlname, vunits)
    10172     newattr = newvar.setncattr('projection','lon lat')
    10173     newvar[:] = varvals
     10257            newvar = objofile.createVariable(varn + 'box', 'f4', ('time','z','y','x'),\
     10258              fill_value=fillValue)
     10259            if searchInlist(varobj.ncattrs(),'standard_name'):
     10260                vsname = varobj.getncattr('standard_name')
     10261            else:
     10262                vsname = variables_values(varn)[1]
     10263            if searchInlist(varobj.ncattrs(),'long_name'):
     10264                vlname = varobj.getncattr('long_name')
     10265            else:
     10266                vlname = variables_values(varn)[4].replace('|',' ')
     10267            if searchInlist(varobj.ncattrs(),'units'):
     10268                vunits = varobj.getncattr('units')
     10269            else:
     10270                vunits = variables_values(varn)[5].replace('|',' ')
     10271
     10272            newattr = basicvardef(newvar, vsname, vlname, vunits)
     10273            newattr = newvar.setncattr('projection','lon lat')
     10274            newvar[:] = varvals
    1017410275
    1017510276# center of the trajectory
    10176     newvar = objofile.createVariable('trj_' + varn, 'f', ('time'),                  \
    10177       fill_value=fillValue)
    10178     newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the trajectory of '+\
    10179       varn, vunits)
    10180     newvar[:] = statvarvals[:,0]
     10277            newvar = objofile.createVariable('trj_' + varn, 'f', ('time','z'),       \
     10278              fill_value=fillValue)
     10279            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +      \
     10280              'trajectory of '+ varn, vunits)
     10281            newvar[:] = statvarvals[:,:,0]
    1018110282
    1018210283# variable box statistics
    10183     ist = 0
    10184 
    10185     statnames = ['minbox', 'maxbox', 'meanbox', 'mean2box', 'stdevbox']
    10186     vstlname = ['minimum value within', 'maximum value within',                      \
    10187      'mean value within', 'squared mean value within',                               \
    10188      'standard deviation value within']
    10189 
    10190     for statn in statnames:
    10191         newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),          \
    10192           fill_value=fillValue)
    10193         newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +          \
    10194           ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + varn, vunits)
    10195         newvar[:] = statvarvals[:,ist+1]
    10196 #        newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
    10197         ist = ist + 1
     10284            ist = 0
     10285            for statn in statnames:
     10286                newvar = objofile.createVariable(statn + '_' + varn,'f',('time','z'),\
     10287                  fill_value=fillValue)
     10288                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
     10289                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + varn, vunits)
     10290                newvar[:] = statvarvals[:,:,ist+1]
     10291#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
     10292                ist = ist + 1
    1019810293       
    10199 
    1020010294# variable circle values
    10201     newvar = objofile.createVariable(varn + 'circle', 'f4', ('time', 'yr', 'xr'),      \
    10202       fill_value=fillValue)
    10203     if searchInlist(varobj.ncattrs(),'standard_name'):
    10204         vsname = varobj.getncattr('standard_name')
    10205     else:
    10206         vsname = variables_values(varn)[1]
    10207     if searchInlist(varobj.ncattrs(),'long_name'):
    10208         vlname = varobj.getncattr('long_name')
    10209     else:
    10210         vlname = variables_values(varn)[4].replace('|',' ')
    10211     if searchInlist(varobj.ncattrs(),'units'):
    10212         vunits = varobj.getncattr('units')
    10213     else:
    10214         vunits = variables_values(varn)[5].replace('|',' ')
    10215 
    10216     newattr = basicvardef(newvar, vsname, vlname, vunits)
    10217     newattr = newvar.setncattr('projection','lon lat')
    10218     newvar[:] = rvarvals
     10295            newvar = objofile.createVariable(varn + 'circle','f4',('time','z','yr',  \
     10296              'xr'), fill_value=fillValue)
     10297            if searchInlist(varobj.ncattrs(),'standard_name'):
     10298                vsname = varobj.getncattr('standard_name')
     10299            else:
     10300                vsname = variables_values(varn)[1]
     10301            if searchInlist(varobj.ncattrs(),'long_name'):
     10302                vlname = varobj.getncattr('long_name')
     10303            else:
     10304                vlname = variables_values(varn)[4].replace('|',' ')
     10305            if searchInlist(varobj.ncattrs(),'units'):
     10306                vunits = varobj.getncattr('units')
     10307            else:
     10308                vunits = variables_values(varn)[5].replace('|',' ')
     10309
     10310            newattr = basicvardef(newvar, vsname, vlname, vunits)
     10311            newattr = newvar.setncattr('projection','lon lat')
     10312            newvar[:] = rvarvals
    1021910313
    1022010314# variable circle statistics
    10221     ist = 0
    10222     statnames = ['mincircle', 'maxcircle', 'meancircle', 'mean2circle', 'stdevcircle']
    10223 
    10224     for statn in statnames:
    10225         newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),          \
    10226           fill_value=fillValue)
    10227         newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +          \
    10228           ' the circle of radius ('+ str(circler)+') of ' + varn, vunits)
    10229         newvar[:] = rstatvarvals[:,ist+1]
    10230 #        newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
    10231         ist = ist + 1
     10315            ist = 0
     10316            for statn in cstatnames:
     10317                newvar = objofile.createVariable(statn + '_' + varn,'f',('time','z'),\
     10318                  fill_value=fillValue)
     10319                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
     10320                  ' the circle of radius ('+ str(circler)+') of ' + varn, vunits)
     10321                newvar[:] = rstatvarvals[:,:,ist+1]
     10322#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
     10323                ist = ist + 1
     10324
     10325        elif Nvardims == 3:
     10326# Too optimistic, but at this stage...
     10327            dimt = varobj.shape[0]
     10328 
     10329            varvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
     10330            lonvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
     10331            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
     10332            rvarvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10333            rlonvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10334            rlatvals = np.ones(tuple([Ttraj,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10335
     10336            statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
     10337            rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
     10338
     10339            for it in range(dimt):
     10340# box values
     10341                slicev.appen(slice(yrangeslice[it]))
     10342                slicev.appen(slice(xrangeslice[it]))
     10343
     10344                slicevnoT.appen(slice(yrangeslice[it]))
     10345                slicevnoT.appen(slice(xrangeslice[it]))
     10346
     10347                slice2D.appen(slice(yrangeslice2D[it]))
     10348                slice2D.appen(slice(xrangeslice2D[it]))
     10349
     10350                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
     10351                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
     10352
     10353                    rvarvalst = np.ones((Nrad*2+1, Nrad*2+1),dtype=np.float)*fillValue
     10354
     10355                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
     10356                    varvals[it,:,:] = varvalst
     10357   
     10358# box stats values
     10359                    maskedvals = ma.masked_values (varvalst, fillValue)
     10360                    statvarvals[it,0] = varvalst[box2,box2]
     10361                    statvarvals[it,1] = maskedvals.min()
     10362                    statvarvals[it,2] = maskedvals.max()
     10363                    statvarvals[it,3] = maskedvals.mean()
     10364                    maskedvals2 = maskedvals*maskedvals
     10365                    statvarvals[it,4] = maskedvals2.mean()
     10366                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                   \
     10367                      statvarvals[it,3]*statvarvals[it,3])
     10368
     10369                    varvalst[tuple(slice2D)] = lonv[tuple(slicevnoT)]
     10370                    lonvals[it,:,:] = varvalst
     10371                    varvalst[tuple(slice2D)] = latv[tuple(slicevnoT)]           
     10372                    latvals[it,:,:] = varvalst
     10373
     10374                else:
     10375                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
     10376                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
     10377                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
     10378                      box2+1))
     10379                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
     10380                      box2+1))
     10381
     10382                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
     10383                      1))
     10384                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
     10385                      1))
     10386
     10387                    varvalst = varobj[tuple(slicev)]
     10388# box values
     10389
     10390                    varvals[it,:,:] = varvalst
     10391#                    print 'values at time t______'
     10392#                    print varvalst
     10393
     10394# box stats values
     10395                    statvarvals[it,0] = varvalst[box2,box2]
     10396                    statvarvals[it,1] = np.min(varvalst)
     10397                    statvarvals[it,2] = np.max(varvalst)
     10398                    statvarvals[it,3] = np.mean(varvalst)
     10399                    statvarvals[it,4] = np.mean(varvalst*varvalst)
     10400                    statvarvals[it,5] = np.sqrt(statvarvals[it,4] -                  \
     10401                      statvarvals[it,3]*statvarvals[it,3])
     10402
     10403                    varvalst = lonv[tuple(slicevnoT)]
     10404                    lonvals[it,:,:] = varvalst
     10405                    varvalst = latv[tuple(slicevnoT)]           
     10406                    latvals[it,:,:] = varvalst
     10407
     10408# Circle values
     10409                cslicev.appen(slice(cyrangeslice[it]))
     10410                cslicev.appen(slice(cxrangeslice[it]))
     10411
     10412                cslicevnoT.appen(slice(cyrangeslice[it]))
     10413                cslicevnoT.appen(slice(cxrangeslice[it]))
     10414
     10415                cslice2D.appen(slice(cyrangeslice2D[it]))
     10416                cslice2D.appen(slice(cxrangeslice2D[it]))
     10417                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1   \
     10418                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
     10419
     10420                    rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
     10421                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
     10422                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslice2D)] >\
     10423                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
     10424
     10425                    rvarvals[it,:,:] = rvarvalst
     10426
     10427# circle stats values
     10428                    maskedvals = ma.masked_values (rvarvalst, fillValue)
     10429                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
     10430                    rstatvarvals[it,1] = maskedvals.min()
     10431                    rstatvarvals[it,2] = maskedvals.max()
     10432                    rstatvarvals[it,3] = maskedvals.mean()
     10433                    maskedvals2 = maskedvals*maskedvals
     10434                    rstatvarvals[it,4] = maskedvals2.mean()
     10435                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
     10436                      rstatvarvals[it,3]*rstatvarvals[it,3])
     10437
     10438                    rvarvalst[tuple(cslice2D)] = lonv[tuple(cslicevnoT)]
     10439                    rlonvals[it,:,:] = rvarvalst
     10440                    rvarvalst[tuple(cslice2D)] = latv[tuple(cslicevnoT)]           
     10441                    rlatvals[it,:,:] = rvarvalst
     10442
     10443                else:
     10444                    slicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
     10445                    slicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
     10446                    slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
     10447                      Nrad+1))
     10448                    slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
     10449                      Nrad+1))
     10450
     10451                    slice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
     10452                    slice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
     10453
     10454                    rvarvalst = varobj[tuple(cslicev)]
     10455                    cdist = circdist[tuple(cslicevnoT)]
     10456# circle values
     10457                    rvarvalst = np.where(cdist > np.float(Nrad),fillValue,rvarvalst)
     10458
     10459                    rvarvals[it,:,:] = rvarvalst
     10460
     10461# circle stats values
     10462                    maskedvals = ma.masked_values (rvarvalst, fillValue)
     10463                    rstatvarvals[it,0] = rvarvalst[Nrad,Nrad]
     10464                    rstatvarvals[it,1] = maskedvals.min()
     10465                    rstatvarvals[it,2] = maskedvals.max()
     10466                    rstatvarvals[it,3] = maskedvals.mean()
     10467                    maskedvals2 = maskedvals*maskedvals
     10468                    rstatvarvals[it,4] = maskedvals2.mean()
     10469                    rstatvarvals[it,5] = np.sqrt(rstatvarvals[it,4] -                \
     10470                      rstatvarvals[it,3]*rstatvarvals[it,3])
     10471
     10472                    rvarvalst = lonv[tuple(cslicevnoT)]
     10473                    rlonvals[it,:,:] = rvarvalst
     10474                    rvarvalst = latv[tuple(cslicevnoT)]           
     10475                    rlatvals[it,:,:] = rvarvalst
     10476
     10477#            print 'statistics:',rstatvarvals[it,:]
     10478
     10479# variable box values
     10480            newvar = objofile.createVariable(varn + 'box', 'f4', ('time', 'y', 'x'), \
     10481              fill_value=fillValue)
     10482            if searchInlist(varobj.ncattrs(),'standard_name'):
     10483                vsname = varobj.getncattr('standard_name')
     10484            else:
     10485                vsname = variables_values(varn)[1]
     10486            if searchInlist(varobj.ncattrs(),'long_name'):
     10487                vlname = varobj.getncattr('long_name')
     10488            else:
     10489                vlname = variables_values(varn)[4].replace('|',' ')
     10490            if searchInlist(varobj.ncattrs(),'units'):
     10491                vunits = varobj.getncattr('units')
     10492            else:
     10493                vunits = variables_values(varn)[5].replace('|',' ')
     10494
     10495            newattr = basicvardef(newvar, vsname, vlname, vunits)
     10496            newattr = newvar.setncattr('projection','lon lat')
     10497            newvar[:] = varvals
     10498
     10499# center of the trajectory
     10500            newvar = objofile.createVariable('trj_' + varn, 'f', ('time'),           \
     10501              fill_value=fillValue)
     10502            newattr = basicvardef(newvar, 'trj_' + vsname, 'value along the ' +      \
     10503              'trajectory of '+ varn, vunits)
     10504            newvar[:] = statvarvals[:,0]
     10505
     10506# variable box statistics
     10507            ist = 0
     10508            for statn in statnames:
     10509                newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),  \
     10510                  fill_value=fillValue)
     10511                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
     10512                  ' the box ('+str(boxs)+'x'+str(boxs)+') of ' + varn, vunits)
     10513                newvar[:] = statvarvals[:,ist+1]
     10514#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
     10515                ist = ist + 1
     10516       
     10517# variable circle values
     10518            newvar = objofile.createVariable(varn + 'circle','f4',('time','yr','xr'),\
     10519              fill_value=fillValue)
     10520            if searchInlist(varobj.ncattrs(),'standard_name'):
     10521                vsname = varobj.getncattr('standard_name')
     10522            else:
     10523                vsname = variables_values(varn)[1]
     10524            if searchInlist(varobj.ncattrs(),'long_name'):
     10525                vlname = varobj.getncattr('long_name')
     10526            else:
     10527                vlname = variables_values(varn)[4].replace('|',' ')
     10528            if searchInlist(varobj.ncattrs(),'units'):
     10529                vunits = varobj.getncattr('units')
     10530            else:
     10531                vunits = variables_values(varn)[5].replace('|',' ')
     10532
     10533            newattr = basicvardef(newvar, vsname, vlname, vunits)
     10534            newattr = newvar.setncattr('projection','lon lat')
     10535            newvar[:] = rvarvals
     10536
     10537# variable circle statistics
     10538            ist = 0
     10539            for statn in cstatnames:
     10540                newvar = objofile.createVariable(statn + '_' + varn, 'f', ('time'),  \
     10541                  fill_value=fillValue)
     10542                newattr = basicvardef(newvar, statn + '_' + vsname, vstlname[ist] +  \
     10543                  ' the circle of radius ('+ str(circler)+') of ' + varn, vunits)
     10544                newvar[:] = rstatvarvals[:,ist+1]
     10545#                newattr = set_attributek(newvar,'_FillValue',fillValue,'npfloat')
     10546                ist = ist + 1
    1023210547
    1023310548# global attributes
    … …  
    1025310568
    1025410569#compute_tevolboxtraj('h', 'wrfout*', 'PSFC')
    10255 #compute_tevolboxtraj('control/trajectory.dat@0,XLONG,XLAT,Times,wrf,3,3',            \
     10570#compute_tevolboxtraj('control/trajectory.dat@0,XLONG,XLAT,ZNU,Times,wrf,3,3',       \
    1025610571#  '../../superstorm/control/wrfout/wrfout_d01_2001-11-09_00:00:00', 'PSFC')
    1025710572#compute_tevolboxtraj('control/trajectory.dat@0,lon,lat,time,cf,5,5',                 \
    1025810573#  'control/wss.nc', 'wss')
     10574# in /home/lluis/etudes/WL_HyMeX/superstorm
     10575compute_tevolboxtraj('control/trajectory.dat@0,XLONG,XLAT,ZNU,Times,wrf,3,3',       \
     10576  'control/wrfout/wrfout_d01_2001-11-09_00:00:00', 'PSFC')
    1025910577
    1026010578def numVector_String(vec,char):
Note: See TracChangeset for help on using the changeset viewer.