Changeset 2520 in lmdz_wrf for trunk


Ignore:
Timestamp:
May 9, 2019, 1:15:58 AM (6 years ago)
Author:
lfita
Message:

Adding 4-variable slicing into `compute_slices_stats_areaweighted'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var_tools.py

    r2518 r2520  
    2922929229                        refblat1D[0,1] = yx
    2923029230                        print '    using variable extremes instead: ', yn, yx
    29231                         dref.append('fake')
     29231                        Nfake = gen.ntimesHval_inlist(list(onewnc.dimensions), 'fake')
     29232                        faken = 'fake' + str(Nfake)
     29233                        dref.append(faken)
    2923229234
    2923329235                    if slicesinf.has_key(dimvars[xdim]):
    … …  
    2926029262                        refblon1D[0,1] = xx
    2926129263                        print '    using variable extremes instead: ', xn, xx
    29262                         dref.append('fake')
     29264                        Nfake = gen.ntimesHval_inlist(list(onewnc.dimensions), 'fake')
     29265                        faken = 'fake' + str(Nfake)
     29266                        dref.append(faken)
    2926329267
    2926429268                    # Slicing following boundaries, dimensions are 1D thus expand to
    … …  
    2952429528            dyget = 1
    2952529529
    29526             dref = ['fake', dn]
     29530            Nfake = gen.ntimesHval_inlist(list(onewnc.dimensions), 'fake')
     29531            faken = 'fake' + str(Nfake)
     29532            dref = [faken, dn]
    2952729533            sref = [1, Nslices]
    29528             dget = ['fake']
    29529             sget = [1]
    29530             if not gen.searchInlist(onewnc.dimensions, 'fake'):
    29531                 newdim = onewnc.createDimension('fake',1)
    29532                 if not gen.searchInlist(onewnc.dimensions, 'slice_fake'):
    29533                     newdim = onewnc.createDimension('slice_fake',1)
    29534                     newdim = onewnc.createDimension('fakebnds',4)
    29535                 if not onewnc.variables.has_key('fake'):
    29536                     newvar = onewnc.createVariable('fake','f',('fake'))
     29534            if not gen.searchInlist(onewnc.dimensions, faken):
     29535                newdim = onewnc.createDimension(faken,1)
     29536                if not gen.searchInlist(onewnc.dimensions, 'slice_'+faken):
     29537                    newdim = onewnc.createDimension('slice_'+faken,1)
     29538                    newdim = onewnc.createDimension(faken+'bnds',4)
     29539                if not onewnc.variables.has_key(faken):
     29540                    newvar = onewnc.createVariable(faken,'f',(faken))
    2953729541                    newvar[:] = [0.]
    29538                     basicvardef(newvar, 'fake', 'fake variable for transforming ' +  \
     29542                    basicvardef(newvar, faken, 'fake variable for transforming ' +   \
    2953929543                      '1D variable without bounds', '-')
    29540                     newvar = onewnc.createVariable('slice_fake','f',('slice_fake'))
     29544                    newvar = onewnc.createVariable('slice_'+faken, 'f', ('slice_' +  \
     29545                      faken))
    2954129546                    newvar[:] = [0.]
    29542                     basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
    29543                       'for transforming variable without bounds', '-')
    29544                     newvar = onewnc.createVariable('slice_fake_bnds','f',            \
    29545                       ('fakebnds','fake'))
     29547                    basicvardef(newvar, 'slice_'+faken, 'slicing along fake ' +      \
     29548                      'variable for transforming variable without bounds', '-')
     29549                    newvar = onewnc.createVariable('slice_'+faken+'_bnds','f',       \
     29550                      (faken+'bnds',faken))
    2954629551                    fakebnds = np.zeros((4,1), dtype=np.float)
    2954729552                    fakebnds[0,0] = -1
    … …  
    2955029555                    fakebnds[3,0] = -1
    2955129556                    newvar[:] = fakebnds
    29552                     basicvardef(newvar, 'slice_fake', 'slicing along fake variable '+\
     29557                    basicvardef(newvar,'slice_'+faken,'slicing along fake variable '+\
     29558                      'for transforming variable without bounds', '-')
     29559                onewnc.sync()
     29560            Nfake = gen.ntimesHval_inlist(list(onewnc.dimensions), 'fake')
     29561            faken = 'fake' + str(Nfake)
     29562            dget = [faken]
     29563            sget = [1]
     29564            if not gen.searchInlist(onewnc.dimensions, faken):
     29565                newdim = onewnc.createDimension(faken,1)
     29566                if not gen.searchInlist(onewnc.dimensions, 'slice_'+faken):
     29567                    newdim = onewnc.createDimension('slice_'+faken,1)
     29568                    newdim = onewnc.createDimension(faken+'bnds',4)
     29569                if not onewnc.variables.has_key(faken):
     29570                    newvar = onewnc.createVariable(faken,'f',(faken))
     29571                    newvar[:] = [0.]
     29572                    basicvardef(newvar, faken, 'fake variable for transforming ' +   \
     29573                      '1D variable without bounds', '-')
     29574                    newvar = onewnc.createVariable('slice_'+faken, 'f', ('slice_' +  \
     29575                      faken))
     29576                    newvar[:] = [0.]
     29577                    basicvardef(newvar, 'slice_'+faken, 'slicing along fake ' +      \
     29578                      'variable for transforming variable without bounds', '-')
     29579                    newvar = onewnc.createVariable('slice_'+faken+'_bnds','f',       \
     29580                      (faken+'bnds',faken))
     29581                    fakebnds = np.zeros((4,1), dtype=np.float)
     29582                    fakebnds[0,0] = -1
     29583                    fakebnds[1,0] = 1
     29584                    fakebnds[2,0] = 1
     29585                    fakebnds[3,0] = -1
     29586                    newvar[:] = fakebnds
     29587                    basicvardef(newvar,'slice_'+faken,'slicing along fake variable '+\
    2955329588                      'for transforming variable without bounds', '-')
    2955429589                onewnc.sync()
    … …  
    2983829873    #quit()
    2983929874
    29840     print 'Lluis Nnewslcs', Nnewslcs
    29841 
    2984229875    if Nnewslcs >= 3:
    29843         islc == 2
    29844         prevdn = dn
    29845         dn = prevdn + '-' + newslcvarns[islc]
    29846 
    29847         osliceNAB = onewnc.variables[prevdn + 'Ngrid']
    29848         osliceinAB = onewnc.variables[prevdn + 'gridin']
    29849         osliceaAB = onewnc.variables[prevdn + 'gridarea']
    29850         oslicepAB = onewnc.variables[prevdn + 'gridpercen']
    29851         osliceNC = onewnc.variables[newslcvarns[islc] + 'Ngrid']
    29852         osliceinC = onewnc.variables[newslcvarns[islc] + 'gridin']
    29853         osliceaC = onewnc.variables[newslcvarns[islc] + 'gridarea']
    29854         oslicepC = onewnc.variables[newslcvarns[islc] + 'gridpercen']
     29876        for islc in range(2,Nnewslcs):
     29877            prevdn = dn + ''
     29878            dn = prevdn + '-' + newslcvarns[islc]
     29879            print '      ' + fname + ' computing slice number:', islc+1, '...'
     29880            print '        joining:', dn
     29881
     29882            osliceNAB = onewnc.variables[prevdn + 'Ngrid']
     29883            osliceinAB = onewnc.variables[prevdn + 'gridin']
     29884            osliceaAB = onewnc.variables[prevdn + 'gridarea']
     29885            oslicepAB = onewnc.variables[prevdn + 'gridpercen']
     29886            osliceNC = onewnc.variables[newslcvarns[islc] + 'Ngrid']
     29887            osliceinC = onewnc.variables[newslcvarns[islc] + 'gridin']
     29888            osliceaC = onewnc.variables[newslcvarns[islc] + 'gridarea']
     29889            oslicepC = onewnc.variables[newslcvarns[islc] + 'gridpercen']
    2985529890   
    29856         sliceNAB = osliceNAB[:]
    29857         sliceinAB = osliceinAB[:]
    29858         sliceaAB = osliceaAB[:]
    29859         sliceNC = osliceNC[:]
    29860         sliceinC = osliceinC[:]
    29861         sliceaC = osliceaC[:]
    29862 
    29863         dxB = osliceNAB.shape[1]
    29864         dyB = osliceNAB.shape[0]
    29865         dxA = osliceNAB.shape[3]
    29866         dyA = osliceNAB.shape[2]
    29867         dxyAB = osliceinAB.shape[1]
    29868         dxC = osliceNC.shape[1]
    29869         dyC = osliceNC.shape[0]
    29870         dxyC = osliceinC.shape[1]
     29891            sliceNAB = osliceNAB[:]
     29892            sliceinAB = osliceinAB[:]
     29893            sliceaAB = osliceaAB[:]
     29894            sliceNC = osliceNC[:]
     29895            sliceinC = osliceinC[:]
     29896            sliceaC = osliceaC[:]
     29897
     29898            NABshape = osliceNAB.shape
     29899            inABshape = osliceinAB.shape
     29900            NCshape = osliceNC.shape
     29901            inCshape = osliceinC.shape
     29902
     29903            shapes = list(NCshape) + list(NABshape)
     29904            shapesjoin = [osliceinC.shape[0], osliceinAB.shape[1]]
     29905            shapesnew = [osliceinAB.shape[1]] + shapes
     29906            shapesptnew = [2, osliceinAB.shape[1]] + shapes
     29907
     29908            rankABC = len(shapes)
     29909
     29910#        dxB = osliceNAB.shape[1]
     29911#        dyB = osliceNAB.shape[0]
     29912#        dxA = osliceNAB.shape[3]
     29913#        dyA = osliceNAB.shape[2]
     29914#        dxyAB = osliceinAB.shape[1]
     29915#        dxC = osliceNC.shape[1]
     29916#        dyC = osliceNC.shape[0]
     29917#        dxyC = osliceinC.shape[1]
     29918
     29919
     29920            NABdims = osliceNAB.dimensions
     29921            inABdims = osliceinAB.dimensions
     29922            NCdims = osliceNC.dimensions
     29923            inCdims = osliceinC.dimensions
     29924
     29925            NABrank = len(NABdims)
     29926
     29927            Srgrid = list(NCdims) + list(NABdims)
     29928            Sigrid = ['coord', dn+'gridin'] + Srgrid
     29929            Spgrid = [dn+'gridin'] + Srgrid
    2987129930   
    29872         Srgrid = list(osliceNC.dimensions) + list(osliceNAB.dimensions)
    29873         Sigrid = ['coord', dn+'gridin'] + Srgrid
    29874         Spgrid = [dn+'gridin'] + Srgrid
    29875    
    29876         NpointsABC = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29877         pointsABC = np.zeros((2,dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29878         inpointsAB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29879         NpointsA = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29880         NpointsB = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29881         inpointsA = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29882         inpointsB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29883         inpointsC = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
    29884         # Following B-slices
    29885         for j in range(dyB):
    29886             for i in range(dxB):
    29887                 sliceNAt = sliceNAB[j,i,:,:].transpose()
    29888                 sliceinAt = sliceinAB[:,:,j,i,:,:].transpose()
     29931#        NpointsABC = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29932#        pointsABC = np.zeros((2,dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29933#        inpointsAB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29934#        NpointsA = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29935#        NpointsB = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29936#        inpointsA = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29937#        inpointsB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29938#        inpointsC = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
     29939            NpointsABC = np.zeros(tuple(shapes), dtype=int)
     29940            pointsABC = np.zeros(tuple(shapesptnew), dtype=int)
     29941            inpointsAB = np.zeros(tuple(shapesnew), dtype=int)
     29942            NpointsA = np.zeros(tuple(shapes), dtype=int)
     29943            NpointsB = np.zeros(tuple(shapes), dtype=int)
     29944            inpointsA = np.zeros(tuple(shapesnew), dtype=int)
     29945            inpointsB = np.zeros(tuple(shapesnew), dtype=int)
     29946            inpointsC = np.zeros(tuple(shapesnew), dtype=int)
     29947
     29948            print '        Lluis shapes dims_______'
     29949            print '        NpointsABC:', shapes, Srgrid
     29950            print '        pointsABC:', shapesptnew, Sigrid
     29951            print '        inpointsABC:', shapesnew, Spgrid
     29952
     29953            # Getting the slices to find the coincident points. It must be 2D
     29954            runNABdims = []
     29955            runinABdims = []
     29956            runrdims = []
     29957            runidims = []
     29958            runpdims = []
     29959            for iid in range(NABrank-2):
     29960                runNABdims.append(NABdims[iid])
     29961                runinABdims.append(NABdims[iid])
     29962                runrdims.append(NABdims[iid])
     29963                runidims.append(NABdims[iid])
     29964                runpdims.append(NABdims[iid])
     29965
     29966            print '         Lluis rundims _______', runpdims
     29967           
     29968            NABslices = gen.provide_slices(NABdims, NABshape, runNABdims)
     29969            inABslices = gen.provide_slices(inABdims, inABshape, runinABdims)
     29970            rslices = gen.provide_slices(Srgrid, shapes, runrdims)
     29971            islices = gen.provide_slices(Sigrid, shapesptnew, runrdims)
     29972            pslices = gen.provide_slices(Spgrid, shapesnew, runrdims)
     29973
     29974            print '    Lluis number of slices: NAB', len(NABslices), 'inAB',         \
     29975              len(inABslices), 'r', len(rslices), 'i', len(islices), 'p', len(pslices)
     29976
     29977#        # Following B-slices
     29978#        for j in range(dyB):
     29979#            for i in range(dxB):
     29980#                sliceNAt = sliceNAB[j,i,:,:].transpose()
     29981#                sliceinAt = sliceinAB[:,:,j,i,:,:].transpose()
     29982#                sliceNBt = sliceNC.transpose()
     29983#                sliceinBt = sliceinC.transpose()
     29984#                NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
     29985#                  npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
     29986#                  dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxC, dyb=dyC, dxyb=dxyC)
     29987#                NpointsABC[:,:,j,i,:,:] = NpointsABt.transpose()
     29988#                pointsABC[:,:,:,:,j,i,:,:] = pointsABt.transpose()
     29989#                inpointsAB[:,:,:,j,i,:,:] = inpA.transpose()
     29990#                inpointsC[:,:,:,j,i,:,:] = inpB.transpose()
     29991
     29992#                # Remembering that it is python (C-like...)
     29993#                inpointsAB[:,:,:,j,i,:,:] = inpointsAB[:,:,:,j,i,:,:]-1
     29994#                inpointsC[:,:,:,j,i,:,:] = inpointsC[:,:,:,j,i,:,:]-1
     29995
     29996            # Following slices
     29997            Nslcs = len(NABslices)
     29998            for isc in range(Nslcs):
     29999                NABslcv = NABslices[isc]
     30000                inABslcv = inABslices[isc]
     30001                rslcv = rslices[isc]
     30002                islcv = islices[isc]
     30003                pslcv = pslices[isc]
     30004
     30005                print isc, '  Lluis slice shape______'
     30006                print '       NAB', NABslcv, sliceNAB[tuple(NABslcv)].shape
     30007                print '       inAB', inABslcv, sliceinAB[tuple(inABslcv)].shape
     30008                print '       r', rslcv, NpointsABC[tuple(rslcv)].shape
     30009                print '       i', islcv, pointsABC[tuple(islcv)].shape
     30010                print '       p', pslcv, inpointsAB[tuple(pslcv)].shape
     30011
     30012                dxA = NABshape[NABrank-1]
     30013                dyA = NABshape[NABrank-2]
     30014                dxyAB = inABshape[1]
     30015                dxC = NCshape[1]
     30016                dyC = NCshape[0]
     30017                dxyC = inCshape[1]
     30018
     30019                sliceNAt = sliceNAB[tuple(NABslcv)].transpose()
     30020                sliceinAt = sliceinAB[tuple(inABslcv)].transpose()
    2988930021                sliceNBt = sliceNC.transpose()
    2989030022                sliceinBt = sliceinC.transpose()
    29891                 NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
    29892                   npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
     30023                NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d(     \
     30024                  npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,        \
    2989330025                  dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxC, dyb=dyC, dxyb=dxyC)
    29894                 NpointsABC[:,:,j,i,:,:] = NpointsABt.transpose()
    29895                 pointsABC[:,:,:,:,j,i,:,:] = pointsABt.transpose()
    29896                 inpointsAB[:,:,:,j,i,:,:] = inpA.transpose()
    29897                 inpointsC[:,:,:,j,i,:,:] = inpB.transpose()
    29898 
    29899                 # Getting respective A & B
    29900                 #for jc in range(dyC):
    29901                 #    for ic in range(dxC):
    29902                 #        print jc, ic, 'Lluis shapes NpointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, \
    29903                 #        'pointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, 'sliceNAt', sliceNAt.shape, \
    29904                 #        'sliceinAt', sliceinAt.shape, 'A dx dy dxy:', dxA, dyA, dxyAB,                   \
    29905                 #        'B dx dy dxy:', dxA, dyA, dxA
    29906                 #       
    29907                 #        NptsAt, inptsAt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
    29908                 #          npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNAt,\
    29909                 #          pointsb=sliceinAt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxA, dyb=dyA, dxyb=dxyA)
    29910                 #        NpointsA[jc,ic,j,i,:,:] = NptsAt.transpose()
    29911                 #        inpointsA[:,jc,ic,j,i,:,:] = inptsAt.transpose()
    29912                 #
    29913                 #        NptsBt, inptsBt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
    29914                 #          npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNBt,\
    29915                 #          pointsb=sliceinBt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxB, dyb=dyB, dxyb=dxyB)
    29916                 #        NpointsB[jc,ic,j,i,:,:] = NptsAt.transpose()
    29917                 #        inpointsB[:,jc,ic,j,i,:,:] = inptsAt.transpose()
     30026                NpointsABC[tuple(rslcv)] = NpointsABt.transpose()
     30027                pointsABC[tuple(islcv)] = pointsABt.transpose()
     30028                inpointsAB[tuple(pslcv)] = inpA.transpose()
     30029                inpointsC[tuple(pslcv)] = inpB.transpose()
    2991830030
    2991930031                # Remembering that it is python (C-like...)
    29920                 inpointsAB[:,:,:,j,i,:,:] = inpointsAB[:,:,:,j,i,:,:]-1
    29921                 inpointsC[:,:,:,j,i,:,:] = inpointsC[:,:,:,j,i,:,:]-1
    29922 
    29923         maxNpointsABC = np.max(NpointsABC)
    29924 
    29925         print "  Joined slice '" + prevdn + "' & '" + newslcvarns[islc] + "':",      \
    29926           NpointsABC.shape, 'maxin:', maxNpointsABC
    29927 
    29928         if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
    29929             newdim = onewnc.createDimension(dn+'gridin', maxNpointsABC)
     30032                inpointsAB[tuple(pslcv)] = inpointsAB[tuple(pslcv)]-1
     30033                inpointsC[tuple(pslcv)] = inpointsC[tuple(pslcv)]-1
    2993030034   
    29931         newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
    29932         newvar[:] = NpointsABC[:]
    29933         basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +        \
    29934           newslcvarns[0] + " laying within slice " + newslcvarns[1], '-')
    29935         newvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
     30035
     30036            maxNpointsABC = np.max(NpointsABC)
     30037
     30038            print "  Joined slice '" + prevdn + "' & '" + newslcvarns[islc] + "':",      \
     30039              NpointsABC.shape, 'maxin:', maxNpointsABC
     30040
     30041            if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
     30042                newdim = onewnc.createDimension(dn+'gridin', maxNpointsABC)
    2993630043   
    29937         innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),            \
    29938           fill_value=gen.fillValueI)
    29939         basicvardef(innewvar, dn+'gridin', "coordinates of the grids cells from slice "+ \
    29940           newslcvarns[0] + " laying within slice " + newslcvarns[1],'-')
    29941         innewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
    29942 
    29943         aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid),                    \
    29944           fill_value=gen.fillValueF)
    29945         basicvardef(aanewvar ,dn+'area', "area of the join slice " + newslcvarns[0] +    \
    29946           " & " + newslcvarns[1], vu)
    29947         aanewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
    29948 
    29949         anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(Spgrid),                 \
    29950           fill_value=gen.fillValueF)
    29951         basicvardef(anewvar ,dn+'gridarea', "area of the join slice " + newslcvarns[0] + \
    29952           " & " + newslcvarns[1], vu)
    29953         anewvar.setncattr('coordinates',' '.join(Spgrid[::-1]))
    29954 
    29955         pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid),               \
    29956           fill_value=gen.fillValueF)
    29957         basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +                    \
    29958           "grids cells from " + newslcvarns[0] + " laying within " + newslcvarns[1], '1')
    29959         pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1]))
    29960         ixA = -1
    29961         iyA = -1
    29962         ixB = -1
    29963         iyB = -1
    29964         ixC = -1
    29965         iyC = -1
    29966         for jA in range(dyA):
    29967             for iA in range(dxA):
    29968                 for jB in range(dyB):
    29969                     for iB in range(dxB):
    29970                         for jC in range(dyC):
    29971                             for iC in range(dxC):
    29972                                 Nin = NpointsABC[jC,iC,jB,iB,jA,iA]
    29973                                 innewvar[:,0:Nin,jC,iC,jB,iB,jA,iA] =                \
    29974                                   pointsABC[:,0:Nin,jC,iC,jB,iB,jA,iA]
    29975                                 aanewvar[jC,iC,jB,iB,jA,iA]= gen.fillValueF
    29976                                 anewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF
    29977                                 pnewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF
    29978                                 slicearea = 0.
    29979                                 for iv in range(Nin):
    29980                                     ivAB = inpointsAB[iv,jC,iC,jB,iB,jA,iA]
    29981                                     ivC = inpointsC[iv,jC,iC,jB,iB,jA,iA]
    29982                                     pA = oslicepAB[ivAB,jB,iB,jA,iA]
    29983                                     pB = oslicepC[ivC,jC,iC]
    29984                                     aA = osliceaAB[ivAB,jB,iB,jA,iA]
    29985                                     aB = osliceaC[ivC,jC,iC]
    29986                                     ixAB = sliceinAB[1,ivAB,jB,iB,jA,iA]
    29987                                     iyAB = sliceinAB[0,ivAB,jB,iB,jA,iA]
    29988                                     ixC = sliceinC[1,ivC,jC,iC]
    29989                                     iyC = sliceinC[0,ivC,jC,iC]
    29990                                     anewvar[iv,jC,iC,jB,iB,jA,iA]= aA*pA*pB
    29991                                     pnewvar[iv,jC,iC,jB,iB,jA,iA]= pA*pB
    29992                                     slicearea = slicearea + aA*pA*pB
    29993                                 if np.isnan(slicearea):
    29994                                     print errormsg
    29995                                     print '  ' + fname + ': Nan slice area for slice:',          \
    29996                                       jC, iC, jB, iB, jA, iA, 'N grids:', Nin,' !!'
    29997                                     print '     #grid AB #grid grid_area ' +         \
    29998                                       'grid_percen C # grid grid_area grid_percen____'
    29999                                     for iv in range(Nin):
    30000                                         print iv, inpointsAB[iv,jC,iC,jB,iB,jA,iA],  \
    30001                                           aA, pA, inpointsC[iv,jC,iC,jB,iB,jA,iA],   \
    30002                                           aB, pB
    30003                                     quit(-1)
    30004                                 aanewvar[jC,iC,jB,iB,jA,iA]= slicearea
    30005         onewnc.sync()
    30006 
     30044            newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
     30045            newvar[:] = NpointsABC[:]
     30046            basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +        \
     30047              newslcvarns[0] + " laying within slice " + newslcvarns[1], '-')
     30048            newvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
     30049   
     30050            innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),            \
     30051              fill_value=gen.fillValueI)
     30052            basicvardef(innewvar, dn+'gridin', "coordinates of the grids cells from slice "+ \
     30053              newslcvarns[0] + " laying within slice " + newslcvarns[1],'-')
     30054            innewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
     30055
     30056            aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid),                    \
     30057              fill_value=gen.fillValueF)
     30058            basicvardef(aanewvar ,dn+'area', "area of the join slice " + newslcvarns[0] +    \
     30059              " & " + newslcvarns[1], vu)
     30060            aanewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
     30061
     30062            anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(Spgrid),                 \
     30063              fill_value=gen.fillValueF)
     30064            basicvardef(anewvar ,dn+'gridarea', "area of the join slice " + newslcvarns[0] + \
     30065              " & " + newslcvarns[1], vu)
     30066            anewvar.setncattr('coordinates',' '.join(Spgrid[::-1]))
     30067
     30068            pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid),               \
     30069              fill_value=gen.fillValueF)
     30070            basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +                    \
     30071              "grids cells from " + newslcvarns[0] + " laying within " + newslcvarns[1], '1')
     30072            pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1]))
     30073
     30074#        ixA = -1
     30075#        iyA = -1
     30076#        ixB = -1
     30077#        iyB = -1
     30078#        ixC = -1
     30079#        iyC = -1
     30080#        for jA in range(dyA):
     30081#            for iA in range(dxA):
     30082#                for jB in range(dyB):
     30083#                    for iB in range(dxB):
     30084#                        for jC in range(dyC):
     30085#                            for iC in range(dxC):
     30086#                                Nin = NpointsABC[jC,iC,jB,iB,jA,iA]
     30087#                                innewvar[:,0:Nin,jC,iC,jB,iB,jA,iA] =                \
     30088#                                  pointsABC[:,0:Nin,jC,iC,jB,iB,jA,iA]
     30089#                                aanewvar[jC,iC,jB,iB,jA,iA]= gen.fillValueF
     30090#                                anewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF
     30091#                                pnewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF
     30092#                                slicearea = 0.
     30093#                                for iv in range(Nin):
     30094#                                    ivAB = inpointsAB[iv,jC,iC,jB,iB,jA,iA]
     30095#                                    ivC = inpointsC[iv,jC,iC,jB,iB,jA,iA]
     30096#                                    pA = oslicepAB[ivAB,jB,iB,jA,iA]
     30097#                                    pB = oslicepC[ivC,jC,iC]
     30098#                                    aA = osliceaAB[ivAB,jB,iB,jA,iA]
     30099#                                    aB = osliceaC[ivC,jC,iC]
     30100#                                    ixAB = sliceinAB[1,ivAB,jB,iB,jA,iA]
     30101#                                    iyAB = sliceinAB[0,ivAB,jB,iB,jA,iA]
     30102#                                    ixC = sliceinC[1,ivC,jC,iC]
     30103#                                    iyC = sliceinC[0,ivC,jC,iC]
     30104#                                    anewvar[iv,jC,iC,jB,iB,jA,iA]= aA*pA*pB
     30105#                                    pnewvar[iv,jC,iC,jB,iB,jA,iA]= pA*pB
     30106#                                    slicearea = slicearea + aA*pA*pB
     30107#                                if np.isnan(slicearea):
     30108#                                    print errormsg
     30109#                                    print '  ' + fname + ': Nan slice area for slice:',          \
     30110#                                      jC, iC, jB, iB, jA, iA, 'N grids:', Nin,' !!'
     30111#                                    print '     #grid AB #grid grid_area ' +         \
     30112#                                      'grid_percen C # grid grid_area grid_percen____'
     30113#                                    for iv in range(Nin):
     30114#                                        print iv, inpointsAB[iv,jC,iC,jB,iB,jA,iA],  \
     30115#                                          aA, pA, inpointsC[iv,jC,iC,jB,iB,jA,iA],   \
     30116#                                          aB, pB
     30117#                                    quit(-1)
     30118#                                aanewvar[jC,iC,jB,iB,jA,iA]= slicearea
     30119            # Getting actual slice values
     30120            rslcvals = gen.provide_slices(Srgrid, shapes, Srgrid)
     30121            iABCdims = np.ones((rankABC), dtype=int)*(-1)
     30122
     30123            Nslcs = len(rslcvals)
     30124            for islc in range(Nslcs):
     30125                slcv = rslcvals[islc]
     30126                Nin = NpointsABC[tuple(slcv)]
     30127                inslcv = [slice(0,2)] + [slice(0,Nin)] + slcv
     30128                pslcv = [slice(0,Nin)] + slcv
     30129                innewvar[tuple(inslcv)] = pointsABC[tuple(inslcv)]
     30130                aanewvar[tuple(slcv)]= gen.fillValueF
     30131                anewvar[tuple(pslcv)]= gen.fillValueF
     30132                pnewvar[tuple(pslcv)]= gen.fillValueF
     30133                slicearea = 0.
     30134                for iv in range(Nin):
     30135                    iislcv = [iv] + slcv
     30136                    ivAB = inpointsAB[tuple(iislcv)]
     30137                    ivC = inpointsC[tuple(iislcv)]
     30138                    pAslcv = [ivAB] + slcv[2:rankABC+1]
     30139                    pA = oslicepAB[tuple(pAslcv)]
     30140                    pCslcv = [ivC] + slcv[0:2]
     30141                    pB = oslicepC[tuple(pCslcv)]
     30142                    aA = osliceaAB[tuple(pAslcv)]
     30143                    aB = osliceaC[tuple(pCslcv)]
     30144                    slcvv = [1] + pAslcv
     30145                    ixAB = sliceinAB[tuple(slcvv)]
     30146                    slcvv = [0] + pAslcv
     30147                    iyAB = sliceinAB[tuple(slcvv)]
     30148                    slcvv = [1] + pCslcv
     30149                    ixC = sliceinC[tuple(slcvv)]
     30150                    slcvv = [0] + pCslcv
     30151                    iyC = sliceinC[tuple(slcvv)]
     30152                    anewvar[tuple(iislcv)]= aA*pA*pB
     30153                    pnewvar[tuple(iislcv)]= pA*pB
     30154                    slicearea = slicearea + aA*pA*pB
     30155                if np.isnan(slicearea):
     30156                    print errormsg
     30157                    print '  ' + fname + ': Nan slice area for slice:', slcv,    \
     30158                      'N grids:', Nin,' !!'
     30159                    print '     #grid AB #grid grid_area grid_percen C # grid ' +\
     30160                      'grid_area grid_percen____'
     30161                    for iv in range(Nin):
     30162                        iislcv = [iv] + slcv
     30163                        print iv, inpointsAB[tuple(iislcv)], aA, pA,             \
     30164                          inpointsC[tuple(iislcv)], aB, pB
     30165                        quit(-1)
     30166                aanewvar[tuple(slcv)]= slicearea
     30167   
     30168            onewnc.sync()
     30169   
    3000730170    slcvarns = slcvarns + newslcvarns
    3000830171
    … …  
    3009930262                                              sgpa[iv,jA,iA,jB,iB,jC,iC]
    3010030263                                        quit(-1)
     30264    elif len(sN.shape)/2 == 4:
     30265        (dyA, dxA, dyB, dxB, dyC, dxC, dyD, dxD) = sN.shape
     30266        for jA in range(dyA):
     30267            for iA in range(dxA):
     30268                for jB in range(dyB):
     30269                    for iB in range(dxB):
     30270                        for jC in range(dyC):
     30271                            for iC in range(dxC):
     30272                                for jD in range(dyD):
     30273                                    for iD in range(dxD):
     30274                                        Nin = sN[jA,iA,jB,iB,jC,iC,jD,iD]
     30275                                        aaval = aanewvar[jA,iA,jB,iB,jC,iC,jD,iD]
     30276                                        if Nin > 0 and aaval != 0.:
     30277                                            sgpa[0:Nin,jA,iA,jB,iB,jC,iC,jD,iD] =    \
     30278                                              anewvar[0:Nin,jA,iA,jB,iB,jC,iC,jD,iD]/\
     30279                                              aanewvar[jA,iA,jB,iB,jC,iC,jD,iD]
     30280                                            panewvar[0:Nin,jA,iA,jB,iB,jC,iC,jD,iD]= \
     30281                                              sgpa[0:Nin,jA,iA,jB,iB,jC,iC,jD,iD]
     30282                                            # Let's check
     30283                                            psum = np.sum(sgpa[0:Nin,jA,iA,jB,iB,jC, \
     30284                                              iC,jD,iD])
     30285                                            if not np.isnan(psum):
     30286                                                if not gen.almost_zero(psum, 1., 4):
     30287                                                    print warnmsg
     30288                                                    print '    '+fname + 'at slide:',\
     30289                                                      jA,iA,jB,iB,jC,iC,jD,iD,       \
     30290                                                      'percentages do NOT add one '+ \
     30291                                                      'down to 0.0001 !!'
     30292                                                    print '    Ngrids:', Nin,        \
     30293                                                      'slice area:', aaval,          \
     30294                                                      'total sum:', psum
     30295                                            else:
     30296                                                print errormsg
     30297                                                print '  ' + fname + ': grid NaN '+  \
     30298                                                  'slice area for final slice:', jA, \
     30299                                                  iA, jB, iB, jC, iC,jD,iD,          \
     30300                                                  'N grids:', Nin,' !!'
     30301                                                print '     slice total area:',      \
     30302                                                  aanewvar[jA,iA,jB,iB,jC,iC,jD,iD]
     30303                                                print '     #grid grid_area ' +      \
     30304                                                  'slice_area grid_percen ______'
     30305                                                for iv in range(Nin):
     30306                                                    print iv, anewvar[iv,jA,iA,jB,iB,\
     30307                                                      jC,iC,jD,iD], \
     30308                                                      aaval,                         \
     30309                                                      sgpa[iv,jA,iA,jB,iB,jC,iC,jD,iD]
     30310                                                quit(-1)
    3010130311    else:
    3010230312        print errormsg
Note: See TracChangeset for help on using the changeset viewer.