Changeset 2290 in lmdz_wrf for trunk

Timestamp:

Jan 25, 2019, 8:31:56 PM (7 years ago)

Author:

lfita

Message:

Getting there !

File:

• trunk/tools/nc_var_tools.py (modified) (12 diffs)

trunk/tools/nc_var_tools.py

@@ -27536 +27536 @@
     getvarybndst = getvarbndsy.transpose()
 
-    Ngridsint, gridsint, percenst = fsci.module_scientific.spacepercen(              \
+    Ngridsint, gridsint, areast, percenst = fsci.module_scientific.spacepercen(      \
       xcavals=reflont, ycavals=reflatt, xbavals=refvarxbndst, ybavals=refvarybndst,  \
       xcbvals=getlont, ycbvals=getlatt, xbbvals=getvarxbndst, ybbvals=getvarybndst,  \
-      strict=strict, dxa=refdx, dya=refdy, navertexmax=refmaxvert,                   \
-      dxb=getdx, dyb=getdy, dxyb=getdx*getdy, nbvertexmax=getmaxvert)
+      strict=strict, dxa=refdx, dya=refdy, navertexmax=refmaxvert, dxb=getdx,        \
+      dyb=getdy, dxyb=getdx*getdy, nbvertexmax=getmaxvert)
 
     Ngridsin = Ngridsint.transpose()
     gridsin = gridsint.transpose()
+    percens = areast.transpose()
     percens = percenst.transpose()
 
@@ -27893 +27894 @@
             else: slcvarnS = slcvarnS + ', ' + varn
 
+        for dimn in rmdims:
+            if not gen.searchInlist(onewnc.dimensions, dimn):
+                newdim = onewnc.createDimension(dimn, len(onc.dimensions[dimn]))
+
+        newvar = onewnc.createVariable(varn, 'f', (rmdims))
+        newvar[:] = varv[:]
+        basicvardef(newvar, varn, 'variable ' + varn + ' used to slice', vu)
+        onewnc.sync()
+
         # dimensions
         newdim = onewnc.createDimension('slice_'+varn, Nslices)
@@ -27903 +27913 @@
         newvar[:] = slcvalsc[:]
         basicvardef(newvar, 'slice_'+varn, 'slices for variable ' + varn, vu)
+        newattr = set_attribute(newvar, 'bounds', 'slice_'+varn+'_bnds')
         onewnc.sync()
 
@@ -27916 +27927 @@
         slcv = np.zeros(tuple(slcshape), dtype=int)
         for islc in range(Nslices):
-            vvslc = ~slcvar[islc,]
+            vvslc =~slcvar[islc,]
             slcv[islc,] = np.where(vvslc, 1, 0)
         newvar = onewnc.createVariable(varn+'sliced', 'i', tuple(slcdims))
@@ -27927 +27938 @@
     vardimbndslice = {}
 
-    print 'Lluis slcvarns:', slcvarns
     newslcvarns = []
     for varn in slcvarns:
@@ -27945 +27955 @@
             rmshape = list(ovar.shape)
 
+        varfinaldims = rmdims + []
         dimbnds = []
         for dn in rmdims:
@@ -28183 +28194 @@
                 getvarybndst = ydimvarbnds2D.transpose()
 
-                Ngridsint, gridsint, percenst =                                      \
-                  fsci.module_scientific.spacepercen(xcavals=reflont,                \
+                Ngridsint, gridsint, areast, percenst =                              \
+                  fsci.module_scientific.grid_spacepercen(xcavals=reflont,           \
                   ycavals=reflatt, xbavals=refvarxbndst, ybavals=refvarybndst,       \
                   xcbvals=getlont, ycbvals=getlatt, xbbvals=getvarxbndst,            \
@@ -28193 +28204 @@
                 Ngridsin = Ngridsint.transpose()
                 gridsin = gridsint.transpose()
+                areas = areast.transpose()
                 percens = percenst.transpose()
 
                 Ngridsinmax = np.max(Ngridsin)
         else:
+            dn = varn + ''
             # Here is only done by 1D variables !!!
             print infmsg
@@ -28214 +28227 @@
 
             Ngridsinmax = np.max(Ngridsin)
+
             gridsin = np.zeros((2,Ngridsinmax,1,Nslices), dtype=int)
+            areas = np.zeros((Ngridsmax,1,Nslices), dtype=np.float)
             percens = np.zeros((Ngridsinmax,1,Nslices), dtype=np.float)
             for islc in range(Nslices):
@@ -28249 +28264 @@
                     percens[0:Ngridsin[0,islc],0,islc] = 1.
 
-        # Let's avoid memory issues...
-        newslcvarns.append(dn)
-        lvalues = [Ngridsin, gridsin, percens]
-        vardimbndslice[dn] = lvalues
+                onewnc.sync()
+
+        if not gen.searchInlist(newslcvarns,dn):
+            # Let's avoid memory issues...
+            newslcvarns.append(dn)
+            lvalues = [Ngridsin, gridsin, percens]
+            vardimbndslice[dn] = lvalues
+
         if not onewnc.variables.has_key(dn+'gridin'):
 
@@ -28283 +28302 @@
             newvar.setncattr('coordinates',' '.join(Scgrid[::-1]))
     
-            innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid))
+            innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),        \
+              fill_value=gen.fillValueI)
             basicvardef(innewvar, dn+'gridin', "coordinates of the grids " +         \
               "cells from .get. laying within .ref.",'-')
             innewvar.setncattr('coordinates',' '.join(Scgrid[::-1]))
