Changeset 2297 in lmdz_wrf

Timestamp:

Jan 29, 2019, 5:29:33 PM (6 years ago)

Author:

lfita

Message:

Implementing the calculation of the weighted statistics

Location:

trunk/tools

Files:

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

trunk/tools/module_scientific.f90

@@ -5634 +5634 @@
   END SUBROUTINE spaceweightstats
 
+  SUBROUTINE multi_spacewightstats_in3DRK3_slc3v(varin, rundim, Ngridsin, gridsin, percentages,       &
+    varout, di1, di2, di3, ds1, ds2, ds3, ds4, ds5, ds6, maxNgridsin, Lstats)
+  ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as 
+  !   running one into a matrix of slices using spatial weights
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: di1, di2, di3, ds1, ds2, ds3, ds4, ds5, ds6
+    INTEGER, INTENT(in)                                  :: rundim, maxNgridsin
+    CHARACTER(len=*), INTENT(in)                         :: stats
+    INTEGER, DIMENSION(ds1,ds2,ds3,ds4,ds5,ds6),                                                      &
+      INTENT(in)                                         :: Ngridsin
+    INTEGER, INTENT(in)                                                                               &
+      DIMENSION(ds1,ds2,ds3,ds4,ds5,ds6,maxNgridsin,2)   :: gridsin
+    REAL(r_k), DIMENSION(di1,di2,di3), INTENT(in)        :: varin
+    REAL(r_k), INTENT(in)                                                                             &
+      DIMENSION(ds1,ds2,ds3,ds4,ds5,ds6,maxNgridsin)     :: percentages
+    REAL(r_k), DIMENSION(ds1,ds2,ds3,ds4,ds5,ds6,di3,7),                                              &
+      INTENT(out)                                        :: varout
+
+! Local
+    INTEGER                                              :: i1, i2, i3, s1, s2, s3, s4, s5, s6
+    INTEGER                                              :: Ncounts, Nin
+    CHARACTER(len=3)                                     :: val1S, val2S
+    CHARACTER(len=30)                                    :: val3S
+    REAL(r_k)                                            :: val1, val2
+    REAL(r_k), DIMENSION(Lstats)                         :: icounts
+    INTEGER, DIMENSION(:), ALLOCATABLE                   :: pin
+    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: gin
+
+!!!!!!! Variables
+! di1, di2, di3: length of dimensions of the 3D matrix of values
+! ds[1-6]: length of dimensions of matrix with the slices
+! maxNgridsin: maximum number of grid points from the 3D matrix in any slice
+! varin: 3D RK variable to be used
+! Ngridsin: number of grids from 3D RK matrix for each slice
+! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices
+! percentages: weights as percentages of space of 3D RK matrix for each slice
+! stats: name of the spatial statistics to compute inside each slice using values from 3D RK matrix
+!   Avaialbe ones:
+!   'min': minimum value
+!   'max': maximum value
+!   'mean': space weighted mean value
+!   'mean2': space weighted quadratic mean value
+!   'stddev': space weighted standard deviation value
+!   'median': median value
+!   'count': percentage of the space of matrix A covered by each different value of matrix B
+! varout: output statistical variable
+
+    fname = 'multi_spacewightstats_in3DRK3_slc3v'
+
+    ! Let's be efficient?
+    varout = fillVal64
+    DO s1 =1, ds1
+      DO s2 =1, ds2
+        DO s3 =1, ds3
+          DO s4 =1, ds4
+            DO s5 =1, ds5
+              DO s5 =1, ds6
+                Nin = Ngridsin(s1,s2,s3,s4,s5,s6)
+                IF (ALLOCATED(gin)) DEALLOCATE(gin)
+                ALLOCATE(gin(Nin,2))
+                IF (ALLOCATED(pin)) DEALLOCATE(pin)
+                ALLOCATE(pin(Nin))
+                gin = gridsin(s1,s2,s3,s4,s5,s6,1:Nin,:)
+                pin = percentages(s1,s2,s3,s4,s5,s6,1:Nin)
+
+                DO i1=1, di1
+                  DO i2=1, di2
+
+        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 DEFAULT
+        msg = "statisitcs '" // TRIM(stats) // "' not ready !!" // CHAR(44) // " available ones: " // &
+          "'min', 'max', 'mean', 'mean2', 'stddev', 'count'" 
+        CALL ErrMsg(msg, fname, -1)
+    END SELECT statn
+
+
+  END SUBROUTINE multi_spacewightstats_in3DRK3_slc3v
+
   SUBROUTINE multi_index_mat2DI(d1, d2, d12, mat, value, Nindices, indices)
   ! Subroutine to provide the indices of the different locations of a value inside a 2D integer matrix

trunk/tools/nc_var_tools.py

