Changeset 2270 in lmdz_wrf

Timestamp:

Dec 27, 2018, 7:30:13 PM (6 years ago)

Author:

lfita

Message:

Getting almost the final version of 'area_weighted'

Location:

trunk/tools

Files:

• module_scientific.f90 (modified) (5 diffs)
• nc_var_tools.py (modified) (13 diffs)

trunk/tools/module_scientific.f90

@@ -62 +62 @@
 !   series of r_k numbers
 ! SwapR_K*: Subroutine swaps the values of its two formal arguments.
+! unique_matrixRK2D: Subroutine to provide the unique values within a 2D RK matrix
 ! write_finaltrack_ascii: Subroutine to read the final trajectories into an ASCII file
 ! write_overlap_polys_ascii: Subroutine to write to an ASCII file the associated polygons at a given time-step
@@ -5011 +5012 @@
   END SUBROUTINE spacepercen
 
+  SUBROUTINE unique_matrixRK2D(dx, dy, dxy, matrix2D, Nunique, unique)
+  ! Subroutine to provide the unique values within a 2D RK matrix
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy, dxy
+    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: matrix2D
+    INTEGER, INTENT(out)                                 :: Nunique
+    REAL(r_k), DIMENSION(dxy), INTENT(out)               :: unique
+
+! Local
+    INTEGER                                              :: ix, iy, iu, minvalv
+    LOGICAL                                              :: single
+
+!!!!!!! Variables
+! dx, dy: dimensions of the matrix
+! dxy: dx*dy, maximum possible amount of different values
+! matrix2D: matgrix of values
+! Nunique: amount of unique values
+! unique: sorted from minimum to maximum vector with the unique values
+
+    fname = 'unique_matrixRK2D'
+
+    minvalv = MINVAL(matrix2D)
+
+    Nunique = 1
+    unique(1) = minvalv
+    DO ix= 1, dx
+      DO iy= 1, dy
+        single = .TRUE.
+        DO iu = 1, Nunique
+          IF (matrix2D(ix,iy) == unique(iu)) THEN
+            single = .FALSE.
+            EXIT
+          END IF
+        END DO
+        IF (single) THEN
+          Nunique = Nunique + 1
+          unique(Nunique) = matrix2D(ix,iy)
+        END IF
+      END DO
+    END DO
+
+    CALL sortR_K(unique, Nunique)
+
+  END SUBROUTINE unique_matrixRK2D
+
   SUBROUTINE spaceweightstats(varin, Ngridsin, gridsin, percentages, stats, varout, dxA, dyA, dxB,    &
     dyB, maxNgridsin, Lstats)
@@ -5018 +5066 @@
 
     INTEGER, INTENT(in)                                  :: dxA, dyA, dxB, dyB, maxNgridsin, Lstats
-    CHARACTER(len=60), INTENT(in)                        :: stats
+    CHARACTER(len=*), INTENT(in)                         :: stats
     INTEGER, DIMENSION(dxA,dyA), INTENT(in)              :: Ngridsin
     INTEGER, DIMENSION(dxA,dyA,maxNgridsin,2), INTENT(in):: gridsin
-    REAL(r_k), DIMENSION(dxB,dyB), INTENT(out)           :: varin
+    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: varin
     REAL(r_k), DIMENSION(dxA,dyA,maxNgridsin), INTENT(in):: percentages
     REAL(r_k), DIMENSION(dxA,dyA,Lstats), INTENT(out)    :: varout
 
 ! Local
-    INTEGER                                              :: ix, iy, iv, ii, jj
+    INTEGER                                              :: ix, iy, iv, ic, iu, ii, jj
+    INTEGER                                              :: Ncounts
+    CHARACTER(len=3)                                     :: val1S, val2S
+    CHARACTER(len=30)                                    :: val3S
+    REAL(r_k)                                            :: val1, val2
+    REAL(r_k), DIMENSION(Lstats)                         :: icounts
 
 !!!!!!! Variables
@@ -5041 +5094 @@
 !   'min': minimum value
 !   'max': maximum value
-!   'mean': mean value
-!   'mean2': quadratic mean value
-!   'stddev': standard deviation value
+!   'mean': space weighted mean value
+!   'mean2': space weighted quadratic mean value
+!   'stddev': space weighted standard deviation value
 !   'count': percentage of the space of matrix A covered by each different value of matrix B
 ! varout: output statistical variable
