Changeset 1368 in lmdz_wrf

Timestamp:

Dec 5, 2016, 6:10:42 PM (9 years ago)

Author:

lfita

Message:

Improving `SpatialWeightedMean?'

File:

• trunk/tools/nc_var_tools.py (modified) (8 diffs)

trunk/tools/nc_var_tools.py

@@ -159 +159 @@
 warnmsg = 'WARNING -- warning -- WARNING -- warning'
 
+## Constants
+# Radius of the Earth (m)
+EarthR = 6378000.
+
 class genericNCvariable(object):
     """ Class to fake a netCDF varibale
@@ -14105 +14109 @@
     if oncf.variables.has_key(varname):
         print warnmsg
-        print '  ' + fname + ": netCDF object already has variable '" + ovar.name +  \
+        print '  ' + fname + ": netCDF object already has variable '" + varname +   \
           "' !!"
         print '    doing noting'
@@ -14144 +14148 @@
         [xdimname]: name of the dimension for the x-axis
         [ydimname]: name of the dimension for the y-axis
-        [addvars]: ':', separetd list of name of variables to add to the output file 
+        [addvars]: ':', separetd list of name of variables to add to the output file ('None' for any)
       filen= name of the netCDF file
       varn= name of the variable
@@ -14407 +14411 @@
         ilatv = iolat[:]
 
-        lonv, latv = lonlat2D(ilonv, ilatv)
-
-        ivarwgtv = np.abs(np.cos(latv*np.pi/180.))
+        lonv, latv = gen.lonlat2D(ilonv, ilatv)
+
+        # Direct and wrong
+        #ivarwgtv = np.abs(np.cos(latv*np.pi/180.))
+
+        # Using shoelace formula
+        dx = lonv.shape[1]
+        dy = lonv.shape[0]
+
+        dlon = np.zeros((lonv.shape), dtype=np.float)
+        dlat = np.zeros((latv.shape), dtype=np.float)
+        ivarwgtv = np.zeros((latv.shape), dtype=np.float)
+
+        dlon[:,0:dx-1] = lonv[:,1:dx] - lonv[:,0:dx-1]
+        dlat[0:dy-1,:] = latv[1:dy,:] - latv[0:dy-1,:]
+        dlon[:,dx-1] = lonv[:,dx-1] - lonv[:,dx-2]
+        dlat[dy-1,:] = latv[dy-1,:] - latv[dy-2,:]
+
+        coslat = np.abs(np.cos(latv*np.pi/180.))
+        dsx = EarthR*coslat*dlon
+        dsy = EarthR*dlat
+        for j in range(dy-1):
+            for i in range(dx-1):
+#                xvals = [dsx[j,i], dsx[j+1,i], dsx[j+1,i+1], dsx[j,i+1]]
+#                yvals = [dsy[j,i], dsy[j+1,i], dsy[j+1,i+1], dsy[j,i+1]]
+                xvals = [0., 0., dsx[j+1,i+1], dsx[j,i+1]]
+                yvals = [0., dsy[j+1,i], dsy[j+1,i+1], 0.]
+                ivarwgtv[j,i] = gen.PolyArea(xvals, yvals)
+        j=dy-1
+        for i in range(dx-1):
+            xvals = [0., 0., dsx[j-1,i+1], dsx[j,i-1]]
+            yvals = [0., dsy[j-1,i], dsy[j-1,i+1], 0.]
+            ivarwgtv[j,i] = gen.PolyArea(xvals, yvals)
+        i=dx-1
+        for j in range(dy-1):
+            xvals = [0., 0., dsx[j+1,i-1], dsx[j,i-1]]
+            yvals = [0., dsy[j+1,i], dsy[j+1,i-1], 0.]
+            ivarwgtv[j,i] = gen.PolyArea(xvals, yvals)
+
+        i=dx-1
+        j=dy-1
+        xvals = [0., 0., dsx[j-1,i-1], dsx[j,i-1]]
+        yvals = [0., dsy[j-1,i], dsy[j-1,i-1], 0.]
+        ivarwgtv[j,i] = gen.PolyArea(xvals, yvals)
+
         TOTsumwgt = np.sum(ivarwgtv)
 
@@ -14449 +14495 @@
                         newsumvals[id1,id2,id3] = np.sum(vals*ivarwgtv)
                         newvals[id1,id2,id3] = newsumvals[id1,id2,id3] / TOTsumwgt
-        outweightvals = ivarwgtv
+        outweightvals = ivarwgtv/TOTsumwgt
 
 # Writting output file
 ##
+    print 'TOTAL sunm of weights:', np.sum(outweightvals)
     onewnc = NetCDFFile(ofile, 'w')
 
@@ -14481 +14528 @@
     onewnc.sync()
 
-    print 'newsumvals:',newsumvals,'newvals:',newvals,'TOTsumwgt:',TOTsumwgt
-
 # Spatial weight
 ##
@@ -14494 +14539 @@
     onewnc.sync()
 
+    newvar = onewnc.createVariable('dlon', 'f4', tuple([ydimname, xdimname]), \
+      fill_value=gen.fillValueF)
+    basicvardef(newvar, 'dlon', 'dlon', 'degrees_East')
+    newvar[:] = dlon[:]
+    newvar = onewnc.createVariable('dlat', 'f4', tuple([ydimname, xdimname]), \
+      fill_value=gen.fillValueF)
+    basicvardef(newvar, 'dlat', 'dlat', 'degrees_North')
+    newvar[:] = dlat[:]
+    newvar = onewnc.createVariable('coslat', 'f4', tuple([ydimname, xdimname]), \
+      fill_value=gen.fillValueF)
+    basicvardef(newvar, 'coslat', 'latitude cosinus', '-')
+    newvar[:] = coslat[:]
+    newvar = onewnc.createVariable('dsx', 'f4', tuple([ydimname, xdimname]), \
+      fill_value=gen.fillValueF)
+    basicvardef(newvar, 'dsx', 'x distance', 'm')
+    newvar[:] = dsx[:]
+    newvar = onewnc.createVariable('dsy', 'f4', tuple([ydimname, xdimname]), \
+      fill_value=gen.fillValueF)
+    basicvardef(newvar, 'dsy', 'y distance', 'm')
+    newvar[:] = dsy[:]
+    onewnc.sync()
+
 # Additional variables
 ##
-    for vn in addvars.split(':'):
-        if not onc.variables.has_key(vn):
-            print errormsg
-            print '  ' + fname + ": file '" + filen + "' does not have variable '" + \
-              vn + "' !!"
-            quit(-1)
-        ovar_onc(onc, onc.variables[vn], onewnc)
-    onewnc.sync()
+    if addvars != 'None':
+        for vn in addvars.split(':'):
+            if not onc.variables.has_key(vn):
+                print errormsg
+                print '  ' + fname + ": file '" + filen + "' does not have variable '" + \
+                  vn + "' !!"
+                quit(-1)
+            ovar_onc(onc, onc.variables[vn], onewnc)
+        onewnc.sync()
 
 # Global attributes
@@ -14654 +14722 @@
 
     # Given constants
-    EarthR = 6378000.
     RadDeg = 180./np.pi