@@ -27759 +27759 @@
     fname = 'compute_slices_stats_areaweighted'
 
+    # Statistics from values within final slice
+    statns = ['min', 'max', 'mean', 'mean2', 'stddev', 'median']
+    Lstatns = {'min': 'space minimum', 'max': 'space maximum', 'mean': 'space mean', \
+      'mean2': 'space quadratic mean', 'stddev': 'spcae standard deviation',         \
+      'median': 'median', 'count': 'amount of grids'}
+    ustatns = {'min': 'vu', 'max': 'vu', 'mean': 'vu', 'mean2': 'vu', 'stddev': 'vu',\
+      'median': 'vu', 'count': 'counts'}
+
     if values == 'h':
         print fname + '_____________________________________________________________'
@@ -28526 +28534 @@
 
     # Remembering that it is python (C-like...)
-    iA = 1
-    jA = 4
-    iB = 0
-    jB = 0
     inpointsA = inpointsA-1
     inpointsB = inpointsB-1
@@ -28597 +28601 @@
     for islc in range(2,Nnewslcs):
         dn = dn + '_' + newslcvarns[iscl+1]
+
         osliceNA = onewnc.variables[newslcvarns[0] + 'Ngrid']
         osliceinA = onewnc.variables[newslcvarns[0] + 'gridin']
@@ -28605 +28610 @@
         osliceaB = onewnc.variables[newslcvarns[1] + 'gridarea']
         oslicepB = onewnc.variables[newslcvarns[1] + 'gridpercen']
-
+    
         sliceNA = osliceNA[:]
         sliceinA = osliceinA[:]
+        sliceaA = osliceaA[:]
         sliceNB = osliceNB[:]
         sliceinB = osliceinB[:]
+        sliceaB = osliceaB[:]
 
         dxA = osliceNA.shape[1]
@@ -28617 +28624 @@
         dyB = osliceNB.shape[0]
         dxyB = osliceinB.shape[1]
-
-        Srgrid = list(osliceNA.dimensions) + list(osliceNB.dimensions)
+    
+        Srgrid = list(osliceNB.dimensions) + list(osliceNA.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,    \
@@ -28635 +28642 @@
         inpointsB = inpB.transpose()
 
+        # Remembering that it is python (C-like...)
+        inpointsA = inpointsA-1
+        inpointsB = inpointsB-1
+
         maxNpointsAB = np.max(NpointsAB)
+        print "  Joined slice '" + newslcvarns[0] + "' & '" + newslcvarns[1] + "':",     \
+          NpointsAB.shape, 'maxin:', maxNpointsAB
 
         if not gen.searchInlist(onewnc.dimensions, dn+'gridin'):
@@ -28652 +28665 @@
         innewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
 
+        aanewvar = onewnc.createVariable(dn+'area','f',tuple(Srgrid),                    \
+          fill_value=gen.fillValueF)
+        basicvardef(aanewvar ,dn+'area', "area of the join slice " + newslcvarns[0] +    \
+          " & " + newslcvarns[1], vu)
+        aanewvar.setncattr('coordinates',' '.join(Srgrid[::-1]))
+
         anewvar = onewnc.createVariable(dn+'gridarea','f',tuple(Srgrid),                 \
           fill_value=gen.fillValueF)
@@ -28661 +28680 @@
           fill_value=gen.fillValueF)
         basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +                    \
-          "grids cells from .get. laying within .ref.", '1')
+          "grids cells from " + newslcvarns[0] + " laying within " + newslcvarns[1], '1')
         pnewvar.setncattr('coordinates',' '.join(Scgrid[::-1]))
 
@@ -28670 +28689 @@
                         Nin = NpointsAB[jB,iB,jA,iA]
                         innewvar[:,0:Nin,jB,iB,jA,iA] = pointsAB[:,0:Nin,jB,iB,jA,iA]
+                        slicearea = 0.
                         for iv in range(Nin):
