Changeset 1740 in lmdz_wrf

Timestamp:

Dec 16, 2017, 8:01:31 PM (7 years ago)

Author:

lfita

Message:

Completting script

File:

• trunk/tools/ORforcing_reconstruct.py (modified) (2 diffs)

trunk/tools/ORforcing_reconstruct.py

@@ -26 +26 @@
 Ymax=15.
 
+
+
+# Projection of the matrix
+matProj = 'latlon'
+
 #######    #######
 ## MAIN
     #######
 onc = NetCDFFile(filen, 'r')
+
+availProj = ['latlon']
+
 ofilen = 'reconstruct_matrix.nc'
 
@@ -61 +69 @@
     iy = int(indv[iid]/dimx)
     ix = indv[iid] - iy*dimx
-    print iid, 'Lluis: indv:', indv[iid], 'iy:', iy, 'ix:', ix
+    veclon[iid] = lonv[iy,ix]
+    veclat[iid] = latv[iy,ix]
+
+if not gen.searchInlist(availProj, matProj):
+    print errormsg
+    print '  ' + fname + ": projection '" + matProj + "' not available !!"
+    print '    available ones:', availProj
     quit(-1)
+
+ncdims = onc.dimensions
+
+if variable == 'all':
+    varns = ncvars
+else:
+    varns = variable.split(',')
+    for vn in varns:
+        if not gen.searchInlist(ncvars, vn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have " +          \
+              " variable '" + vn + "' !!"
+            print '  available ones:', ncvars
+            quit(-1)
+
+if gen.searchInlist(olon.ncattrs(), 'units'):
+    xunits = olon.getncattr('units')
+else:
+    xunits = '-'
+if gen.searchInlist(olat.ncattrs(), 'units'):
+    yunits = olat.getncattr('units')
+else:
+   yunits = '-'
+
+if type(Xmin) == type('2') and Xmin == 'all':
+    Xmin = np.min(xposv)
+    Xmax = np.max(xposv)
+    Ymin = np.min(yposv)
+    Ymax = np.max(yposv)
+
+# 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,     \
+  dmatx=dimx, dmaty=dimy, matproj=matProj, maxdiff=maxdiff)
+
+matind = matindt.transpose()
+
+# Creation of file
+onewnc = NetCDFFile(ofilen, 'w')
+
+# Dimensions
+newdim = onewnc.createDimension('x', dimx)
+newdim = onewnc.createDimension('y', dimy)
+
+# Variable-dimension
+newvar = onewnc.createVariable(varXpos, 'f8', ('y', 'x'))
+newvar[:] = matXt.transpose()
+ncvar.basicvardef(newvar, varXpos, varXpos, xunits)   
+
+newvar = onewnc.createVariable(varYpos, 'f8', ('y', 'x'))
+newvar[:] = matYt.transpose()
+ncvar.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.getncattr(atn)
+                set_attribute(newvar, atn, atv)
+        onewnc.sync()
+   
+# Global attributes
+for atn in onc.ncattrs():
+    atv = ovar.getncattr(atn)
+    set_attribute(onewnc, atn, atv)
+onewnc.sync()
+add_global_PyNCplot(onewnc, pyscript, fname, '0.1')
+
+onc.close()
+onewnc.sync()
+onewnc.close()
+
+print fname + ": Successful writing of file 'reconstruct_matrix.nc' !!"