Changeset 2315 in lmdz_wrf for trunk/tools

Timestamp:

Feb 4, 2019, 11:40:15 PM (7 years ago)

Author:

lfita

Message:

Working on the triple slicing!

File:

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

trunk/tools/nc_var_tools.py

@@ -28890 +28890 @@
         dxyC = osliceinC.shape[1]
     
-        Srgrid = list(osliceNB.dimensions) + list(osliceNAB.dimensions)
+        Srgrid = list(osliceNC.dimensions) + list(osliceNAB.dimensions)
         Sigrid = ['coord', dn+'gridin'] + Srgrid
         Spgrid = [dn+'gridin'] + Srgrid
@@ -28897 +28897 @@
         pointsABC = np.zeros((2,dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
         inpointsAB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
+        NpointsA = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
+        NpointsB = np.zeros((dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
+        inpointsA = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
+        inpointsB = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
         inpointsC = np.zeros((dxyAB,dyC,dxC,dyB,dxB,dyA,dxA), dtype=int)
-        print 'Lluis shapes sliceNAB:', sliceNAB.shape, 'sliceinAB:', sliceinAB.shape
         # Following B-slices 
         for j in range(dyB):
@@ -28906 +28909 @@
                 sliceNBt = sliceNC.transpose()
                 sliceinBt = sliceinC.transpose()
-    
-                print 'LLuis j i:', j,i, 'shapes sliceNAt:', sliceNAt.shape, 'sliceinAt:', sliceinAt.shape, \
-                  'sliceNBt:', sliceNBt.shape, 'sliceinBt:', sliceinBt.shape
-                print '  A dx dy dxy:', dxA, dyA, dxyAB, 'B dx dy dxy:', dxC, dyC, dxyC
                 NpointsABt, pointsABt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
                   npointsa=sliceNAt, pointsa=sliceinAt, npointsb=sliceNBt, pointsb=sliceinBt,    \
@@ -28918 +28917 @@
                 inpointsC[:,:,:,j,i,:,:] = inpB.transpose()
 
+                # Getting respective A & B
+                for jc in range(dyC):
+                    for ic in range(dxC):
+                        print jc, ic, 'Lluis shapes NpointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, \
+                        'pointsABt[jc,ic,:,:]', NpointsABt[jc,ic,:,:].shape, 'sliceNAt', sliceNAt.shape, \
+                        'sliceinAt', sliceinAt.shape, 'A dx dy dxy:', dxA, dyA, dxyAB,                   \
+                        'B dx dy dxy:', dxA, dyA, dxA
+                        
+                        NptsAt, inptsAt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
+                          npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNAt,\
+                          pointsb=sliceinAt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxA, dyb=dyA, dxyb=dxyA)
+                        NpointsA[jc,ic,j,i,:,:] = NptsAt.transpose()
+                        inpointsA[:,jc,ic,j,i,:,:] = inptsAt.transpose()
+
+                        NptsBt, inptsBt, inpA, inpB = fsci.module_scientific.coincident_gridsin2d( \
+                          npointsa=NpointsABt[jc,ic,:,:], pointsa=pointsABt[jc,ic,:,:], npointsb=sliceNBt,\
+                          pointsb=sliceinBt, dxa=dxA, dya=dyA, dxya=dxyAB, dxb=dxB, dyb=dyB, dxyb=dxyB)
+                        NpointsB[jc,ic,j,i,:,:] = NptsAt.transpose()
+                        inpointsB[:,jc,ic,j,i,:,:] = inptsAt.transpose()
+
                 # Remembering that it is python (C-like...)
                 inpointsAB = inpointsAB-1
                 inpointsC = inpointsC-1
 
-        maxNpointsABC = np.max(NpointsAB)
+        maxNpointsABC = np.max(NpointsABC)
 
         print "  Joined slice '" + prevdn + "' & '" + newslcvarns[islc] + "':",      \
@@ -28928 +28947 @@
 
         if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
-            newdim = onewnc.createDimension(dn+'gridin', maxNpointsAB)
+            newdim = onewnc.createDimension(dn+'gridin', maxNpointsABC)
     
         newvar = onewnc.createVariable(dn+'Ngrid','i', tuple(Srgrid))
+        print 'Lluis shapes: newvar', newvar.shape, 'NpointsABC:', NpointsABC.shape, 'Srgrid:', Srgrid
         newvar[:] = NpointsABC[:]
         basicvardef(newvar, dn+'Ngrid', "number of grids cells from slice " +        \
@@ -28960 +28980 @@
         pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1]))
 
+        print 'Lluis shapes NpointsANC:', NpointsABC.shape, 'pointsABC:', pointsABC.shape, \
+          'inpoints AB:', inpointsAB.shape, 'sliceinAB:', sliceinAB.shape, 'oslicepAB:',   \
+          oslicepAB.shape, 'osliceaAB:', osliceaAB.shape, 'inpoints C:', inpointsC.shape,  \
+          'sliceinC:', sliceinC.shape, 'oslicepC:', oslicepC.shape, 'osliceaC:', osliceaC.shape
+        ixA = -1
+        iyA = -1
+        ixB = -1
+        iyB = -1
+        ixC = -1
+        iyC = -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]
-                        aanewvar[jB,iB,jA,iA]= gen.fillValueF
-                        anewvar[:,jB,iB,jA,iA]= gen.fillValueF
-                        pnewvar[:,jB,iB,jA,iA]= gen.fillValueF
-                        slicearea = 0.
-                        for iv in range(Nin):
-                            pA = oslicepA[inpointsA[iv,jB,iB,jA,iA],jA,iA]
-                            pB = oslicepB[inpointsB[iv,jB,iB,jA,iA],jB,iB]
-                            pnewvar[iv,jB,iB,jA,iA]= pA*pB
-                            ixA = sliceinA[1,inpointsA[iv,jB,iB,jA,iA],jA,iA]
-                            iyA = sliceinA[0,inpointsA[iv,jB,iB,jA,iA],jA,iA]
-                            aA = osliceaA[iyA,ixA]
-                            ixB = sliceinB[1,inpointsB[iv,jB,iB,jA,iA],jB,iB]
-                            iyB = sliceinB[0,inpointsB[iv,jB,iB,jA,iA],jB,iB]
-                            aB = osliceaB[iyB,ixB]
-                            anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB)
-                            slicearea = slicearea + (aA*pA)*(aB*pB)
-                        aanewvar[jB,iB,jA,iA]= slicearea
+                        for jC in range(dyC):
+                            for iC in range(dxC):
+                                Nin = NpointsABC[jC,iC,jB,iB,jA,iA]
+                                innewvar[:,0:Nin,jC,iC,jB,iB,jA,iA] = pointsABC[:,0:Nin,jC,iC,jB,iB,jA,iA]
+                                aanewvar[jC,iC,jB,iB,jA,iA]= gen.fillValueF
+                                anewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF
+                                pnewvar[:,jC,iC,jB,iB,jA,iA]= gen.fillValueF
+                                slicearea = 0.
+                                for iv in range(Nin):
+                                    pA = oslicepAB[inpointsAB[iv,jC,iC,jB,iB,jA,iA],jB,iB,jA,iA]
+                                    pB = oslicepC[inpointsC[iv,jC,iC,jB,iB,jA,iA],jC,iC]
+                                    pnewvar[iv,jC,iC,jB,iB,jA,iA]= pA*pB
+                                    ixAB = sliceinAB[1,inpointsAB[iv,jC,iC,jB,iB,jA,iA],jB,iB,jA,iA]
+                                    iyAB = sliceinAB[0,inpointsAB[iv,jC,iC,jB,iB,jA,iA],jB,iB,jA,iA]
+                                    ixA = inpointsA[1,NpointsA[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]
+                                    iyA = inpointsA[0,NpointsA[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]
+                                    ixB = inpointsB[1,NpointsB[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]
+                                    iyB = inpointsB[0,NpointsB[iv,jC,iC,jB,iB,jA,iA],jC,iC,jB,iB,jA,iA]
+                                    ixC = sliceinC[1,inpointsC[iv,jC,iC,jB,iB,jA,iA],jC,iC]
+                                    iyC = sliceinC[0,inpointsC[iv,jC,iC,jB,iB,jA,iA],jC,iC]
+                                    print jA,iA,jB,iB,jC,iC,':',iyA,ixA,iyB,ixB,iyC,ixC
+                                    aA = osliceaAB[iyB,ixB,iyA,ixA]
+                                    aB = osliceaC[iyC,ixC]
+                                    anewvar[iv,jC,iC,jB,iB,jA,iA]= (aA*pA)*(aB*pB)
+                                    slicearea = slicearea + (aA*pA)*(aB*pB)
+                                aanewvar[jC,iC,jB,iB,jA,iA]= slicearea
         onewnc.sync()