-                            pA = oslicepA[inpA[iv],jA,iA]
-                            pB = oslicepB[inpB[iv],jB,iB]
+                            pA = oslicepA[inpointsA[iv,jA,iA],jA,iA]
+                            pB = oslicepB[inpointsB[iv,jB,iB],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
+                            ixA = sliceinA[1,inpointsA[iv,jA,iA],jA,iA]
+                            iyA = sliceinA[0,inpointsA[iv,jA,iA],jA,iA]
+                            aA = osliceaA[iyA,ixA]
+                            ixB = sliceinB[1,inpointsB[iv,jB,iB],jB,iB]
+                            iyB = sliceinB[0,inpointsB[iv,jB,iB],jB,iB]
+                            aB = osliceaB[iyB,ixB]
+                            # I do not understand why this is needed !!
+                            if aA != gen.mamat[1] and aB != gen.mamat[1]: 
+                                anewvar[jB,iB,jA,iA]= (aA*pA)*(aB*pB)
+                                slicearea = slicearea + (aA*pA)*(aB*pB)
+                            aanewvar[jB,iB,jA,iA]= slicearea
 
         onewnc.sync()
 
     slcvarns = slcvarns + newslcvarns
+
+    # Slicing with the final combined slice:', dn
+    print 'Slicing with the final combined slice:', dn
+    oNslice = onewnc.variables[dn+'Ngrid']
+    oslicein = onewnc.variables[dn+'gridin']
+    oslicearea = onewnc.variables[dn+'area']
+    oslicegridarea = onewnc.variables[dn+'gridarea']
+    oslicegridpercen = onewnc.variables[dn+'gridpercen']
+
+    sN = oNslice[:]
+    sin = oslicein[:]
+    sa = oslicearea[:]
+    sga = oslicegridarea[:]
+    sgp = oslicegridpercen[:]
+
+    dimsslice = list(oNslice.dimensions)
+    shapeslice = list(oNslice.shape)
 
     for varn in varns:
@@ -28690 +28735 @@
         vu = ovar.units
 
-        varbnds = []
-        for dn in vdimns:
-            if not gen.searchInlist(onewnc.dimensions,dn): add_dims(onc, onewnc, [dn])
-
-            odn = variables[dn]
-            dv = odn[:]
-
-            # Looking for boundaries
-            vdarattrs = odn.ncattrs()
-            if vdarattrs.has_key('bounds'):
-                print '    ' + infmsg
-                print '    ' +fname+ ": dimension variable '" +dn+ "' with bounds !!"
-                print '      bounds found:', varattrs['bounds']
-                print '      getting them to retrieve the slices'
-                bndvarns = varattrs['bounds']
-                varbnds.append(dn)
-                ovarbnds = onc.variables[varbnds]
-
-                # Slicing following boundaries, dimensions are 1D thus expand to 2D
-                #    filling as NW, NE, SE, SW
-                if gen.searchInlist(slcvarns, dn) and not vardimbndslice.has_key(dn):
-                    print infmsg
-                    print '    ' + fname + ": slicing dimension '" + dn + "' ... "
-                    varslcv = slicesinf[dn]
-                    Nslices = varslcv[0]
-                    ovardims = varslcv[1]
-                    slcvar = varslcv[2]
-
-                    # Shapping 2D slices
-                    reflat = np.zeros((Nslices), dtype=np.float)
-                    xslice2D = np.zeros((Nslices,4), dtype=np.float)
-                    xslice2D[:,0] = slcvar[1,:]
-                    xslice2D[:,1] = slcvar[1,:]
-                    xslice2D[:,2] = slcvar[0,:]
-                    xslice2D[:,3] = slcvar[0,:]
-                    yslice2D = np.zeros((Nslices,4), dtype=np.float)
-                    yslice2D[:,0] = -1.
-                    yslice2D[:,1] = 1.
-                    yslice2D[:,2] = 1.
-                    yslice2D[:,3] = -1.
-
-                    dd = len(onc.dimensions[dn])
-                    getlat = np.zeros((dd), dtype=np.float)
-                    xdimvarbnds2D = np.zeros((dd,4), dtype=np.float)
-                    xdimvarbnds2D[:,0] = ovarbnds[1,:]
-                    xdimvarbnds2D[:,1] = ovarbnds[1,:]
-                    xdimvarbnds2D[:,2] = ovarbnds[0,:]
-                    xdimvarbnds2D[:,3] = ovarbnds[0,:]
-                    ydimvarbnds2D = np.zeros((dd,4), dtype=np.float)
-                    ydimvarbnds2D[:,0] = -1.
-                    ydimvarbnds2D[:,1] = 1.
-                    ydimvarbnds2D[:,2] = 1.
-                    ydimvarbnds2D[:,3] = -1.
-
-                    reflont = slcvar.transpose()
-                    reflatt = reflat.transpose()
-                    refvarxbndst = xslice2D.transpose()
-                    refvarybndst = yslice2D.transpose()
-                    getlont = dv.transpose()
-                    getlatt = getlat.transpose()
-                    getvarxbndst = xdimvarbnds2D.transpose()
-                    getvarybndst = ydimvarbnds2D.transpose()
-
-                    Ngridsint, gridsint, percenst =                                  \
-                      fsci.module_scientific.spacepercen(xcavals=reflont,            \
-                      ycavals=reflatt, xbavals=refvarxbndst, ybavals=refvarybndst,   \
-                      xcbvals=getlont, ycbvals=getlatt, xbbvals=getvarxbndst,        \
-                      ybbvals=getvarybndst, strict=strict, dxa=Nslices, dya=1,       \
-                      navertexmax=4, dxb=dd, dyb=1, dxyb=dd, nbvertexmax=4)
-
-                    Ngridsin = Ngridsint.transpose()
-                    gridsin = gridsint.transpose()
-                    percens = percenst.transpose()
-
-                    Ngridsinmax = np.max(Ngridsin)
-
-                    # Let's avoid memory issues...
-                    lvalues = [Ngridsin, gridsin, percens]
-                    vardimbndslice[dn] = lvalues
-
-                    # Writting it to the file variables space-weight
-                    if not gen.searchInlist(onewnc.dimensions, dn+'bnds'):
-                        newdim = onewnc.createDimension(dn+'bnds',4)
-                        newdim = onewnc.createDimension(dn+'gridin',Ngridsinmax)
-                    if not gen.searchInlist(onewnc.dimensions, 'coord'):
-                        newdim = onewnc.createDimension('coord',2)
-
-                    newvar = onewnc.createVariable(dn + 'Ngrid','i',(dn))
-                    newvar[:] = Ngridsin[:]
-                    basicvardef(newvar, dn+'Ngrid', "number of grids cells from " +  \
-                      ".get. laying within .ref.", '-')
-                    newvar.setncattr('coordinates',dn)
-
-                    innewvar = onewnc.createVariable(dn+'gridin', 'i', ('coord',     \
-                      'gridin',dn))
-                    basicvardef(innewvar, dn+'gridin', "coordinates of the grids " + \
-                      "cells from .get. laying within .ref.",'-')
-                    innewvar.setncattr('coordinates',dn)
-
-                    pnewvar = onewnc.createVariable(dn+'gridpercen','f',('gridin',dn))
-                    basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +    \
-                      "grids cells from .get. laying within .ref.", '1')
-                    pnewvar.setncattr('coordinates',dn)
-
-                    for j in range(1):
-                        for i in range(dd):
-                            innewvar[:,0:Ngridsin[j,i],j,i] =                        \
-                              gridsin[:,0:Ngridsin[j,i],j,i]
-                            pnewvar[0:Ngridsin[j,i],j,i]= percens[0:Ngridsin[j,i],j,i]
-                    onewnc.sync()
+        # Looking for possible new dimensions and dim-variables
+        vardims = []
+        varshape = []
+        for dnn in vdimns:
+            vdimn = dimvars[dnn]
+            if vdimn == 'WRFtime': 
+                vdimn = 'time'
+                odimvar = onc.variables['Times']
+                timewrfv = odimvar[:]
+                tvals, urefvals = compute_WRFtime(timewrfv, refdate='19491201000000',\
+                  tunitsval='minutes')
+                tvals = np.array(tvals)
+                dimtwrf = len(tvals)
+                if gen.searchInlist(slicestatsdim, dnn):
+                    vardims.append(vdimn)
+                    varshape.append(dimtwrf)
+                if not gen.searchInlist(onewnc.dimensions,vdimn):
+                    newdim = onewnc.createDimension(vdimn, None)
+                if not gen.searchInlist(onewnc.variables,vdimn):
+                    newvar = onewnc.createVariable(vdimn, 'f', ('time'))
+                    newvar[:] = tvals
+                    basicvardef(newvar, 'time', 'Time', urefvals)
+                    newvar.setncattr('calendar','gregorian')
+
+            else:
+                if gen.searchInlist(slicestatsdim, dnn):
+                    vardims.append(vdimn)
+                    varshape.append(dimtwrf)
+                    if not gen.searchInlist(onewnc.dimensions,dnn): 
+                        add_dims(onc, onewnc, [dnn])
+                        vardims.append(dn)
+                        varshape.append(len(onc.dimensions[dn]))
+                    if not gen.searchInlist(onewnc.variables,vdimn): add_vars(onc,   \
+                      onewnc, [vdimn])
+
+        newvardims = vardims + dimsslice
+        newvarshape = varshape + shapeslice
+
+        for statn in statns:
+            Lstn = Lstatns[statn]
+            if Ustatns[statn] == 'vu': vvu = vu
+            else: 
+                if statn == 'stddev': vvu = vu + '2'
+                else: vvu = Ustatns[statn]
+            newvar= onewnc.createVariable(varn+'sliced'+statn,'f', tuple(newvardims),\
+              fill_value = gen.fillValueF)
+            basicvardef(newvar, varn+'sliced'+statn, Lstn + ' ' + varn +             \
+              ' within slice by '+ dn, vvu)
+            onewnc.sync()
+
+        # getting all the slices
+        allslices = gen.provide_slices(dimsslice, shapeslice, dimsslice)
+        print allslices
+
+        quit()
 
         # mask in var slice