Changeset 1728 in lmdz_wrf

Timestamp:

Dec 15, 2017, 4:18:27 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `econstruct_matrix_from_vector': Function to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x and y coordinates

Location:

trunk/tools

Files:

• module_scientific.f90 (modified) (2 diffs)
• nc_var_tools.py (modified) (2 diffs)

trunk/tools/module_scientific.f90

@@ -36 +36 @@
 ! read_overlap_polys_ascii: Subroutine to read from an ASCII file the associated polygons at a given time-step
 ! read_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step
+! reconstruct_matrix: Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x 
+!   and y coordinates
 ! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
 ! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a 
@@ -51 +53 @@
 
   CONTAINS
+
+  SUBROUTINE reconstruct_matrix(vectorXpos, vectorYpos, dvec, Xmin, Xmax, Ymin, Ymax, dmatx, dmaty,   &
+    matProj, maxdiff, matind, matX, matY)
+! Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x 
+!   and y coordinates
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dvec, dmatx, dmaty
+    REAL(r_k), DIMENSION(dvec), INTENT(in)               :: vectorXpos, vectorYpos
+    REAL(r_k), INTENT(in)                                :: Xmin, Xmax, Ymin, Ymax, maxdiff
+    CHARACTER(len=50), INTENT(in)                        :: matPRoj
+    INTEGER, DIMENSION(dmatx, dmaty), INTENT(out)        :: matind
+    REAL(r_k), DIMENSION(dmatx, dmaty), INTENT(out)      :: matX, matY
+
+! Local
+    INTEGER                                              :: i,j, idiff, jdiff
+    REAL(r_k)                                            :: ddx, ddy, Xpos, Ypos, mindiff
+    REAL(r_k), DIMENSION(dvec)                           :: diff
+    INTEGER, DIMENSION(1)                                :: mindiffloc
+    CHARACTER(LEN=50)                                    :: RSa, RSb
+
+!!!!!!! Variables
+! vectorXpos, vectorYpos: syncronized vectors with the x,y coordinates from which the matrix will be reconstructed
+! dvec: quantitiy of coordinates to use
+! Xmin,Xmax,Ymin,Ymax: Location range of the matrix to be reconstructed
+! maxdiff: Authorized maximum distance between matrix location and vector location
+! dmatx, dmaty: Size in grid points of the matrix to be reconstructed
+! matProj: Assumed projection of the values of the matrix
+!   'latlon': regular lat/lon projection
+! matind: matrix with the correspondant indiced from the vector of positions
+! matX, matY: matrix with the coordinates of the elements of the matrix
+
+    fname = 'reconstruct_matrix'
+
+    matind = zeroRK
+
+    Projection: SELECT CASE(TRIM(matProj))
+
+    CASE('latlon')
+      ! Resolution along matrix axes
+      ddx = (Xmax - Xmin)/(dmatx-1)
+      ddy = (Ymax - Ymin)/(dmaty-1)
+
+      ! Filling matrix positions
+      DO i=1,dmatx
+        Xpos = Xmin + ddx*(i-1)
+        DO j=1,dmaty
+          Ypos = Ymin + ddy*(j-1)
+          matX(i,j) = Xpos
+          matY(i,j) = Ypos
+
+          diff = SQRT((Xpos - vectorXpos)**2 + (Ypos - vectorYpos)**2)
+          mindiffloc = MINLOC(diff)
+          mindiff = MINVAL(diff)
+          IF (mindiff > maxdiff) THEN
+            WRITE(RSa, '(f20.10)')mindiff
+            WRITE(RSb, '(f20.10)')maxdiff
+            msg = 'Difference: ' // TRIM(RSa) // ' overpasses the maximum authorized distance: ' //   &
+              TRIM(RSb)
+            CALL ErrMsg(msg, fname, -1)
+          END IF
+
+          matind(i,j) = mindiffloc(1)
+
+        END DO
+      END DO
+
+    CASE DEFAULT
+      msg = "Projection of matrix '" // TRIM(matProj) // "' not ready " // CHAR(10) //                &
+        "  Available ones: 'latlon'"
+      CALL ErrMsg(msg, fname, -1)
+    END SELECT Projection
+
+  END SUBROUTINE reconstruct_matrix
 
   SUBROUTINE read_finaltrack_ascii(funit, dt, itrack, ftrack)

trunk/tools/nc_var_tools.py

@@ -111 +111 @@
 # pinterp: Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
 # put_variable_slice: Function to add a variable from a netcdf object to a new one following a slice
+# reconstruct_matrix_from_vector: Function to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x and y coordinates
 # remapnn: Function to remap to the nearest neightbor a variable using projection from another file
 # remapnn_old: Function to remap to the nearest neightbor a variable using projection from another fil
@@ -20693 +20694 @@
     return WRFtime, tunits
 
+def reconstruct_matrix_from_vector(values, ncfile, variable):
+    """ Function to reconstruct a 2D matrix from a pair of syncronized vectors with 
+      the positions on x and y coordinates
+        values=[varXpos]:[varYpos]:[Xmin]:[Xmax]:[Ymin]:[Ymax]:[matProj]:[maxdiff]
+          [varXpos]: variable with the vector of X positions
+          [varYpos]: variable with the vector of Y positions
+          [Xmin]: minimum X-value of the range of the matrix to be reconstructed ('all' for all range)
+          [Xmax]: maximum X-value of the range of the matrix to be reconstructed ('all' for all range)
+          [Ymin]: minimum Y-value of the range of the matrix to be reconstructed ('all' for all range)
+          [Ymax]: maximum Y-value of the range of the matrix to be reconstructed ('all' for all range)
+          [matProj]: assumed projection of the final matrix
+            'latlon': regular lon/lat projection
+          [maxdiff]: maximum allowed distance between matrix grid point and vector location
+        ncfile= netCDF file to use
+        variable= ',' list of variables to re-matrix ('all' for all variables)
+    """
+    import module_ForSci as Sci
+    fname = 'reconstruct_matrix_from_vector'
+
+    availProj = ['latlon']
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print reconstruct_matrix_from_vector.__doc__
+        quit()
+
+    expectargs = '[varXpos]:[varYpos]:[Xmin]:[Xmax]:[Ymin]:[Ymax]:[matProj]:' +      \
+      '[maxdiff]'
+    gen.check_arguments(fname, values, expectargs, ':')
+
+    varXpos = values.split(':')[0]
+    varYpos = values.split(':')[1]
+    if values.split(':')[2] != 'all':
+        Xmin = np.float(values.split(':')[2])
+        Xmax = np.float(values.split(':')[3])
+        Ymin = np.float(values.split(':')[4])
+        Ymax = np.float(values.split(':')[5])
+    else:
+        Xmin = values.split(':')[2]
+    matProj = values.split(':')[6]
+    maxdiff = np.float(values.split(':')[7])
+
+    if not gen.searchInlist(availProj, matProj):
+        print errormsg
+        print '  ' + fname + ": projection '" + matProj + "' not available !!"
+        print '    available ones:', availProj
+        uit(-1)
+
+    onc = NetCDFFile(ncfile, 'r')
+    ncvars = onc.variables.keys()
+    ncdims = onc.dimensions
+
+    # checking
+    if not gen.serarchInlist(ncvars, varXpos):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' does not have variable '" +    \
+          varXpos + "' !!"
+        print '  available ones:', ncvars
+        quit(-1)
+    if not gen.serarchInlist(ncvars, varYpos):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' does not have variable '" +    \
+          varYpos + "' !!"
+        print '  available ones:', ncvars
+        quit(-1)
+
+    if variable == 'all':
+        varns = ncvars
+    else:
+        varns = variable.split(',')
+        for vn in varns:
+            if not gen.serarchInlist(ncvars, vn):
+                print errormsg
+                print '  ' + fname + ": file '" + ncfile + "' does not have " +      \
+                  " variable '" + vn + "' !!"
+                print '  available ones:', ncvars
+                quit(-1)
+
+    # Getting values
+    oxpos = onc.variables[varXpos]
+    oypos = onc.variables[varYpos]
+
+    xposv = oxpos[:]
+    yposv = oypos[:]
+
+    gen.same_shape(xposv, yposv)
+    vecdim = oxpos.dimensions[0]
+
+    if gen.searchInlist(oxpos.ncattrs(), 'units'):
+        xunits = oxpos.getattribute('units')
+    else:
+        xunits = '-'
+    if gen.searchInlist(oypos.ncattrs(), 'units'):
+        yunits = oypos.getattribute('units')
+    else:
+        yunits = '-'
+
+    if Xmin == 'all':
+        Xmin = np.min(xposv)
+        Xmax = np.max(xposv)
+        Ymin = np.min(yposv)
+        Ymax = np.max(yposv)
+
+    # Matrix values
+    if matProj == 'latlon':
+        ddx = xposv[1]- xposv[0]
+        ddy = yposv[1]- yposv[0]
+        dimx = (Xmax - Xmin+ddx)/ddx
+        dimy = (Ymax - Ymin+ddy)/ddy
+
+    matindt, matXt, matYt = Sci.reconstruct_matrix(vectorxpos=xposv, vectorypos=yposv,  \
+      dvec=len(xposv), xmin=Xmin, xmax=Xmax, ymin=Ymin, ymax=Ymax, dmatx=dimx,          \
+      dmaty=dimy, matProj, maxdiff)
+
+    matind = matindt.transpose()
+
+    # Creation of file
+    onewnc = NetCDFFile('reconstruct_matrix.nc', 'w')
+
+    # Dimensions
+    newdim = onewnc.createDimension('x', dimx)
+    newdim = onewnc.createDimension('y', dimy)
+
+    # Variable-dimension
+    newvar = onewnc.createVariable(varXpos, 'f8', ('y', 'x'))
+    newvar[:] = matXt.transpose()
+    basicvardef(newvar, varXpos, varXpos, xunits)   
+
+    newvar = onewnc.createVariable(varYpos, 'f8', ('y', 'x'))
+    newvar[:] = matYt.transpose()
+    basicvardef(newvar, varYpos, varYpos, yunits)
+    onewnc.sync()
+
+    # Getting variables
+    for vn in varns:
+        if not onewnc.variables.has_key(vn):
+            ovar = onc.variables[vn]
+            indn = ovar.dimensions
+            vardims = []
+            shapevar = []
+            for dn in indn:
+                if not gen.searchInlist(onewnc.dimensions, dn) and dn != Xdimn and   \
+                  dn != Ydimn:
+                    if onc.dimensions[dn].isunlimited():
+                        newdim = nc.createDimension(dn, None)
+                    else:
+                        newdim = nc.createDimension(dn, len(onc.dimensions[dn]))
+                
+                if dn == vecdim: 
+                    vardims.append('y')
+                    vardims.append('x')
+                    shapevar.append(dimy)
+                    shapevar.append(dimx)
+                else: 
+                    vardims.append(dn)
+                    shapevar.append(len(onc.dimensions[dn]))
+
+            if ovar.dtype = type(int(2)):
+                newvar= onewnc.createVariable(vn, nctype(ovar.dtype), tuple(vardims),\
+                  fill_value=gen.fillValueI)
+                varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI
+            else:
+                newvar= onewnc.createVariable(vn, nctype(ovar.dtype), tuple(vardims),\
+                  fill_value=gen.fillValueF)
+                varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF
+
+            slices = gen.provide_slices(vardims, shapevar, ['x','y'])
+            # Filling variable. It would be faster if we can avoid this loop... I'm feeling lazy!
+            iix = vardims.index('x')
+            iiy = vardims.index('y')
+            for slc in slices:
+                ival = slc[iix]
+                jval = slc[iiy]
+                oldvarslice = []
+                for dn in indn:
+                    iid = vardims.index(dn)
+                    if dn == vecdimn: oldvarslice.append(matind(jval,ival))
+                    else: oldvarslice.append(slc[iid])
+
+                newvar[tuple(slc)] = onc[tuple(oldvarslice)]
+
+            # Attributes
+            for atn in ovar.ncattrs():
+                if atn != 'fill_Value':
+                    atv = ovar.getattribute(atn)
+                    set_attribute(newvar, atn, atv)
+            onewnc.sync()
+   
+    # Global attributes
+    for atn in onc.ncattrs():
+        atv = ovar.getattribute(atn)
+        set_attribute(onewnc, atn, atv)
+    onewnc.sync()
+    add_global_PyNCplot(onewnc, pyscript, fname, '0.1')
+
+    return
+
 #quit()