Changeset 1744 in lmdz_wrf

Timestamp:

Dec 16, 2017, 11:29:53 PM (8 years ago)

Author:

lfita

Message:

Improving an getting there...

Location:

trunk/tools

Files:

• ORforcing_reconstruct.py (modified) (11 diffs)
• nc_var_tools.py (modified) (1 diff)

trunk/tools/ORforcing_reconstruct.py

@@ -9 +9 @@
 # Importing generic tools file 'generic_tools.py'
 import generic_tools as gen
+import nc_var_tools as ncvar
 import subprocess as sub
+import module_ForSci as Sci
 
 filen = 'cruncep_halfdeg_1950.nc'
@@ -21 +23 @@
 
 # Range to retrieve
-Xmin=-90.
-Xmax=-7.
-Ymin=-67.
-Ymax=15.
-
-
+Xmin=-90.25
+#Xmin='all'
+Xmax=-5.25
+Ymin=-67.25
+Ymax=15.25
+
+# Variables to reconstruct
+#variable = 'all'
+variable = 'LWdown'
+
+# Resolution
+resX = 0.5
+resY = 0.5
+
+# Minimum difference from matrix to localized point
+maxdiff = 0.05
 
 # Projection of the matrix
@@ -34 +46 @@
 ## MAIN
     #######
+main = 'ORforcing_reconstruct.py'
+fname = 'ORforcing_reconstruct.py'
+
 onc = NetCDFFile(filen, 'r')
 
@@ -40 +55 @@
 ofilen = 'reconstruct_matrix.nc'
 
-ovarns = onc.variables.keys()
+ncvars = onc.variables.keys()
 
 varstocheck = [indvar, lonvar, latvar]
-for vn in vastocheck:
-    if not gen.searchInlist(ovarns, vn):
+for vn in varstocheck:
+    if not gen.searchInlist(ncvars, vn):
         print gen.errormsg
         print '  ' + main + ": file '" + filen + "' does not have variable '" + vn + \
           "' !!"
-        print '    available ones:', ovarns
+        print '    available ones:', ncvars
         quit(-1)
 
@@ -55 +70 @@
 olat = onc.variables[latvar]
 
+vecdimn = oind.dimensions[0]
+Xdimn = olon.dimensions[1]
+Ydimn = olon.dimensions[0]
+
 indv = oind[:]
 lonv = olon[:]
@@ -62 +81 @@
 dimy = lonv.shape[0]
 
-veclon = np.zeros((len(indv), dtype=np.float)
-veclat = np.zeros((len(indv), dtype=np.float)
+veclon = np.zeros((len(indv)), dtype=np.float)
+veclat = np.zeros((len(indv)), dtype=np.float)
 
 # Constuct veclon, veclat
@@ -102 +121 @@
 
 if type(Xmin) == type('2') and Xmin == 'all':
-    Xmin = np.min(xposv)
-    Xmax = np.max(xposv)
-    Ymin = np.min(yposv)
-    Ymax = np.max(yposv)
+    Xmin = np.min(lonv)
+    Xmax = np.max(lonv)
+    Ymin = np.min(latv)
+    Ymax = np.max(latv)
 
 # Matrix values
 if matProj == 'latlon':
-    dimx = (Xmax - Xmin+resX)/resX
-    dimy = (Ymax - Ymin+resY)/resY
-
-print dir(Sci)
-matindt, matXt, matYt = Sci.module_scientific.reconstruct_matrix(vectorxpos=xposv,   \
-  vectorypos=yposv, dvec=len(xposv), xmin=Xmin, xmax=Xmax, ymin=Ymin, ymax=Ymax,     \
+    dimx = int((Xmax - Xmin+resX)/resX)
+    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)
 
@@ -127 +145 @@
 
 # Variable-dimension
-newvar = onewnc.createVariable(varXpos, 'f8', ('y', 'x'))
+newvar = onewnc.createVariable('lon', 'f8', ('y', 'x'))
 newvar[:] = matXt.transpose()
-ncvar.basicvardef(newvar, varXpos, varXpos, xunits)   
-
-newvar = onewnc.createVariable(varYpos, 'f8', ('y', 'x'))
+ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_east')   
+
+newvar = onewnc.createVariable('lat', 'f8', ('y', 'x'))
 newvar[:] = matYt.transpose()
-ncvar.basicvardef(newvar, varYpos, varYpos, yunits)
+ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north')   
 onewnc.sync()
 
@@ -147 +165 @@
               dn != Ydimn:
                 if onc.dimensions[dn].isunlimited():
-                    newdim = nc.createDimension(dn, None)
+                    newdim = onewnc.createDimension(dn, None)
                 else:
-                    newdim = nc.createDimension(dn, len(onc.dimensions[dn]))
+                    newdim = onewnc.createDimension(dn, len(onc.dimensions[dn]))
                 
-            if dn == vecdim: 
+            if dn == vecdimn: 
                 vardims.append('y')
                 vardims.append('x')
@@ -161 +179 @@
 
         if ovar.dtype == type(int(2)):
-            newvar= onewnc.createVariable(vn, nctype(ovar.dtype), tuple(vardims),\
+            newvar= onewnc.createVariable(vn,ncvar.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),\
+        elif ovar.dtype == type(np.int32(2)):
+            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
+              fill_value=gen.fillValueI)
+            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI
+        elif ovar.dtype == type(np.int64(2)):
+            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
+              fill_value=gen.fillValueI)
+            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI
+        elif ovar.dtype == type(np.float(2.)):
+            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
               fill_value=gen.fillValueF)
             varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF
