Changeset 2308 in lmdz_wrf for trunk

Timestamp:

Jan 31, 2019, 3:38:37 PM (6 years ago)

Author:

lfita

Message:

Working!!
Adding:

• `compute_slices_stats_areaweighted': Function to compute stats of variables of a file splitting variables by slices along sets of ranges for a series of variables weighting by the area covered by each grid (defined as polygon) within the slice

Location:

trunk/tools

Files:

• module_scientific.f90 (modified) (1 diff)
• nc_var.py (modified) (1 diff)
• nc_var_tools.py (modified) (15 diffs)

trunk/tools/module_scientific.f90

@@ -5719 +5719 @@
           END DO
 
-          IF (s1 == ss1 .AND. s2 == ss2 .AND. s3 == ss3) THEN
-            PRINT *, ' Lluis FOR Nin:', Nin
-            PRINT *, ' Lluis FOR values:', vin(:,ii3)
-            DO iv=1, Nin
-              PRINT *,' c:', iv, ':', gin(iv,:), ' p:', pin(iv)
-            END DO
-          END IF
-
           ! Computing along d3
           DO i3=1, di3

trunk/tools/nc_var.py

@@ -69 +69 @@
 ## e.g. # nc_var.py -o area_weighted -f 'reference_data.nc:lon;lon;lon_bnds;-1;lat;lat;lat_bnds;-1,get_data.nc:lon;lon;lon_bnds;-1;lat;lat;lat_bnds;-1' -S 'yes:min,max,mean,stddev,count' -v ct_values,xband_values,box_values,mosaic_values
 ## e.g. # nc_var.py -o area_weighted -f '/media/lluis/ExtDiskC_ext4/bkup/llamp/estudios/dominios/SA150k/geo_em.d01.nc:west_east;XLONG_M;WRFxbnds;-1;south_north;XLAT_M;WRFybnds;-1,/media/lluis/ExtDiskC_ext4/bkup/llamp/estudios/dominios/SA50k/geo_em.d01.nc:west_east;XLONG_M;WRFxbnds;-1;south_north;XLAT_M;WRFybnds;-1' -S 'no:mean' -v HGT_M
-## e.g. # nc_var.py -o compute_slices_stats_areaweighted -S 'XLONG,-74.,-36.4,4.;XLAT,-63.,19.,4.;HGT,500.,7000.,500.@Time|WRFtime:west_east|XLONG:south_north|XLAT@Time@south_north|lat_bnds,west_east|lon_bnds@Time' -f /home/lluis/PY/wrfout_d01_1995-01-01_00:00:00 -V T2
+## e.g. # nc_var.py -o compute_slices_stats_areaweighted -S 'XLONG,-74.,-36.4,4.;XLAT,-63.,19.,4.;HGT,500.,7000.,500.@Time|WRFtime:west_east|XLONG:south_north|XLAT@Time@west_east|lon_bnds,south_north|lat_bnds@XLONG|lat_bnds;lon_bnds,XLAT|lat_bnds;lon_bnds@Time' -f /home/lluis/estudios/Andes/space_weighted/wrfout_bnds.nc -v T2,U10,V10,Q2
 
 from optparse import OptionParser

trunk/tools/nc_var_tools.py

@@ -26476 +26476 @@
     for islc in range(Nslices1):
         vvslc = slcvar1[islc,]
-        slcv[islc,] = np.where(vvslc, 1, 0)
-    newvar = onewnc.createVariable(varn1+'sliced', 'i', tuple(slcdims))
+        slcv[islc,] = np.where(vvslc, 1, gen.fillValueI)
+    newvar = onewnc.createVariable(varn1+'sliced', 'i', tuple(slcdims),              \
+      fill_value=gen.fillValueI)
     newvar[:] = slcv[:]
     basicvardef(newvar, varn1+'sliced', 'sliced variable ' + varn1, v1u)
@@ -26499 +26500 @@
     for islc in range(Nslices2):
         vvslc = slcvar2[islc,]
-        slcv[islc,] = np.where(vvslc, 1, 0)
-    newvar = onewnc.createVariable(varn2+'sliced', 'i', tuple(slcdims))
+        slcv[islc,] = np.where(vvslc, 1, gen.fillValueI)
+    newvar = onewnc.createVariable(varn2+'sliced', 'i', tuple(slcdims),              \
+      fill_value=gen.fillValueI)
     newvar[:] = slcv[:]
     basicvardef(newvar, varn2+'sliced', 'sliced variable ' + varn2, v2u)
@@ -27769 +27771 @@
     if values == 'h':
         print fname + '_____________________________________________________________'
-        print compute_slices_stats.__doc__
+        print compute_slices_stats_areaweighted.__doc__
         quit()
 
@@ -27953 +27955 @@
         for islc in range(Nslices):
             vvslc =~slcvar[islc,]
-            slcv[islc,] = np.where(vvslc, 1, 0)
-        newvar = onewnc.createVariable(varn+'sliced', 'i', tuple(slcdims))
+            slcv[islc,] = np.where(vvslc, 1, gen.fillValueI)
+        newvar = onewnc.createVariable(varn+'sliced', 'i', tuple(slcdims),           \
+          fill_value=gen.fillValueI)
         newvar[:] = slcv[:]
         basicvardef(newvar, varn+'sliced', 'sliced variable ' + varn, vu)
@@ -28228 +28231 @@
                   nbvertexmax=dvget)
 
+                # Comming from Fortran (dx, dy) !!!
                 Ngridsin = Ngridsint.transpose()
-                gridsin = gridsint.transpose()
+                gridsin0 = gridsint.transpose()
+                gridsin = np.zeros(gridsin0.shape, dtype=int)
+                gridsin[0,...] = gridsin0[1,...]
+                gridsin[1,...] = gridsin0[0,...]
+                gridsin = gridsin - 1
+
                 areas = areast.transpose()
                 percens = percenst.transpose()
 
                 # Remember we are in C-like!
-                gridsin = gridsin - 1
 
                 Ngridsinmax = np.max(Ngridsin)
@@ -28335 +28343 @@
 
                 if Nt > 0:
+                    # Results come from Fortran... (dx, dy)
                     gridsin[1,0:Nt,0,islc] = ainds[0:Nt,0]
                     gridsin[0,0:Nt,0,islc] = ainds[0:Nt,1]
@@ -28486 +28495 @@
 
     onewnc.sync()
+
     # Getting full slice
     # NOTE: for simplicity let's assume that there is only one slice with bounds !!
@@ -28587 +28597 @@
     pnewvar.setncattr('coordinates',' '.join(Spgrid[::-1]))
 
+    jjA = 10
+    iiA = 2
+    jjB = 0
+    iiB = 4
     for jA in range(dyA):
         for iA in range(dxA):
@@ -28612 +28626 @@
 
     onewnc.sync()
+    #onewnc.close()
+    #quit()
 
     for islc in range(2,Nnewslcs):
@@ -28758 +28774 @@
                         # Let's check
                         psum = np.sum(sgpa[0:Nin,jA,iA,jB,iB])
-                        if psum != 1.:
-                            print warnmsg
-                            print '    ' + fname + 'at slide:', jA,iA,jB,iB,         \
-                              'percentages do NOT add one !!'
-                            print '    Ngrids:', Nin, 'slice area:',                 \
-                              aanewvar[jA,iA,jB,iB], 'total sum:', psum
-                            #for iv in range(Nin):
-                            #    print iv, sgpa[iv,jA,iA,jB,iB]
+                        if not np.isnan(psum):
+                            if not gen.almost_zero(psum, 1., 4):
+                                print warnmsg
+                                print '    ' + fname + 'at slide:', jA,iA,jB,iB,     \
+                                  'percentages do NOT add one down to 0.0001 !!'
+                                print '    Ngrids:', Nin, 'slice area:',             \
+                                  aanewvar[jA,iA,jB,iB], 'total sum:', psum
+                        else:
+                            print errormsg
+                            print '  ' + fname + ': grid Nan slice area for final '+ \
+                              ' slice:', jA, iA, jB, iB, 'N grids:', Nin,' !!'
+                            print '     slice total area:', aanewvar[jA,iA,jB,iB]
+                            print '     #grid grid_area grid_percen _______'
+                            for iv in range(Nin): 
+                                print iv, anewvar[iv,jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB]
+                            quit(-1)
     onewnc.sync()
                      
@@ -28785 +28809 @@
     #sint = sin.transpose()
     sint0 = sin.transpose()
-    ## We are in Fortran thus (dx,dy) and 1st:1 !!!!
+    ## We are going to use them in Fortran thus (dx,dy) and 1st:1 !!!!
     sint = np.ones((sint0.shape), dtype=np.float)*gen.fillValueF
     sint[...,0] = sint0[...,1]
     sint[...,1] = sint0[...,0]
+    sint = sint + 1
     spt = sgpa.transpose()
 
@@ -28798 +28823 @@
     it = 1
     slicesvals = []
-    print 'Lluis Ndims:', Ndims
     idd = 0
     for dnn in Ndims:
         scvv = slicesinf[dnn.split('_')[1]]
         sccvvv = scvv[3]
-        print 'sccvvv:', sccvvv
-        print 'Lluis dnn:', dnn, sccvvv[chijk[idd]]
         slicesvals.append(sccvvv[chijk[idd]])
         idd = idd + 1
@@ -28856 +28878 @@
         newvardims = vardims + dimsslice
         newvarshape = varshape + shapeslice
-        print 'Lluis newvardims:', newvardims, 'newvarshape:', newvarshape, 'varv:', varv.shape
-        print 'Lluis shapes sint:', sint.shape
 
 #        newvarslice = []
@@ -28884 +28904 @@
         varvt = None
         varvt = varv.transpose()
-        print 'vs:', varvt[7, 69, :]
         if len(vshape) == 3 and len(sN.shape) == 4:
             voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4(    \
@@ -28906 +28925 @@
             quit(-1)
 
-        print "Lluis Let's check !!", vout.shape
-        print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals
-        for iv in range(sN[i,j,k]):
-            print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:',             \
-              varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]]
-        print 'min max mean mean2 std median count _______'
-        print vout[:,it,i,j,k]
+        #print "Lluis Let's check !!", vout.shape
+        #print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals
+        #for iv in range(sN[i,j,k]):
+        #    print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:',             \
+        #      varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]]
+        #print 'min max mean mean2 std median count _______'
+        #print vout[:,it,i,j,k]
 
         ist = 0
         for statn in statns:
             newvar = onewvars[statn]
-            print 'Lluis stn:', statn, 'shape:', newvar.shape
             rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist])
             newvar[:] = rightvals