@@ -5062 +5115 @@
           END DO
         END DO
+      CASE('max')
+        varout = -fillVal64
+        DO ix=1, dxA
+          DO iy=1, dyA
+            DO iv=1, Ngridsin(ix,iy)
+              ii = gridsin(ix,iy,iv,1)
+              jj = gridsin(ix,iy,iv,2)
+              IF (varin(ii,jj) > varout(ix,iy,Lstats)) varout(ix,iy,1) = varin(ii,jj)
+            END DO
+          END DO
+        END DO
+      CASE('mean')
+        varout = zeroRK
+        DO ix=1, dxA
+          DO iy=1, dyA
+            DO iv=1, Ngridsin(ix,iy)
+              ii = gridsin(ix,iy,iv,1)
+              jj = gridsin(ix,iy,iv,2)
+              varout(ix,iy,1) = varout(ix,iy,1) + varin(ii,jj)*percentages(ix,iy,iv)
+            END DO
+          END DO
+        END DO
+      CASE('mean2')
+        varout = zeroRK
+        DO ix=1, dxA
+          DO iy=1, dyA
+            DO iv=1, Ngridsin(ix,iy)
+              ii = gridsin(ix,iy,iv,1)
+              jj = gridsin(ix,iy,iv,2)
+              varout(ix,iy,1) = varout(ix,iy,1) + percentages(ix,iy,iv)*(varin(ii,jj))**2
+            END DO
+            varout(ix,iy,1) = varout(ix,iy,1) / Ngridsin(ix,iy)
+          END DO
+        END DO
+      CASE('stddev')
+        varout = zeroRK
+        DO ix=1, dxA
+          DO iy=1, dyA
+            val1 = zeroRK
+            val2 = zeroRK
+            DO iv=1, Ngridsin(ix,iy)
+              ii = gridsin(ix,iy,iv,1)
+              jj = gridsin(ix,iy,iv,2)
+              val1 = val1 + varin(ii,jj)*percentages(ix,iy,iv)
+              val2 = val2 + percentages(ix,iy,iv)*(varin(ii,jj))**2
+            END DO
+            varout(ix,iy,1) = SQRT(val2 - val1**2)
+          END DO
+        END DO
+      CASE('count')
+        CALL unique_matrixRK2D(dxA, dyA, dxA*dyA, varin, Ncounts, icounts)
+        IF (Lstats /= Ncounts) THEN
+          WRITE(val1S,'(I3)')Lstats
+          WRITE(val2S,'(I3)')Ncounts
+          msg = "for 'count' different amount of passed categories: " // TRIM(val1S) //               &
+            ' and found ' // TRIM(val2S) 
+          CALL ErrMsg(msg, fname, -1)
+        END IF
+        varout = zeroRK
+        DO ix=1, dxA
+          DO iy=1, dyA
+            DO iv=1, Ngridsin(ix,iy)
+              ii = gridsin(ix,iy,iv,1)
+              jj = gridsin(ix,iy,iv,2)
+              ic = Index1DArrayR_K(icounts, Ncounts, varin(ii,jj))
+              IF (ic == -1) THEN
+                WRITE(val3S,'(f30.20)')varin(ii,jj)
+                msg = "value '" // val3S // "' for 'count' not found"
+                CALL ErrMSg(msg, fname, -1)
+              ELSE
+                varout(ix,iy,ic) = varout(ix,iy,ic) + percentages(ix,iy,iv)
+              END IF
+            END DO
+          END DO
+        END DO
       CASE DEFAULT
         msg = "statisitcs '" // TRIM(stats) // "' not ready !!" // CHAR(44) // " available ones: " // &

trunk/tools/nc_var_tools.py

