Changeset 2314 in lmdz_wrf for trunk/tools


Ignore:
Timestamp:
Feb 4, 2019, 5:47:35 PM (6 years ago)
Author:
lfita
Message:

Getting there with the 1D variables and the 3 un bounded variables!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var_tools.py

    r2313 r2314  
    244244        for idd in dimlist: dv.append(dicdim[idd])
    245245        self.values = values[:]
    246         self.dimensions = tuple(list(dicdim.keys()))
     246        self.dimensions = tuple(dimlist)
    247247        self.shape = tuple(dv)
    248248        self.name = name
    … …  
    2773027730        covered by each grid (defined as polygon) within the slice
    2773127731      values=[slicevarsinf]@[dimvars]@[sliceremovedim]@[slicebndsdim]@
    27732           [slicebndsvar]@[slicestatsdim]
     27732          [slicebndsvar]@[slicespacedim]@[slicestatsdim]
    2773327733        [slicevarsinf] ';' separated list of [varn1],[slcvar1];...;[varnN],[slcvarN]
    2773427734          varn[i]: name of the variable i
    … …  
    2776027760          bounds variables (as ';' list when more than 1) from the slicing variables
    2776127761          ('None' for any)
     27762        [slicespacedim]: ',' [ydimn],[xdimn] ordered list of dimensions from which 
     27763          construct the 2D horizontal space shape of the slices. When sliced variables
     27764          do not have them (e.g. 1D) they are expanded to them
    2776227765        [slicestatsdim]: ',' list of name of dimensions to do not use for computing
    2776327766          statistics of variables at each slice ('any' for not leaving any
    … …  
    2778227785
    2778327786    expectargs = '[slicevarsinf]@[dimvars]@[sliceremovedim]@[slicebndsdim]@' +       \
    27784       '[slicebndsvar]@[slicestatsdim]'
     27787      '[slicebndsvar]@[slicespacedim]@[slicestatsdim]'
    2778527788    gen.check_arguments(fname, values, expectargs, '@')
    2778627789
    … …  
    2779027793    slicebndsdim0 = values.split('@')[3]
    2779127794    slicebndsvar0 = values.split('@')[4]
    27792     slicestatsdim = values.split('@')[5].split(',')
     27795    slicespacedim = values.split('@')[5].split(',')
     27796    slicestatsdim = values.split('@')[6].split(',')
    2779327797
    2779427798    varvalues = varvalues0.split(';')
    … …  
    2780927813            slicebndsdim[dn] = bndn
    2781027814    else:
    27811         sliceremovdim = None
     27815        slicebndsdim = None
    2781227816
    2781327817    # Variables with boundaries from slicing variables
    … …  
    2783727841    onc = NetCDFFile(ncfile, 'r')
    2783827842
     27843    # Space shape
     27844    for dimn in slicespacedim:
     27845        if not gen.searchInlist(onc.dimensions,dimn):
     27846            print errormsg
     27847            print '  ' +fname+ ": file '" + ncfile + "' does not have space-slice "+ \
     27848              "dimension: '" + dn + "' !!"
     27849            dimns = list(onc.dimensions)
     27850            dimns.sort()
     27851            print '    available ones:', dimns
     27852            quit(-1)
     27853    sdimy = len(onc.dimensions[slicespacedim[0]])
     27854    sdimx = len(onc.dimensions[slicespacedim[1]])
     27855    print '    slicespacedim:', slicespacedim, 'sdimy, sdimx:', sdimy, sdimx
     27856
     27857    olon1D = onc.variables[gen.dictionary_key(dimvars, slicespacedim[1])]
     27858    olat1D = onc.variables[gen.dictionary_key(dimvars, slicespacedim[0])]
     27859    # Re-shaping them in case of 1D dimensions...
     27860    if len(olon1D.shape) == 1 and len(olat1D.shape) == 1:
     27861        print warnmsg
     27862        print '  ' + fname + ": 1D space dimensions !!"
     27863        print '    creation of a 2D version of them '
     27864        lon2D, lat2D = gen.lonlat2D(olon1D[:], olat1D[:])
     27865
     27866        # Boundaries
     27867        olonbnds1D = onc.variables[slicebndsdim[slicespacedim[1]]]
     27868        olatbnds1D = onc.variables[slicebndsdim[slicespacedim[0]]]
     27869        newlonbndsdim = ['bnds'] + slicespacedim[:]
     27870        newlatbndsdim = ['bnds'] + slicespacedim[:]
     27871        newlonbndsshape = [4,sdimy,sdimx]
     27872        newlatbndsshape = [4,sdimy,sdimx]
     27873        newlonbndsv = np.zeros((4,sdimy,sdimx), dtype=np.float)
     27874        newlatbndsv = np.zeros((4,sdimy,sdimx), dtype=np.float)
     27875        print 'Lluis shapes olonbnds1D:', olonbnds1D.shape
     27876        for j in range(sdimy):
     27877            for i in range(sdimx):
     27878                newlonbndsv[0,j,i] = olonbnds1D[i,0]
     27879                newlatbndsv[0,j,i] = olatbnds1D[j,1]
     27880                newlonbndsv[1,j,i] = olonbnds1D[i,1]
     27881                newlatbndsv[1,j,i] = olatbnds1D[j,1]
     27882                newlonbndsv[2,j,i] = olonbnds1D[i,1]
     27883                newlatbndsv[2,j,i] = olatbnds1D[j,0]
     27884                newlonbndsv[3,j,i] = olonbnds1D[i,0]
     27885                newlatbndsv[3,j,i] = olatbnds1D[j,0]
     27886
     27887        newvdimvs = {slicespacedim[1]: sdimx, slicespacedim[0]: sdimy}
     27888        newbdimvs = {slicespacedim[1]: sdimx, slicespacedim[0]: sdimy, 'bnds': 4}
     27889        newbdims = ['bnds'] + slicespacedim
     27890        olon2D = genericNCvariable_Dict(newvdimvs, slicespacedim, 'lon', 'lon',      \
     27891          'Longitude', 'degrees_east', lon2D)
     27892        olat2D = genericNCvariable_Dict(newvdimvs, slicespacedim, 'lat', 'lat',      \
     27893          'Latitude', 'degrees_north', lat2D)
     27894        olonbnds2D = genericNCvariable_Dict(newbdimvs, newbdims, slicespacedim[1] +  \
     27895          '_bnds', slicespacedim[1]+'lon_bnds', 'Longitude boundaries',              \
     27896          'degrees_east', newlonbndsv)
     27897        olatbnds2D = genericNCvariable_Dict(newbdimvs, newbdims, slicespacedim[0] +  \
     27898          '_bnds', slicespacedim[0]+'_bnds', 'Latitude boundaries','degrees_north',  \
     27899          newlatbndsv)
     27900
     27901        print "    new shape for space variable-dimensions '", olon2D.dimensions,    \
     27902          "' :", olon2D.shape, 'related boundaries:', olonbnds2D.shape
     27903
     27904        new2Dvars = {slicespacedim[1]: olon2D, slicespacedim[0]: olat2D,             \
     27905          slicespacedim[1]+'_bnds': olonbnds2D, slicespacedim[0]+'_bnds': olatbnds2D}
     27906        slicebndsvar[slicespacedim[1]] = [slicespacedim[1]+'_bnds',                  \
     27907          slicespacedim[0]+'_bnds']
     27908        slicebndsvar[slicespacedim[0]] = [slicespacedim[1]+'_bnds',                  \
     27909          slicespacedim[0]+'_bnds']
     27910   
    2783927911    # Varibales to slice
    2784027912    if variable == 'all':
    … …  
    2787327945            quit(-1)
    2787427946
    27875         ovar = onc.variables[varn]
     27947        if new2Dvars.has_key(varn):
     27948            print infmsg
     27949            print "  replacing 1D varibale '" + varn + "' by its 2D counterpart "
     27950            ovar = new2Dvars[varn]
     27951        else:
     27952            ovar = onc.variables[varn]
    2787627953        vu = ovar.units
    2787727954        dimnv = list(ovar.dimensions)
    … …  
    2799528072    vardimbndslice = {}
    2799628073
     28074    print 'Slicing following areas ...'
    2799728075    newslcvarns = []
    2799828076    for varn in slcvarns:
    27999         ovar = onc.variables[varn]
     28077        if gen.searchInlist(new2Dvars, varn):
     28078            print infmsg
     28079            print "  replacing 1D varibale '" + varn + "' by its 2D counterpart "
     28080            ovar = new2Dvars[varn]
     28081        else:
     28082            ovar = onc.variables[varn]
    2800028083        vu = ovar.units
    2800128084        dimnv = list(ovar.dimensions)
    … …  
    2803728120                    quit(-1)
    2803828121
    28039                 ovarbnds = onc.variables[vn]
     28122                if gen.searchInlist(new2Dvars, vn):
     28123                    print infmsg
     28124                    print "  replacing 1D bounded varibale '" + vn + "' by its 2D "+ \
     28125                      "counterpart "
     28126                    ovarbnds = new2Dvars[vn]
     28127                else:
     28128                    ovarbnds = onc.variables[vn]
    2804028129
    2804128130                # Removing undesired dimensions from bounds slicing variable
    … …  
    2805128140                dimsbnds[dnn] = rmshape
    2805228141
    28053             dn = '_'.join(varbnds)
    28054             if not gen.searchInlist(slcvarns, dn) and not                            \
    28055               gen.searchInlist(newslcvarns, dn) and not vardimbndslice.has_key(dn):
     28142            # Slicing variable name after its dimensions...
     28143            #dn = '_'.join(varbnds)
     28144            # Slicing variable name after its values ...
     28145            dn = varn + ''
     28146            if not gen.searchInlist(newslcvarns, dn) and                             \
     28147              not vardimbndslice.has_key(dn):
    2805628148                if len(dv.shape) == 1:
    2805728149                    print infmsg
    28058                     print '    ' + fname + ": slicing dimension '" + dn + "' ... "
     28150                    print '    ' + fname + ": slicing 1D dimension '" + dn + "' ... "
    2805928151
    2806028152                    # Checking for the existence of slicing information for the 
    … …  
    2811828210       
    2811928211                elif len(dv.shape) == 2:
    28120                     print "  2D slicing  bounded variable: '" + varn + "': ", dv.shape
     28212                    print infmsg
     28213                    print '    ' + fname + ": slicing 2D dimension '" + dn + "' ... "
    2812128214                    if len(rmdims) == 3:
    2812228215                        xdim = rmdims[2]
    … …  
    2812628219                        xdim = ovar.dimensions[rankv-1]
    2812728220                        ydim = ovar.dimensions[rankv-2]
     28221                    print "  2D slicing  bounded variable: '" + varn + "' rmdims:", rmdims
    2812828222                    getdims = dimdbnds[xdim]
    2812928223
    2813028224                    # ref values
     28225                    dref  = []
     28226                    if slicesinf.has_key(dimvars[ydim]):
     28227                        varslcv = slicesinf[dimvars[ydim]]
     28228                        dyref = varslcv[0]
     28229                        reflat1D = varslcv[3]
     28230                        refblat1D = varslcv[4]
     28231                        dref.append(varn)
     28232                    else:
     28233                        print infmsg
     28234                        print '  ' + fname + ": 2D slicing varibale '" + varn +      \
     28235                          "' with y-dimension '" + dimvars[ydim] + "' without " +    \
     28236                          "related slicing-variable !!"
     28237                        oyvar = onc.variables[dimvars[ydim]]
     28238
     28239                        # Removing undesired dimensions from bounds slicing variable
     28240                        if sliceremovedim is not None:
     28241                            yvar, rmd, rms = ovar_reducedims(oyvar, sliceremovedim)
     28242                        else:
     28243                            yvar = oyvar[:]
     28244                            rmd = oyvar.dimensions
     28245                            rms = list(oxvar.shape)
     28246
     28247                        yn = np.min(yvar)
     28248                        yx = np.max(yvar)
     28249                        dyref = 1
     28250                        reflat1D = [(yn+yx)/2.]
     28251                        refblat1D = np.zeros((1,2), dtype=np.float)
     28252                        refblat1D[0,0] = nx
     28253                        refblat1D[0,1] = yx
     28254                        print '    using variable extremes instead: ', yn, yx
     28255                        dref.append('fake')
     28256
    2813128257                    if slicesinf.has_key(dimvars[xdim]):
    2813228258                        varslcv = slicesinf[dimvars[xdim]]
    … …  
    2813428260                        reflon1D = varslcv[3]
    2813528261                        refblon1D = varslcv[4]
     28262                        dref.append(varn)
    2813628263                    else:
    2813728264                        print infmsg
    … …  
    2815228279                        xx = np.max(xvar)
    2815328280                        dxref = 1
    28154                         reflon1D = (xn+xx)/2.
    28155                         refblon1D = np.zeros((1,4), dtype=np.float)
     28281                        reflon1D = [(xn+xx)/2.]
     28282                        refblon1D = np.zeros((1,2), dtype=np.float)
    2815628283                        refblon1D[0,0] = xn
    2815728284                        refblon1D[0,1] = xx
    28158                         refblon1D[0,2] = xx
    28159                         refblon1D[0,3] = xn
    2816028285                        print '    using variable extremes instead: ', xn, xx
    28161 
    28162                     if slicesinf.has_key(dimvars[ydim]):
    28163                         varslcv = slicesinf[dimvars[ydim]]
    28164                         dyref = varslcv[0]
    28165                         reflat1D = varslcv[3]
    28166                         refblat1D = varslcv[4]
    28167                     else:
    28168                         print infmsg
    28169                         print '  ' + fname + ": 2D slicing varibale '" + varn +      \
    28170                           "' with y-dimension '" + dimvars[ydim] + "' without " +    \
    28171                           "related slicing-variable !!"
    28172                         oyvar = onc.variables[dimvars[ydim]]
    28173 
    28174                         # Removing undesired dimensions from bounds slicing variable
    28175                         if sliceremovedim is not None:
    28176                             yvar, rmd, rms = ovar_reducedims(oyvar, sliceremovedim)
    28177                         else:
    28178                             yvar = oyvar[:]
    28179                             rmd = oyvar.dimensions
    28180                             rms = list(oxvar.shape)
    28181 
    28182                         yn = np.min(yvar)
    28183                         yx = np.max(yvar)
    28184                         dyref = 1
    28185                         reflat1D = (yn+yx)/2.
    28186                         refblat1D = np.zeros((1,4), dtype=np.float)
    28187                         refblat1D[0,0] = yx
    28188                         refblat1D[0,1] = yx
    28189                         refblat1D[0,2] = yn
    28190                         refblat1D[0,3] = yn
    28191                         print '    using variable extremes instead: ', yn, yx
    28192 
    28193                     # Wrong 3D reflon, reflat!!! lon, lat 1D on 2D var !!!
    28194  
     28286                        dref.append('fake')
     28287
    2819528288                    # Slicing following boundaries, dimensions are 1D thus expand to
    2819628289                    #    2D filling as NW, NE, SE, SW
    … …  
    2821428307                            yslice2D[3,j,i] = refblat1D[j,0]
    2821528308
     28309                    sref = [dyref, dxref]
     28310
    2821628311                    # get values
    28217                     if dimvars.has_key(xdim):
    28218                         ogetlon = onc.variables[dimvars[xdim]]
     28312                    if dimvars.has_key(xdim):
     28313                        if gen.searchInlist(new2Dvars.keys(), dimvars[xdim]):
     28314                            ogetlon = new2Dvars[xdim]
     28315                        else:
     28316                            ogetlon = onc.variables[dimvars[xdim]]
    2821928317                        # Removing undesired dimensions from bounds slicing variable
    2822028318                        if sliceremovedim is not None:
    … …  
    2823128329                        quit(-1)
    2823228330                    if dimvars.has_key(ydim):
    28233                         ogetlat = onc.variables[dimvars[ydim]]
     28331                        if gen.searchInlist(new2Dvars.keys(), dimvars[ydim]):
     28332                            ogetlat = new2Dvars[ydim]
     28333                        else:
     28334                            ogetlat = onc.variables[dimvars[ydim]]
    2823428335                        # Removing undesired dimensions from bounds slicing variable
    2823528336                        if sliceremovedim is not None:
    … …  
    2824728348
    2824828349                    getshape = dimsbnds[xdim]
     28350
    2824928351                    if len(getshape) == 3:
    2825028352                        dxget = getshape[2]
    … …  
    2825628358                        dyget = getshape[0]
    2825728359                        dvget = 4
     28360
    2825828361                    xdimvarbnds2D = dimvbnds[xdim]
    2825928362                    ydimvarbnds2D = dimvbnds[ydim]
    … …  
    2828228385
    2828328386                    # Values for file
    28284                     dref = [dimvars[ydim], dimvars[xdim]]
    28285                     sref = [dyref, dxref]
    2828628387                    dget = [ydim, xdim]
    2828728388                    sget = [dyget, dxget]
    2828828389
    28289                 # Issue with WRF's XLONG(time, sout_north, west_east) ....
    2829028390                else:
    2829128391                    print errormsg
    28292                     print '  '+fname+": rank ", dv.shape, " of slicing varibale '" + \
     28392                    print '  '+fname+": rank ",dv.shape, " of slicing varibale '" +\
    2829328393                      dn + "' not ready !!"
    2829428394                    print '    function only works with 1D and 2D slicing variables'
    … …  
    2830428404                getvarybndst = ydimvarbnds2D.transpose()
    2830528405
    28306                 print 'Lluis dxref, dyref, dvref:', dxref, dyref, dvref, 'shapes: reflon', reflon.shape, 'reflat:', reflat.shape, \
    28307                   'refvarxbnds:', xslice2D.shape, 'refvarybnds:', yslice2D.shape
    28308                 print 'Lluis dxget, dyget, dvget:', dxget, dyget, dvget, 'shapes: getlon', getlon.shape, 'getlat:', getlat.shape, \
    28309                   'getvarxbnds:', xdimvarbnds2D.shape, 'getvarybnds:', ydimvarbnds2D.shape
    28310 
    2831128406                Ngridsint, gridsint, areast, percenst =                              \
    2831228407                  fsci.module_scientific.grid_spacepercen(xcavals=reflont,           \
    … …  
    2833128426
    2833228427                Ngridsinmax = np.max(Ngridsin)
    28333                 print 'Lluis shapes: Ngridsin', Ngridsin.shape, 'gridsin:', gridsin.shape, \
    28334                   'areas:', areas.shape, 'percens:', percens.shape
     28428                if Ngridsinmax == 0:
     28429                    print errormsg
     28430                    print '  ' +fname + ": slices for '" + varn + "' without any " + \
     28431                      "value !!"
     28432                    quit(-1)
    2833528433        else:
    2833628434            dn = varn + ''
    … …  
    2836528463                    print "    bounded dimension '" + dimn + "' by '" + slicebndsdim[dimn] + "' !!"
    2836628464                    if iid == 0:
    28367                         ogetblat = onc.variables[slicebndsdim[dimn]]
     28465                        if new2Dvars.has_key(slicebndsdim[dimn]):
     28466                            ogetblat = new2Dvars[slicebndsdim[dimn]]
     28467                        else:
     28468                            ogetblat = onc.variables[slicebndsdim[dimn]]
    2836828469                        # Removing undesired dimensions from bounds slicing variable
    2836928470                        if sliceremovedim is not None:
    … …  
    2837428475                            rms = list(ogetblat.shape)
    2837528476                    elif iid == 1:
    28376                         ogetblon = onc.variables[slicebndsdim[dimn]]
     28477                        if new2Dvars.has_key(slicebndsdim[dimn]):
     28478                            ogetblon = new2Dvars[slicebndsdim[dimn]]
     28479                        else:
     28480                            ogetblon = onc.variables[slicebndsdim[dimn]]
    2837728481                        # Removing undesired dimensions from bounds slicing variable
    2837828482                        if sliceremovedim is not None:
    … …  
    2845228556            if not gen.searchInlist(onewnc.dimensions, 'fake'):
    2845328557                newdim = onewnc.createDimension('fake',1)
    28454                 newdim = onewnc.createDimension('slice_fake',1)
    28455                 newdim = onewnc.createDimension('fakebnds',4)
    28456                 newvar = onewnc.createVariable('fake','f',('fake'))
    28457                 newvar[:] = [0.]
    28458                 basicvardef(newvar, 'fake', 'fake variable for transforming ' +  \
    28459                   '1D variable without bounds', '-')
    28460                 newvar = onewnc.createVariable('slice_fake','f',('slice_fake'))
    28461                 newvar[:] = [0.]
    28462                 basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
    28463                   'for transforming variable without bounds', '-')
    28464                 newvar = onewnc.createVariable('slice_fake_bnds','f',            \
    28465                   ('fakebnds','fake'))
    28466                 fakebnds = np.zeros((4,1), dtype=np.float)
    28467                 fakebnds[0,0] = -1
    28468                 fakebnds[1,0] = 1
    28469                 fakebnds[2,0] = 1
    28470                 fakebnds[3,0] = -1
    28471                 newvar[:] = fakebnds
    28472                 basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
    28473                   'for transforming variable without bounds', '-')
     28558                if not gen.searchInlist(onewnc.dimensions, 'slice_fake'):
     28559                    newdim = onewnc.createDimension('slice_fake',1)
     28560                    newdim = onewnc.createDimension('fakebnds',4)
     28561                if not onewnc.variables.has_key('fake'):
     28562                    newvar = onewnc.createVariable('fake','f',('fake'))
     28563                    newvar[:] = [0.]
     28564                    basicvardef(newvar, 'fake', 'fake variable for transforming ' +  \
     28565                      '1D variable without bounds', '-')
     28566                    newvar = onewnc.createVariable('slice_fake','f',('slice_fake'))
     28567                    newvar[:] = [0.]
     28568                    basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
     28569                      'for transforming variable without bounds', '-')
     28570                    newvar = onewnc.createVariable('slice_fake_bnds','f',            \
     28571                      ('fakebnds','fake'))
     28572                    fakebnds = np.zeros((4,1), dtype=np.float)
     28573                    fakebnds[0,0] = -1
     28574                    fakebnds[1,0] = 1
     28575                    fakebnds[2,0] = 1
     28576                    fakebnds[3,0] = -1
     28577                    newvar[:] = fakebnds
     28578                    basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
     28579                      'for transforming variable without bounds', '-')
    2847428580                onewnc.sync()
    2847528581
    28476         if not gen.searchInlist(newslcvarns,dn):
     28582        #if not gen.searchInlist(newslcvarns,dn):
     28583        #    # Let's avoid memory issues...
     28584        #    newslcvarns.append(dn)
     28585        #    lvalues = [Ngridsin, gridsin, percens]
     28586        #    vardimbndslice[dn] = lvalues
     28587
     28588        if not gen.searchInlist(newslcvarns,varn):
    2847728589            # Let's avoid memory issues...
    28478             newslcvarns.append(dn)
     28590            newslcvarns.append(varn)
    2847928591            lvalues = [Ngridsin, gridsin, percens]
    28480             vardimbndslice[dn] = lvalues
     28592            vardimbndslice[varn] = lvalues
     28593            #dn = varn
    2848128594
    2848228595        if not onewnc.variables.has_key(dn+'gridin'):
    … …  
    2850528618
    2850628619            newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
    28507             newvar[:] = Ngridsin[:]
     28620            if Ngridsin.shape[1] == 1:
     28621                newvar[:,0] = Ngridsin[:,0]
     28622            else:
     28623                newvar[0,:] = Ngridsin[0,:]
    2850828624            basicvardef(newvar, dn+'Ngrid', "number of grids cells from " +          \
    2850928625              dn + " laying within slice", '-')
    … …  
    2853528651            onewnc.sync()
    2853628652
    28537             print 'Lluis dv shape:', dv.shape, 'Ngridsin:', Ngridsin.shape,          \
    28538               'gridsin:', gridsin.shape
    2853928653            if len(dv.shape) == 1:
    2854028654                j = 0
    … …  
    2864728761    #sliceinBt[...,1] = sliceinBt0[...,0]+1
    2864828762
     28763    iiB = 3
     28764    jjB = 0
     28765    iiA = 0
     28766    jjA = 3
    2864928767    NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
    2865028768      npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
    … …  
    2866228780    print "  Joined slice '" + newslcvarns[0] + "' & '" + newslcvarns[1] + "':",     \
    2866328781      NpointsAB.shape, 'maxin:', maxNpointsAB
     28782
     28783    if maxNpointsAB == 0:
     28784        print errormsg
     28785        print '  '  + fname + ": slice without any conincidence !!"
     28786        quit(-1)
    2866428787
    2866528788    if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
    … …  
    2873728860    #quit()
    2873828861
    28739     for islc in range(2,Nnewslcs):
    28740         dn = dn + '_' + newslcvarns[iscl+1]
    28741 
    28742         osliceNA = onewnc.variables[newslcvarns[0] + 'Ngrid']
    28743         osliceinA = onewnc.variables[newslcvarns[0] + 'gridin']
    28744         osliceaA = onewnc.variables[newslcvarns[0] + 'gridarea']
    28745         oslicepA = onewnc.variables[newslcvarns[0] + 'gridpercen']
    28746         osliceNB = onewnc.variables[newslcvarns[1] + 'Ngrid']
    28747         osliceinB = onewnc.variables[newslcvarns[1] + 'gridin']
    28748         osliceaB = onewnc.variables[newslcvarns[1] + 'gridarea']
    28749         oslicepB = onewnc.variables[newslcvarns[1] + 'gridpercen']
     28862    if Nnewslcs == 3:
     28863        islc == 2
     28864        prevdn = dn
     28865        dn = prevdn + '_' + newslcvarns[islc]
     28866
     28867        osliceNAB = onewnc.variables[prevdn + 'Ngrid']
     28868        osliceinAB = onewnc.variables[prevdn + 'gridin']
     28869        osliceaAB = onewnc.variables[prevdn + 'gridarea']
     28870        oslicepAB = onewnc.variables[prevdn + 'gridpercen']
     28871        osliceNC = onewnc.variables[newslcvarns[islc] + 'Ngrid']
     28872        osliceinC = onewnc.variables[newslcvarns[islc] + 'gridin']
     28873        osliceaC = onewnc.variables[newslcvarns[islc] + 'gridarea']
     28874        oslicepC = onewnc.variables[newslcvarns[islc] + 'gridpercen']
    2875028875   
    28751         sliceNA = osliceNA[:]
    28752         sliceinA = osliceinA[:]
    28753         sliceaA = osliceaA[:]
    28754         sliceNB = osliceNB[:]
    28755         sliceinB = osliceinB[:]
    28756         sliceaB = osliceaB[:]
    28757 
    28758         dxA = osliceNA.shape[1]
    28759         dyA = osliceNA.shape[0]
    28760         dxyA = osliceinA.shape[1]
    28761         dxB = osliceNB.shape[1]
    28762         dyB = osliceNB.shape[0]
    28763         dxyB = osliceinB.shape[1]
     28876        sliceNAB = osliceNAB[:]
     28877        sliceinAB = osliceinAB[:]
     28878        sliceaAB = osliceaAB[:]
     28879        sliceNC = osliceNC[:]
     28880        sliceinC = osliceinC[:]
     28881        sliceaC = osliceaC[:]
     28882
     28883        dxB = osliceNAB.shape[1]
     28884        dyB = osliceNAB.shape[0]
     28885        dxA = osliceNAB.shape[3]
     28886        dyA = osliceNAB.shape[2]
     28887        dxyAB = osliceinAB.shape[1]
     28888        dxC = osliceNC.shape[1]
     28889        dyC = osliceNC.shape[0]
     28890        dxyC = osliceinC.shape[1]
    2876428891   
    28765         Srgrid = list(osliceNB.dimensions) + list(osliceNA.dimensions)
     28892        Srgrid = list(osliceNB.dimensions) + list(osliceNAB.dimensions)
    2876628893        Sigrid = ['coord', dn+'gridin'] + Srgrid
    2876728894        Spgrid = [dn+'gridin'] + Srgrid
    2876828895   
    28769         sliceNAt = sliceNA.transpose()
    28770         sliceinAt = sliceinA.transpose()
    28771         sliceNBt = sliceNB.transpose()
    28772         sliceinBt = sliceinB.transpose()
     28896        NpointsABC = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     28897        pointsABC = np.zeros((2,dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     28898        inpointsAB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     28899        inpointsC = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     28900        print 'Lluis shapes sliceNAB:', sliceNAB.shape, 'sliceinAB:', sliceinAB.shape
     28901        # Following B-slices
     28902        for j in range(dyB):
     28903            for i in range(dxB):
     28904                sliceNAt = sliceNAB[j,i,:,:].transpose()
     28905                sliceinAt = sliceinAB[:,:,j,i,:,:].transpose()
     28906                sliceNBt = sliceNC.transpose()
     28907                sliceinBt = sliceinC.transpose()
    2877328908   
    28774         NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
    28775           npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
    28776           dxa=dxA, dya=dyA, dxya=dxyA, dxb=dxB, dyb=dyB, dxyb=dxyB)
    28777         NpointsAB = NpointsABt.transpose()
    28778         pointsAB = pointsABt.transpose()
    28779         inpointsA = inpA.transpose()
    28780         inpointsB = inpB.transpose()
    28781 
    28782         # Remembering that it is python (C-like...)
    28783         inpointsA = inpointsA-1
    28784         inpointsB = inpointsB-1
    28785 
    28786         maxNpointsAB = np.max(NpointsAB)
    28787         print "  Joined slice '" + newslcvarns[0] + "' & '" + newslcvarns[1] + "':",     \
    28788           NpointsAB.shape, 'maxin:', maxNpointsAB
     28909                print 'LLuis j i:', j,i, 'shapes sliceNAt:', sliceNAt.shape, 'sliceinAt:', sliceinAt.shape, \
     28910                  'sliceNBt:', sliceNBt.shape, 'sliceinBt:', sliceinBt.shape
     28911                print '  A dx dy dxy:', dxA, dyA, dxyAB, 'B dx dy dxy:', dxC, dyC, dxyC
     28912                NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
     28913                  npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
     28914                  dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxC, dyb=dyC, dxyb=dxyC)
     28915                NpointsABC[:,:,j,i,:,:] = NpointsABt.transpose()
     28916                pointsABC[:,:,:,:,j,i,:,:] = pointsABt.transpose()
     28917                inpointsAB[:,:,:,j,i,:,:] = inpA.transpose()
     28918                inpointsC[:,:,:,j,i,:,:] = inpB.transpose()
     28919
     28920                # Remembering that it is python (C-like...)
     28921                inpointsAB = inpointsAB-1
     28922                inpointsC = inpointsC-1
     28923
     28924        maxNpointsABC = np.max(NpointsAB)
     28925
     28926        print "  Joined slice '" + prevdn + "' & '" + newslcvarns[islc] + "':",      \
     28927          NpointsABC.shape, 'maxin:', maxNpointsABC
    2878928928
    2879028929        if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
    … …  
    2879228931   
    2879328932        newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
    28794         newvar[:] = NpointsAB[:]
    28795         basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +            \
     28933        newvar[:] = NpointsABC[:]
     28934        basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +        \
    2879628935          newslcvarns[0] + " laying within slice " + newslcvarns[1], '-')
    2879728936        newvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
    2879828937   
    28799         innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),                \
     28938        innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),            \
    2880028939          fill_value=gen.fillValueI)
    2880128940        basicvardef(innewvar, dn+'gridin', "coordinates of the grids cells from slice "+ \
    … …  
    2892529064    spt = sgpa.transpose()
    2892629065
     29066    print '  ' + infmsg
     29067    print '    final slice:', Ndims, 'sape:', sN.shape
     29068
    2892729069    # Checking
    2892829070    i = 3
    … …  
    2897529117                if gen.searchInlist(slicestatsdim, dnn):
    2897629118                    vardims.append(vdimn)
    28977                     varshape.append(dimtwrf)
     29119                    varshape.append(len(onc.dimensions[dnn]))
    2897829120                    if not gen.searchInlist(onewnc.dimensions,dnn):
    2897929121                        add_dims(onc, onewnc, [dnn])
    28980                         vardims.append(dn)
    28981                         varshape.append(len(onc.dimensions[dn]))
     29122                        vardims.append(dnn)
     29123                        varshape.append(len(onc.dimensions[dnn]))
    2898229124                    if not gen.searchInlist(onewnc.variables,vdimn): add_vars(onc,   \
    2898329125                      onewnc, [vdimn])
Note: See TracChangeset for help on using the changeset viewer.