+        elif ovar.dtype == type(np.float32(2.)):
+            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
+              fill_value=gen.fillValueF)
+            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF
+        else:
+            print gen.errormsg
+            print '  ' + fname + ': variable type:', ovar.dtype, ' not ready !!'
+            quit(-1)
 
         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)]
+        if not gen.searchInlist(vardims,'x') and not gen.searchInlist(vardims,'y'):
+            newvar[:] = ovar[:] 
+        else:
+            iix = vardims.index('x')
+            iiy = vardims.index('y')
+            for slc in slices:
+                ival = slc[iix]
+                jval = slc[iiy]
+                # Avoiding not found grid points
+                if matind[jval,ival] != -1:
+                    oldvarslice = []
+                    for dn in indn:
+                        if dn == vecdimn: oldvarslice.append(matind[jval,ival])
+                        else: 
+                            iid = vardims.index(dn)
+                            oldvarslice.append(slc[iid])
+
+                    print vn, 'Slice new:', slc, 'old:', oldvarslice, 'shape:', newvar.shape, 'oldshape:', ovar.shape
+                    newvar[tuple(slc)] = ovar[tuple(oldvarslice)]
+                else:
+                    if type(newvar.dtype) == type(float(2.)) or                      \
+                      type(newvar.dtype) == type(np.float(2.)) or                    \
+                      type(newvar.dtype) == type(np.float32(2)) or                   \
+                      type(newvar.dtype) == type(np.float64(2)):
+                        newvar[tuple(slc)] = gen.fillValuerF
+                    else:
+                        newvar[tuple(slc)] = gen.fillValueI
 
         # Attributes
         for atn in ovar.ncattrs():
-            if atn != 'fill_Value':
+            if atn != 'fill_Value' and atn != 'units':
                 atv = ovar.getncattr(atn)
-                set_attribute(newvar, atn, atv)
+                ncvar.set_attribute(newvar, atn, atv)
+        ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
         onewnc.sync()
    
@@ -194 +244 @@
 for atn in onc.ncattrs():
     atv = ovar.getncattr(atn)
-    set_attribute(onewnc, atn, atv)
+    ncvar.set_attribute(onewnc, atn, atv)
 onewnc.sync()
-add_global_PyNCplot(onewnc, pyscript, fname, '0.1')
+ncvar.add_global_PyNCplot(onewnc, pyscript, fname, '0.1')
 
 onc.close()

trunk/tools/nc_var_tools.py

@@ -20783 +20783 @@
     yposv = oypos[:]
 
+    if len(xposv.shape) != 1:
+        print errormsg
+        print '  ' + fname + ": vector with X locations (from variable '" + varXpos +\
+          "') has a rank larger than 1!!"
+        print '    current shape:', xposv.shape
+        quit(-1)
+
+    if len(yposv.shape) != 1:
+        print errormsg
+        print '  ' + fname + ": vector with Y locations (from variable '" + varXpos +\
+          "') has a rank larger than 1!!"
+        print '    current shape:', yposv.shape
+        quit(-1)
+
     gen.same_shape(xposv, yposv)
     vecdim = oxpos.dimensions[0]