Changeset 1782 in lmdz_wrf

Timestamp:

Feb 26, 2018, 6:33:54 PM (7 years ago)

Author:

lfita

Message:

New version for matrix reconstruction taking as reference vector values

Location:

trunk/tools

Files:

• ORforcing_reconstruct.py (modified) (10 diffs)
• module_scientific.f90 (modified) (5 diffs)

trunk/tools/ORforcing_reconstruct.py

@@ -15 +15 @@
 
 parser = OptionParser()
+parser.add_option("-F", "--filename", dest="fn", help="name of files (overwrites -f option)", \
+  metavar="VALUE")
 parser.add_option("-f", "--fileHEader", dest="fh", help="header of files",           \
   metavar="VALUE")
+parser.add_option("-i", "--indices", dest="indn", help="name of the variable indices", \
+  metavar="VALUE")
 parser.add_option("-L", "--latitude", dest="latn", help="name of the variable latitiude", \
   metavar="VALUE")
 parser.add_option("-l", "--longitude", dest="lonn", help="name of the variable longitiude", \
   metavar="VALUE")
+parser.add_option("-t", "--TransposeVariables", dest="tvars", help="whether variables should be trasposed", \
+  metavar="VALUE")
 parser.add_option("-v", "--Variables", dest="varns", help="',' separated list of variables", \
   metavar="VALUE")
@@ -30 +36 @@
 ####### ###### ##### #### ### ## #
 
-filen = opts.fh + opts.year + '.nc'
+if opts.fn is None:
+    filen = opts.fh + opts.year + '.nc'
+else:
+    filen = opts.fn
 
 # Variable whcih provides the indices of a 1D vector from the dimy, dimx space
-indvar = 'land'
+indvar = opts.indn 
 
 # 2D longitude, latitude matrices
@@ -41 +50 @@
 # Range to retrieve
 Xmin=-90.25
-#Xmin='all'
+Xmin='all'
 Xmax=-33.25
 Ymin=-67.25
@@ -91 +100 @@
 latv = olat[:]
 
-vecdimn = oind.dimensions[0]
 if len(olon.dimensions) == 2:
+    dimx = lonv.shape[1]
+    dimy = lonv.shape[0]
+    veclon = lonv.reshape(dimx*dimy)
+    veclat = latv.reshape(dimx*dimy)
     Xdimn = olon.dimensions[1]
     Ydimn = olon.dimensions[0]
-    dimx = lonv.shape[1]
-    dimy = lonv.shape[0]
-elif len(olon.dimensions) == 1:
-    Xdimn = olon.dimensions[0]
-    Ydimn = olat.dimensions[0]
-    dimx = lonv.shape[0]
-    dimy = latv.shape[0]
-    lonv, latv = np.meshgrid(lonv, latv)
-
-# Fortran indices, first 1
-indv = indv - 1
-
-veclon = np.zeros((len(indv)), dtype=np.float)
-veclat = np.zeros((len(indv)), dtype=np.float)
-
-# Construct veclon, veclat
-for iid in range(len(indv)):
-    iy = int((indv[iid]-1)/dimx)
-    ix = int(indv[iid] - iy*dimx - 1)
-    veclon[iid] = lonv[iy,ix]
-    veclat[iid] = latv[iy,ix]
+else:
+    veclon = lonv.copy()
+    veclat = latv.copy()
+    Xdimn = 'x'
+    Ydimn = 'y'
+
+vecdimn = oind.dimensions[0]
 
 if not gen.searchInlist(availProj, matProj):
@@ -157 +155 @@
     dimy = int((Ymax - Ymin+resY)/resY)
 
-matindt, matXt, matYt = Sci.module_scientific.reconstruct_matrix(vectorxpos=veclon,  \
-  vectorypos=veclat, dvec=veclon.shape[0], xmin=Xmin, xmax=Xmax, ymin=Ymin,ymax=Ymax,\
-  dmatx=dimx, dmaty=dimy, matproj=matProj, maxdiff=maxdiff)
+print 'Xmin:', Xmin, 'Xmax:', Xmax, 'Ymin:', Ymin, 'Ymax:', Ymax, 'maxdiff:', maxdiff
+
+matindt, matXt, matYt, matdifft = Sci.module_scientific.reconstruct_matrix(          \
+  vectorxpos=veclon, vectorypos=veclat, dvec=veclon.shape[0], xmin=Xmin, xmax=Xmax,  \
+  ymin=Ymin,ymax=Ymax, dmatx=dimx, dmaty=dimy, matproj=matProj, maxdiff=maxdiff)
 
 matind = matindt.transpose()