@@ -27186 +27186 @@
         [stats]: ',' separated list of space-statistics within each `ref' grid to 
             provide for each grid cell from the 'get' data
-          'mean': weghted mean (always provided)
-          'min': minimum value and is associated weight
-          'max': maximum value and is associated weight
-          'mean2': weghted quadratic mean
+          'min': minimum value
+          'max': maximum value
+          'mean': weighted mean
+          'mean2': weighted quadratic mean
           'stddev': standard deviation
+          'count': percentage of the space of matrix `ref' covered by each different 
+            value of matrix `get'
       [ncfiles]= ',' separated list of [netcdfref]:[dimvalsref],[netcdfget]:[dimvalsget]
         [netcdfref]: netCDF file name to be used as reference or destination        
@@ -27221 +27223 @@
 
     statsavail = ['mean', 'min', 'max', 'mean2', 'stddev']
+    # Characteristics of the statistics
+    statnschars = {'mean': ['spcmean', 'spatial mean', 1],                           \
+      'min': ['spcmin', 'spatial minimum', 1],                                       \
+      'max': ['spcmax', 'spatial maximum', 1],                                       \
+      'mean2': ['spc2mean', 'spatial quadratic mean', 1],                            \
+      'stddev': ['spcstddev', 'spatial standard deviation', 1],                      \
+      'count': ['spccount', 'spatial coverage by value', 1]}
 
     expectargs = '[strict]:[stats]'
@@ -27226 +27235 @@
 
     strict = gen.Str_Bool(values.split(':')[0])
-    stats = values.split(':')[1].split(',')
-
-    for stn in stats:
+    statns = gen.str_list(values.split(':')[1], ',')
+
+    for stn in statns:
         if not gen.searchInlist(statsavail, stn):
             print errormsg
             print '  ' + fname + ": statistics '" + stn + "' not ready !!"
-            print '    available ones:', statsavail
+            print '    available ones:', statnschar.keys()
             quit(-1)
 
@@ -27437 +27446 @@
                 quit(-1)
 
-            slicedict = {getdimxn: slicedimx, getdimyn: slicedimy}
+            getslicedict = {getdimxn: slicedimx, getdimyn: slicedimy}
             ovarx = oncget.variables[getvarxn]
             ovarxbnds = oncget.variables[getvarxbndsn]
@@ -27443 +27452 @@
             ovarybnds = oncget.variables[getvarybndsn]
 
-            sgetvarx, dims = SliceVarDict(ovarx,slicedict)
-            sgetvarbndsx, dims = SliceVarDict(ovarxbnds,slicedict)
-            sgetvary, dims = SliceVarDict(ovary,slicedict)
-            sgetvarbndsy, dims = SliceVarDict(ovarybnds,slicedict)
+            sgetvarx, dims = SliceVarDict(ovarx,getslicedict)
+            sgetvarbndsx, dims = SliceVarDict(ovarxbnds,getslicedict)
+            sgetvary, dims = SliceVarDict(ovary,getslicedict)
+            sgetvarbndsy, dims = SliceVarDict(ovarybnds,getslicedict)
 
             getlon, getlat =  gen.lonlat2D(ovarx[tuple(sgetvarx)], ovary[tuple(sgetvary)])
@@ -27476 +27485 @@
     percens = percenst.transpose()
 
+    Ngridsinmax = np.max(Ngridsin)
+
     # Opening new file
     ofilen = fname + '.nc'
@@ -27484 +27495 @@
     newdim = onewnc.createDimension('lat',refdy)
     newdim = onewnc.createDimension('bnds',refmaxvert)
-    newdim = onewnc.createDimension('gridin',np.max(Ngridsin))
+    newdim = onewnc.createDimension('gridin',Ngridsinmax)
     newdim = onewnc.createDimension('coord',2)
 
@@ -27509 +27520 @@
     newvar = onewnc.createVariable('Ngrid','i',('lat','lon'))
     newvar[:] = Ngridsin[:]
-    basicvardef(newvar, 'Ngrid', "number of grids cells grom 'get' laying within" +  \
-      " 'ref'",'-')
+    basicvardef(newvar, 'Ngrid', "number of grids cells from .get. laying within" +  \
+      " .ref.",'-')
     newvar.setncattr('coordinates','lon lat')
 
     innewvar = onewnc.createVariable('gridin','i',('coord', 'gridin', 'lat', 'lon'))
-    basicvardef(innewvar,'gridin',"coordinates of the grids cells grom 'get' laying "+ \
-      "within 'ref'",'-')
+    basicvardef(innewvar,'gridin',"coordinates of the grids cells from .get. laying "+ \
+      "within .ref.",'-')
     innewvar.setncattr('coordinates','lon lat')
 
     pnewvar = onewnc.createVariable('gridpercen','f',('gridin', 'lat', 'lon'))
-    basicvardef(pnewvar,'gridpercen',"percentages of the grids cells grom 'get' " +   \
-      "laying within 'ref'", '1')
+    basicvardef(pnewvar,'gridpercen',"percentages of the grids cells from .get. " +   \
+      "laying within .ref.", '1')
     pnewvar.setncattr('coordinates','lon lat')
     for j in range(refdy):
@@ -27528 +27539 @@
     onewnc.sync()
 
+    # Getting the right sizes
+    gridsin = innewvar[:]
+    gridsint = gridsin.transpose()
+    percens = pnewvar[:]
+    percenst = percens.transpose()
+
     # Getting values
     if variable == 'all':
@@ -27536 +27553 @@
 
     # Removing variables
-    varns = ''
+    varns = []
     for vn in variables:
         if not oncget.variables.has_key(vn):
@@ -27549 +27566 @@
         ovar = oncget.variables[vn]
         vardims = ovar.dimensions
-        if not searchInlist(vardims, getdimxn) or searchInlist(vardims, getdimyn):
+        if not gen.searchInlist(vardims, getdimxn) or                                \
+          not gen.searchInlist(vardims, getdimyn):
             print '  ' + fname + ": get variable '" + vn + "' does not have space "+ \
-              " dimensions '" + getdimxn + "', '" +  getdimxn + "' not using it !!"
+              " dimensions '" + getdimxn + "', '" +  getdimyn + "' not using it !!"
             print "    variables' dimensions:", vardims
         else:
@@ -27557 +27575 @@
 
     for vn in varns:
-        ovar = onget.variables[vn]
-        
+        print '  ' + fname + ": using '" + vn + "' ..."
+        ovar = oncget.variables[vn]
+        svar, dimvar = SliceVarDict(ovar,getslicedict)
+        varvals = ovar[tuple(svar)]
+
+        varvalst = varvals.transpose()
+
+        for stn in statns:
+            print "    computing '" + stn + "' "
+            statchar = statnschars[stn]
+
+            if stn == 'count': 
+                svals = np.unique(varin)
+                statchar[2] = len(svals)
+
+            varout = fsci.module_scientific.spaceweightstats(varin=varvalst,         \
+              ngridsin=Ngridsint, gridsin=gridsint, percentages=percenst,            \
+              stats=stn, dxa=refdx, dya=refdy, dxb=getdx, dyb=getdy,                 \
+              maxngridsin=Ngridsinmax, lstats=statchar[2])
+
+            varout = varout.transpose()
+            if statchar[2] == 1:
+                newvar = onewnc.createVariable(vn + statchar[0], 'f', ('lat','lon'))
+                newvar[:] = varout[0,:,:]
+            else:
+                if stn == 'count':
+                    if not gen.searchInlist(onewnc.dimensions, 'count'): 
+                        newdim = onewnc.createDimension('count', statchar[2])
+                        newvar = onewnc.createVariable(vn + 'count', 'f', ('count'))
+                        newvar[:] = svals
+                        basicvardef(newvar, vn+'count', 'unique .get. values of ' +  \
+                          vn, ovar.units)
+
+                    newvar = onewnc.createVariable(vn + statchar[0], 'f', ('count',  \
+                      'lat','lon'))
+                    newvar[:] = varout[:]
+            basicvardef(newvar, vn + statchar[0], vn + ' ' + statchar[1] + " from "+ \
+              " .get. values", ovar.units)
+            newvar.setncattr('coordinates','lon lat')
+            onewnc.sync()
 
     add_global_PyNCplot(onewnc, 'nc_var_tools', fname, '1.0')
+    onewnc.setncattr('ref_file',netcdfrefn)
+    onewnc.setncattr('get_file',netcdfgetn)
+
+    onewnc.sync()
+    onewnc.close()
 
     print fname + ": Successfull witting of file '" + ofilen + "' !!"
@@ -27567 +27628 @@
 
 fvals= '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'
-area_weighted('yes:mean,min,max',fvals,'ct_values')
+area_weighted('yes:min,max,mean,stddev',fvals,'ct_values')