Changeset 1782 in lmdz_wrf for trunk/tools/ORforcing_reconstruct.py

Timestamp:

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

Author:

lfita

Message:

New version for matrix reconstruction taking as reference vector values

File:

• trunk/tools/ORforcing_reconstruct.py (modified) (10 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()