+Nfound = np.sum(matind != -1)
+Nstations = veclon.shape[0]
+print '  Nfound:', Nfound, ' number of stations:', Nstations
+
+if Nfound*1. / Nstations < 0.8:
+    print gen.errormsg
+    print '  '+main + ': only ', '{:.2f}'.format(Nfound*100./Nstations),\
+      '% of points ' + 'have been found !!'
+    print '    this is not enough. Something must went wrong!'
+    print '    Longitudes Latitudes _______'
+    for i in range(veclon.shape[0]):
+        print '    ', veclon[i], veclat[i]
+    dx2 = dimx/2
+    dy2 = dimy/2
+    print '    dx2, dy2 -/+ 5 fraction of destiny longitudes _______'
+    for j in range(-5,5):
+        print matXt[dx2-5:dx2+5,dy2+j]
+    print '    dx2, dy2 -/+ 5 fraction of destiny latitudes _______'
+    for j in range(-5,5):
+        print matYt[dx2-5:dx2+5,dy2+j]
+    print '    dx2, dy2 -/+ 5 fraction of indices equivalency _______'
+    for j in range(-5,5):
+        print matindt[dx2-5:dx2+5,dy2+j]
+    print '     min distance lon(dy2,dx2)=', matXt[dx2,dy2], ':',                    \
+      np.min(veclon - matXt[dx2,dy2])
+    print '     min distance lat(dy2,dx2)=', matYt[dx2,dy2], ':',                    \
+      np.min(veclat - matYt[dx2,dy2])
+    print '     longitude borders:', Xmin, Xmax, 'dX:', resX
+    print '     latitude borders:', Ymin, Ymax, 'dY:', resY
+    #quit(-1)
+
 # Fortran like, First 1
 matind = np.where(matind != -1, matind - 1, matind)
@@ -189 +220 @@
 ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
 
+# Variable differences
+newvar = onewnc.createVariable('vec1D_matdiff', 'i', ('y', 'x'), fill_value=-1)
+newvar[:] = matdifft.transpose()
+ncvar.basicvardef(newvar,'vec1D_matdiff', 'matrix differences respect 1D ', 'degrees')
+ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
+
 # Looking for equivalencies in the 1D vector
 matlonlat = matind.copy()
@@ -207 +244 @@
     if not onewnc.variables.has_key(vn):
         ovar = onc.variables[vn]
-        indn = ovar.dimensions
+        if gen.Str_Bool(opts.tvars):
+            indn0 = ovar.dimensions
+            indn = list(indn0)[::-1]
+        else:
+            indn = ovar.dimensions
         vardims = []
         shapevar = []
@@ -252 +293 @@
             quit(-1)
 
-        print '  reconstructing:', vn, '...'
+        print '  reconstructing:', vn, ' shape:', newvar.shape, '...'
         # Filling variable. It would be faster if we can avoid this loop... I'm feeling lazy!
         if not gen.searchInlist(vardims,'x') and not gen.searchInlist(vardims,'y'):
-            newvar[:] = ovar[:] 
+            if gen.Str_Bool(opts.tvars):
+                newvar[:] = ovar[:].transpose() 
+            else:
+                newvar[:] = ovar[:] 
         else:
-            ovart = ovar[:].transpose()
+            if gen.Str_Bool(opts.tvars):
+                ovart = ovar[:]
+            else:
+                ovart = ovar[:].transpose()
+            print '  Lluis shapes ovart:', ovart.shape, 'newvar:', newvar.shape
             if newvar.dtype == type(float(2.)) or newvar.dtype == type(np.float(2.)) \
               or newvar.dtype == type(np.float32(2)) or                              \
@@ -280 +328 @@
 # Global attributes
 for atn in onc.ncattrs():
