Changeset 283 in lmdz_wrf


Ignore:
Timestamp:
Feb 25, 2015, 2:40:29 PM (10 years ago)
Author:
lfita
Message:

Fixing variable 4dims

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/tools/nc_var_tools.py

    r282 r283  
    1040410404            trajvals[it,2] = latv[gtrajvals[it,2],gtrajvals[it,1]]
    1040510405
    10406 #            print iline, it,'time:',trajvals[it,0],'lon lat:', trajvals[it,1], trajvals[it,2]
     10406#            print iline, it,'time:',trajvals[it,0],'lon lat:', trajvals[it,1],       \
     10407#              trajvals[it,2]
    1040710408
    1040810409# Assuming time as the first dimension in a fortran like way (t=0 as 1)
     
    1041110412
    1041210413# box values
    10413             if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1 \
     10414            if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1       \
    1041410415              or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
    1041510416
     
    1049510496                cyrangeslice.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1])
    1049610497                cxrangeslice.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1])
    10497                 cyrangeslice2D.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad + 1])
    10498                 cxrangeslice2D.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad + 1])
     10498                cyrangeslice2D.append([gtrajvals[it,2]-Nrad, gtrajvals[it,2]+Nrad+1])
     10499                cxrangeslice2D.append([gtrajvals[it,1]-Nrad, gtrajvals[it,1]+Nrad+1])
    1049910500            it = it + 1
    1050010501            iline = iline + 1
     
    1053910540            dimz = varobj.shape[1]
    1054010541 
    10541             if not objofile.dimensions.has_key('z'):
    10542                 newdim = objofile.createDimension('z', dimz)
    10543 
    10544                 varzobj = objfile.variables[zn]
    10545                 newvar = objofile.createVariable(zn, 'f8', ('z'))
    10546 
    10547                 if searchInlist(varzobj.ncattrs(),'standard_name'):
    10548                     vsname = varzobj.getncattr('standard_name')
    10549                 else:
    10550                     vsname = variables_values(zn)[1]
    10551                 if searchInlist(varzobj.ncattrs(),'long_name'):
    10552                     vlname = varzobj.getncattr('long_name')
    10553                 else:
    10554                     vlname = variables_values(zn)[4].replace('|',' ')
    10555                 if searchInlist(varzobj.ncattrs(),'units'):
    10556                     vunits = varzobj.getncattr('units')
    10557                 else:
    10558                     vunits = variables_values(zn)[5].replace('|',' ')
    10559 
    10560                 newattr = basicvardef(newvar, vsname, vlname, vunits)
    10561                 if len(varzobj.shape) == 1:
    10562                     newvar[:] = varzobj[:]
    10563                 elif len(varzobj.shape) == 2:
    10564                     newvar[:] = varzobj[0,:]
    10565                 elif len(varzobj.shape) == 3:
    10566                     newvar[:] = varzobj[0,0,:]
    10567 
    1056810542            varvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
     10543            lonvals = np.ones(tuple([Ttraj,dimz,boxs,boxs]), dtype=np.float)
     10544            latvals = np.ones(tuple([Ttraj,boxs,boxs]), dtype=np.float)
    1056910545            rvarvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
    10570 
    10571             statvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
    10572             rstatvarvals = np.ones(tuple([Ttraj,6]), dtype=np.float)
    10573 
    10574             for it in range(dimt):
    10575 
     10546            rlonvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10547            rlatvals = np.ones(tuple([Ttraj,dimz,Nrad*2+1,Nrad*2+1]), dtype=np.float)
     10548
     10549            statvarvals = np.ones(tuple([Ttraj,dimz,6]), dtype=np.float)
     10550            rstatvarvals = np.ones(tuple([Ttraj,dimz,6]), dtype=np.float)
     10551
     10552            for it0 in range(Ttraj):
     10553                it = Tbeg + it0
     10554
     10555                slicev = []
     10556                slice2D = []
     10557                slicevnoT = []
     10558                cslicev = []
     10559                cslice2D = []
     10560                cslicevnoT = []
     10561
     10562                slicev.append(gtrajvals[it,0])
    1057610563                if gtrajvals[it,2]-box2 < 0 or gtrajvals[it,2]+box2 + 1 > dimy + 1   \
    1057710564                  or gtrajvals[it,1]-box2 < 0 or gtrajvals[it,1]+box2 + 1 > dimx + 1:
    10578 
    10579                     varvalst = np.ones((dimz,Nrad*2+1,Nrad*2+1),dtype=np.float)*    \
     10565# box values
     10566                    slicev.append(slice(0,dimz)
     10567                    slicev.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
     10568                    slicev.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
     10569
     10570                    slicevnoT.append(slice(0,dimz)
     10571                    slicevnoT.append(slice(yrangeslice[it][0],yrangeslice[it][1]))
     10572                    slicevnoT.append(slice(xrangeslice[it][0],xrangeslice[it][1]))
     10573
     10574                    slice2D.append(slice(0,dimz)
     10575                    slice2D.append(slice(0,yrangeslice[it][1]-yrangeslice[it][0]))
     10576                    slice2D.append(slice(0,xrangeslice[it][1]-xrangeslice[it][0]))
     10577
     10578                    rvarvalst = np.ones((dimz, Nrad*2+1, Nrad*2+1),dtype=np.float)*  \
    1058010579                      fillValue
    10581 
    10582 # box values
    10583                     slicev.append(slice(0,dimz))
    10584                     slicev.append(slice(yrangeslice[it]))
    10585                     slicev.append(slice(xrangeslice[it]))
    10586 
    10587                     slicevnoT.append(slice(0,dimz))
    10588                     slicevnoT.append(slice(yrangeslice[it]))
    10589                     slicevnoT.append(slice(xrangeslice[it]))
    10590 
    10591                     slice2D.append(slice(0,dimz))
    10592                     slice2D.append(slice(yrangeslice2D[it]))
    10593                     slice2D.append(slice(xrangeslice2D[it]))
    1059410580
    1059510581                    varvalst[tuple(slice2D)] = varobj[tuple(slicev)]
     
    1060810594
    1060910595                else:
     10596                    slicev.append(slice(0,dimz)
    1061010597                    slicev.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2+1))
    1061110598                    slicev.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2+1))
     10599                    slicevnoT.append(slice(0,dimz)
    1061210600                    slicevnoT.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+    \
    1061310601                      box2+1))
    1061410602                    slicevnoT.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+    \
    1061510603                      box2+1))
    10616 
     10604                    slice2D.append(slice(0,dimz)
    1061710605                    slice2D.append(slice(gtrajvals[it,2]-box2, gtrajvals[it,2]+box2 +\
    1061810606                      1))
    1061910607                    slice2D.append(slice(gtrajvals[it,1]-box2, gtrajvals[it,1]+box2 +\
    1062010608                      1))
    10621                     print 'Lluis',it,'shapes varvalst:',varvals.shape,'varobj:',slicev
     10609
    1062210610                    varvalst = varobj[tuple(slicev)]
    1062310611# box values
     
    1063710625
    1063810626# Circle values
    10639                 cslicev.append(slice(0,dimz))
    10640                 cslicev.append(slice(cyrangeslice[it]))
    10641                 cslicev.append(slice(cxrangeslice[it]))
    10642 
    10643                 cslicevnoT.append(slice(0,dimz))
    10644                 cslicevnoT.append(slice(cyrangeslice[it]))
    10645                 cslicevnoT.append(slice(cxrangeslice[it]))
    10646 
    10647                 cslice2D.append(slice(0,dimz))
    10648                 cslice2D.append(slice(cyrangeslice2D[it]))
    10649                 cslice2D.append(slice(cxrangeslice2D[it]))
    10650 
     10627                cslicev.append(gtrajvals[it,0])
    1065110628                if gtrajvals[it,2]-Nrad < 0 or gtrajvals[it,2]+Nrad + 1 > dimy + 1   \
    1065210629                  or gtrajvals[it,1]-Nrad < 0 or gtrajvals[it,1]+Nrad + 1 > dimx + 1:
    10653 
    10654                     rvarvalst = np.ones((Nrad*2+1,Nrad*2+1),dtype=np.float)*fillValue
     10630                    cslicev.append(slice(0,dimz)
     10631                    cslicev.append(slice(cyrangeslice[it][0],cyrangeslice[it][1]))
     10632                    cslicev.append(slice(cxrangeslice[it][0],cxrangeslice[it][1]))
     10633
     10634                    cslicevnoT.append(slice(0,dimz)
     10635                    cslicevnoT.append(slice(cyrangeslice[it][0],cyrangeslice[it][1]))
     10636                    cslicevnoT.append(slice(cxrangeslice[it][0],cxrangeslice[it][1]))
     10637
     10638                    cslice2D.append(slice(0,dimz)
     10639                    cslice2D.append(slice(0,cyrangeslice[it][1]-cyrangeslice[it][0]))
     10640                    cslice2D.append(slice(0,cxrangeslice[it][1]-cxrangeslice[it][0]))
     10641
     10642                    rvarvalst = np.ones((dimz,Nrad*2+1,Nrad*2+1),dtype=np.float)*    \
     10643                      fillValue
    1065510644                    rvarvalst[tuple(cslice2D)] = varobj[tuple(cslicev)]
    10656                     rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslice2D)] >\
     10645                    rvarvalst[tuple(cslice2D)] = np.where(circdist[tuple(cslice2D)]>\
    1065710646                      np.float(Nrad), fillValue, rvarvalst[tuple(cslice2D)])
    1065810647
     
    1067110660
    1067210661                else:
    10673                     slicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
    10674                     slicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
    10675                     slicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+    \
     10662                    cslicev.append(0,dimz)
     10663                    cslicev.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
     10664                    cslicev.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
     10665                    cslicevnoT.append(0,dimz)
     10666                    cslicevnoT.append(slice(gtrajvals[it,2]-Nrad, gtrajvals[it,2]+   \
    1067610667                      Nrad+1))
    10677                     slicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+    \
     10668                    cslicevnoT.append(slice(gtrajvals[it,1]-Nrad, gtrajvals[it,1]+   \
    1067810669                      Nrad+1))
    10679 
    10680                     slice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+1))
    10681                     slice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+1))
     10670                    cslice2D.append(0,dimz)
     10671                    cslice2D.append(slice(gtrajvals[it,2]-Nrad,gtrajvals[it,2]+Nrad+ \
     10672                      1))
     10673                    cslice2D.append(slice(gtrajvals[it,1]-Nrad,gtrajvals[it,1]+Nrad+ \
     10674                      1))
    1068210675
    1068310676                    rvarvalst = varobj[tuple(cslicev)]
     
    1069610689                    maskedvals2 = maskedvals*maskedvals
    1069710690                    rstatvarvals[it,:,4] = maskedvals2.mean()
    10698                     rstatvarvals[it,:,5] = np.sqrt(rstatvarvals[it,4] -              \
    10699                       rstatvarvals[it,:,3]*rstatvarvals[it,3])
     10691                    rstatvarvals[it,:,5] = np.sqrt(rstatvarvals[it,:,4] -            \
     10692                      rstatvarvals[it,:,3]*rstatvarvals[it,:,3])
    1070010693
    1070110694#            print 'statistics:',rstatvarvals[it,:]
Note: See TracChangeset for help on using the changeset viewer.