-        
-            pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid))
+
+            anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(varfinaldims))
+            basicvardef(anewvar ,dn+'gridarea', "area of the grids cells from " +    \
+             ".get.", vu)
+            anewvar.setncattr('coordinates',' '.join(varfinaldims[::-1]))
+            anewvar[:] = areas[:]
+
+            pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid),       \
+              fill_value=gen.fillValueF)
             basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +            \
              "grids cells from .get. laying within .ref.", '1')
@@ -28307 +28334 @@
     
             onewnc.sync()
+
+    # Getting full slice
+    # NOTE: for simplicity let's assume that there is only one slice with bounds !!
+    # Otherwise we need to re-calculate the polygons taking the outcome from each 
+    #   slice
+    print 'joining slices of:', newslcvarns
+
+    Nnewslcs = len(newslcvarns)
+
+    dn = newslcvarns[0] + '_' + newslcvarns[1]
+    
+    osliceNA = onewnc.variables[newslcvarns[0] + 'Ngrid']
+    osliceinA = onewnc.variables[newslcvarns[0] + 'gridin']
+    osliceaA = onewnc.variables[newslcvarns[0] + 'gridarea']
+    oslicepA = onewnc.variables[newslcvarns[0] + 'gridpercen']
+    osliceNB = onewnc.variables[newslcvarns[1] + 'Ngrid']
+    osliceinB = onewnc.variables[newslcvarns[1] + 'gridin']
+    osliceaB = onewnc.variables[newslcvarns[1] + 'gridarea']
+    oslicepB = onewnc.variables[newslcvarns[1] + 'gridpercen']
+
+    sliceNA = osliceNA[:]
+    sliceinA = osliceinA[:]
+    sliceNB = osliceNB[:]
+    sliceinB = osliceinB[:]
+
+    dxA = osliceNA.shape[1]
+    dyA = osliceNA.shape[0]
+    dxyA = osliceinA.shape[1]
+    dxB = osliceNB.shape[1]
+    dyB = osliceNB.shape[0]
+    dxyB = osliceinB.shape[1]
+
+    Srgrid = list(osliceNA.dimensions) + list(osliceNB.dimensions)
+    Sigrid = ['coord', dn+'gridin'] + Srgrid
+    Spgrid = [dn+'gridin'] + Srgrid
+
+    sliceNAt = sliceNA.transpose()
+    sliceinAt = sliceinA.transpose()
+    sliceNBt = sliceNB.transpose()
+    sliceinBt = sliceinB.transpose()
+
+    NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
+      npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
+      dxa=dxA, dya=dyA, dxya=dxyA, dxb=dxB, dyb=dyB, dxyb=dxyB)
+    NpointsAB = NpointsABt.transpose()
+    pointsAB = pointsABt.transpose()
+    inpointsA = inpA.transpose()
+    inpointsB = inpB.transpose()
+
+    print "  Joined slice '" + newslcvarns[0] + "' & '" + newslcvarns[1] + "':",     \
+      NpointsAB.shape
+    quit(-1)
+
+    maxNpointsAB = np.max(NpointsAB)
+
+    if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
+        newdim = onewnc.createDimension(dn+'gridin', maxNpointsAB)
+
+    newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
+    newvar[:] = NpointsAB[:]
+    basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +            \
+      newslcvarns[0] + " laying within slice " + newslcvarns[1], '-')
+    newvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+    
+    innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),                \
+      fill_value=gen.fillValueI)
+    basicvardef(innewvar, dn+'gridin', "coordinates of the grids cells from slice "+ \
+      newslcvarns[0] + " laying within slice " + newslcvarns[1],'-')
+    innewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+
+    anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(Srgrid),                 \
+      fill_value=gen.fillValueF)
+    basicvardef(anewvar ,dn+'gridarea', "area of the join slice " + newslcvarns[0] + \
+      " & " + newslcvarns[1], vu)
+    anewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+
+    pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid),               \
+      fill_value=gen.fillValueF)
+    basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +                    \
+      "grids cells from .get. laying within .ref.", '1')
+    pnewvar.setncattr('coordinates',' '.join(Scgrid[::-1]))
+
+    for jA in range(dyA):
+        for iA in range(dxA):
+            for jB in range(dyB):
+                for iB in range(dxB):
+                    Nin = NpointsAB[jB,iB,jA,iA]
+                    innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA]
+                    print Nin, 'Lluis shapes sliceaA:', osliceaA.shape, 'sliceaB:', osliceaB.shape
+                    for iv in range(Nin):
+                        pA = oslicepA[inpA[jA,iA,iv],jA,iA]
+                        pB = oslicepB[inpB[jB,iB,iv],jB,iB]
+                        pnewvar[iv,jB,iB,jA,iA]= pA*pB
+                        ixA = sliceinA[0,inpA[jA,iA,iv],jA,iA]
+                        iyA = sliceinA[1,inpA[jA,iA,iv],jA,iA]
+                        print '  A:', ixA, iyA
+                        aA = osliceaA[ixA,iyA]
+                        ixB = sliceinB[0,inpB[jB,iB,iv],jB,iB]
+                        iyB = sliceinB[1,inpB[jB,iB,iv],jB,iB]
+                        print '  B:', ixB, iyB
+                        aB = osliceaB[ixB,iyB]
+                        anewvar[jB,iB,jA,iA]= aA*aB
+
+    onewnc.sync()
+
+    for islc in range(2,Nnewslcs):
+        dn = dn + '_' + newslcvarns[iscl+1]
+        osliceNA = onewnc.variables[newslcvarns[0] + 'Ngrid']
+        osliceinA = onewnc.variables[newslcvarns[0] + 'gridin']
+        osliceaA = onewnc.variables[newslcvarns[0] + 'gridarea']
+        oslicepA = onewnc.variables[newslcvarns[0] + 'gridpercen']
+        osliceNB = onewnc.variables[newslcvarns[1] + 'Ngrid']
+        osliceinB = onewnc.variables[newslcvarns[1] + 'gridin']
+        osliceaB = onewnc.variables[newslcvarns[1] + 'gridarea']
+        oslicepB = onewnc.variables[newslcvarns[1] + 'gridpercen']
+
+        sliceNA = osliceNA[:]
+        sliceinA = osliceinA[:]
+        sliceNB = osliceNB[:]
+        sliceinB = osliceinB[:]
+
+        dxA = osliceNA.shape[1]
+        dyA = osliceNA.shape[0]
+        dxyA = osliceinA.shape[1]
+        dxB = osliceNB.shape[1]
+        dyB = osliceNB.shape[0]
+        dxyB = osliceinB.shape[1]
+
+        Srgrid = list(osliceNA.dimensions) + list(osliceNB.dimensions)
+        Sigrid = ['coord', dn+'gridin'] + Srgrid
+        Spgrid = [dn+'gridin'] + Srgrid
+
+        sliceNAt = sliceNA.transpose()
+        sliceinAt = sliceinA.transpose()
+        sliceNBt = sliceNB.transpose()
+        sliceinBt = sliceinB.transpose()
+
+        NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
+          npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
+          dxa=dxA, dya=dyA, dxya=dxyA, dxb=dxB, dyb=dyB, dxyb=dxyB)
+        NpointsAB = NpointsABt.transpose()
+        pointsAB = pointsABt.transpose()
+        inpointsA = inpA.transpose()
+        inpointsB = inpB.transpose()
+
+        maxNpointsAB = np.max(NpointsAB)
+
+        if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
+            newdim = onewnc.createDimension(dn+'gridin', maxNpointsAB)
+    
+        newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
+        newvar[:] = NpointsAB[:]
+        basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +            \
+          newslcvarns[0] + " laying within slice " + newslcvarns[1], '-')
+        newvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+    
+        innewvar = onewnc.createVariable(dn+'gridin', 'i', tuple(Sigrid),                \
+          fill_value=gen.fillValueI)
+        basicvardef(innewvar, dn+'gridin', "coordinates of the grids cells from slice "+ \
+          newslcvarns[0] + " laying within slice " + newslcvarns[1],'-')
+        innewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+
+        anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(Srgrid),                 \
+          fill_value=gen.fillValueF)
+        basicvardef(anewvar ,dn+'gridarea', "area of the join slice " + newslcvarns[0] + \
+          " & " + newslcvarns[1], vu)
+        anewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+
+        pnewvar = onewnc.createVariable(dn+'gridpercen','f',tuple(Spgrid),               \
+          fill_value=gen.fillValueF)
+        basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +                    \
+          "grids cells from .get. laying within .ref.", '1')
+        pnewvar.setncattr('coordinates',' '.join(Scgrid[::-1]))
+
+        for jA in range(dyA):
+            for iA in range(dxA):
+                for jB in range(dyB):
+                    for iB in range(dxB):
+                        Nin = NpointsAB[jB,iB,jA,iA]
+                        innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA]
+                        for iv in range(Nin):
+                            pA = oslicepA[inpA[iv],jA,iA]
+                            pB = oslicepB[inpB[iv],jB,iB]
+                            pnewvar[iv,jB,iB,jA,iA]= pA*pB
+                            aA = osliceaA[inpA[iv],jA,iA]
+                            aB = osliceaB[inpB[iv],jB,iB]
+                            anewvar[iv,jB,iB,jA,iA]= aA*aB
+
+        onewnc.sync()
 
     slcvarns = slcvarns + newslcvarns