-    atv = ovar.getncattr(atn)
+    atv = onc.getncattr(atn)
     ncvar.set_attribute(onewnc, atn, atv)
 onewnc.sync()
@@ -287 +335 @@
 
 # Reconstructing times
-otime = onewnc.variables['time']
-ncvar.set_attribute(otime, 'units', 'seconds since ' + opts.year + '-01-01 00:00:00')
+#otime = onewnc.variables['time']
+#ncvar.set_attribute(otime, 'units', 'seconds since ' + opts.year + '-01-01 00:00:00')
 
 onewnc.sync()

trunk/tools/module_scientific.f90

@@ -127 +127 @@
 
   SUBROUTINE reconstruct_matrix(vectorXpos, vectorYpos, dvec, Xmin, Xmax, Ymin, Ymax, dmatx, dmaty,   &
-    matProj, maxdiff, matind, matX, matY)
+    matProj, maxdiff, matind, matX, matY, matdiff)
 ! Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x 
 !   and y coordinates
@@ -138 +138 @@
     CHARACTER(len=50), INTENT(in)                        :: matPRoj
     INTEGER, DIMENSION(dmatx, dmaty), INTENT(out)        :: matind
-    REAL(r_k), DIMENSION(dmatx, dmaty), INTENT(out)      :: matX, matY
+    REAL(r_k), DIMENSION(dmatx, dmaty), INTENT(out)      :: matX, matY, matdiff
 
 ! Local
-    INTEGER                                              :: i,j, idiff, jdiff
+    INTEGER                                              :: i,j,iv, idiff, jdiff
     REAL(r_k)                                            :: ddx, ddy, Xpos, Ypos, mindiff
-    REAL(r_k), DIMENSION(dvec)                           :: diff
-    INTEGER, DIMENSION(1)                                :: mindiffloc
+    REAL(r_k), DIMENSION(dmatx,dmaty)                    :: diff
+    REAL(r_k)                                            :: nXpos, xXpos, nYpos, xYpos
+    INTEGER, DIMENSION(2)                                :: mindiffloc
     CHARACTER(LEN=50)                                    :: RSa, RSb
 
@@ -157 +158 @@
 ! matind: matrix with the correspondant indiced from the vector of positions
 ! matX, matY: matrix with the coordinates of the elements of the matrix
+! matdiff: matrix with the differences at each grid point
 
     fname = 'reconstruct_matrix'
 
-    matind = zeroRK
+    matind = -oneRK
+    matdiff = -oneRK
+
+    nXpos = MINVAL(vectorXpos)
+    xXpos = MAXVAL(vectorXpos)
+    nYpos = MINVAL(vectorYpos)
+    xYpos = MAXVAL(vectorYpos)
+    PRINT *, 'Limits of the positions to localize in a matrix: longitudes:', nXpos,  &
+     ' ,',xXpos, ' latitudes:', nYpos, ' ,', xYpos
 
     Projection: SELECT CASE(TRIM(matProj))
@@ -176 +186 @@
           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
-            !Assuming no values
-            matind(i,j) = -1
-            !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)
-          ELSE
-            matind(i,j) = mindiffloc(1)
-          END IF
-
         END DO
       END DO
@@ -200 +194 @@
       CALL ErrMsg(msg, fname, -1)
     END SELECT Projection
+
+    ! Let's do it properly...
+    DO iv=1,dvec
+      diff = SQRT((matX - vectorXpos(iv))**2 + (matY - vectorYpos(iv))**2)
+      mindiffloc = MINLOC(diff)
+      mindiff = MINVAL(diff)
+      IF (mindiff > maxdiff) THEN
+        PRINT *,'  vectorXpos:', vectorXpos(iv), 'vectorYpos:', vectorYpos(iv)
+        !PRINT *,'  Xpos:', Xpos, 'Ypos:', Ypos
+        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)
+      ELSE
+        i = mindiffloc(1)
+        j = mindiffloc(2)
+        matind(i,j) = iv
+        matdiff(i,j) = mindiff
+      END IF
+    END DO
+
+    RETURN
 
   END SUBROUTINE reconstruct_matrix