Changeset 2181 in lmdz_wrf for trunk/tools

Timestamp:

Oct 12, 2018, 7:32:27 PM (7 years ago)

Author:

lfita

Message:

ADding:

• `topography_shadow': Function to compute shadow topography

File:

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

trunk/tools/drawing_tools.py

@@ -198 +198 @@
 # plot_ZQradii: Function to plot following radial averages only at exact grid poins
 # transform: Function to transform the values and the axes
+# topography_shadow: Function to compute shaow topography
 
 # From nc_var_tools.py
@@ -13827 +13828 @@
     return mapenv, merids, parals, xmap, ymap
 
+def topography_shadow(geogf, hgtvn, lonvn, latvn, mapfmtS)
+    """ Function to compute shaow topography
+      geogf: file with the topography values
+      hgtvn: name of the variable with the topography
+      lonvn: name of the variable with the longitudes
+      latvn: name of the variable with the latitudes
+      mapfmt: [cbarmap]|[mapkind]|[lonlatbox] values for the map 
+        cbarmap: name of the colorbar for the map
+        mapkind: kind of map to use in the plot
+          'direct': values as they come
+          'shadow',[pos],[enhance]: pseudo-shadding from a given location of the sun
+            [pos]: 'N', 'NW'
+            [enhance]: enhance factor for the shading ('auto' for 1./5.)
+        lonlatbox: [lonSW],[latSW],[lonNE],[latNE] to plot only a lon,lat box
+          'full': for the whole domain
+    """
+    fname = 'topography_shadow'
+
+    # topography
+    if not os.paths.isfile(geogf):
+       print errormsg
+       print '  ' + fname + ": topography file '" + geogf + "' does not exist !!"
+       quit(-1)
+
+    ong = NetCDFFile(geogf,'r')
+
+    if not ong.variables.has_key(hgtvn):
+        print errormsg
+        print '  ' + fname + ": topography file '" + geogf + "' does not have " +    \
+          "topgraphy as variable named: '" + hgtvn + "' !!"
+        varns = ong.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+    if not ong.variables.has_key(lonvn):
+        print errormsg
+        print '  ' + fname + ": topography file '" + geogf + "' does not have " +    \
+          "longitudes as variable named: '" + lonvn + "' !!"
+        varns = ong.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+    if not ong.variables.has_key(hgtvn):
+        print errormsg
+        print '  ' + fname + ": topography file '" + geogf + "' does not have " +    \
+          "latitudes as variable named: '" + latvn + "' !!"
+        varns = ong.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+
+    otopo = ong.variables[hgtvn]
+    olon = ong.variables[lonvn]
+    olat = ong.variables[latvn]
+
+    if len(otopo.shape) == 3: valsmap = otopo[0,:,:]
+    elif len(otopo.shape) == 2: valsmap = otopo[:]
+    else:
+        print errormsg
+        print '  ' + fname + ": wring rank for topography variable '" + hgtvn + "' !!"
+        print '    provided shape:', otopo.shape
+        quit(-1)
+
+    lon,lat = gen.lonlat2D(olon[:], olat[:])
+    ong.close()
+
+    dx = valsmap.shape[1]
+    dy = valsmap.shape[0]
+
+    mapfmt = mapfmtS.split
+
+    mapkind = mapfmt[1]
+    if mapfmt[2] == 'full':
+        lonlatbox = None
+    elif len(mapfmt[2].split(',')) == 4:
+        lonlatbox = []
+        for lL in mapfmt[2].split(','):
+            lonlatbox.append(np.float(lL))
+    else:
+        print errormsg
+        print '  ' + fname + ": lonlatbox value '" + mapfmt[2] + "' not ready !!"
+        print '    for a lonlatbox 4 values are required. Values passed:',           \
+          mapfmt[2].splot(',')
+        quit(-1)
+    
+    if mapkind == 'direct':
+        mapv = valsmap[:]
+    elif mapkind[0:6] == 'shadow':
+        mapv = np.zeros((dy,dx), dtype=np.float)
+        sunpos = mapkind.split(',')[1]
+        enhance = gen.auto_val(mapkind.split(',')[2], 1./5.)
+        if sunpos == 'N':
+            ## Linear shadow from N
+            mapv[0:dy-1,:] = valsmap[1:dy,:] - valsmap[0:dy-1,:]
+        elif sunpos == 'NW':
+            ## Diagonal shadow from NW
+            for i in range(dx-1):
+                for j in range(dy-1):
+                    if (i-1 >= 0) and (j+1 <= dy-1): 
+                        mapv[j,i] = valsmap[j+1,i-1] - valsmap[j,i]
+        else:
+            print errormsg
+            print '  ' + fname + ": sun location for map sahdding '" + sunpos +      \
+              "' not ready !!"
+            print '    available ones:', availsunpos
+            quit(-1)
+        xmapv = np.max(mapv)*enhance
+        mapv = np.where(mapv > xmapv, xmapv, mapv)
+        mapv = ma.masked_less(mapv, 0.)
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of map '" + mapkind + "' not ready !!"
+        print '    available ones:', availmapkind
+        quit(-1)
+
+    return mapv