Changeset 2751 in lmdz_wrf for trunk

Timestamp:

Nov 12, 2019, 1:20:22 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `draw_geogrid_landuse': plotting land-use data from geo_em.d[nn].nc WPS files

Location:

trunk/tools

Files:

• drawing.py (modified) (5 diffs)
• drawing_tools.py (modified) (2 diffs)

trunk/tools/drawing.py

@@ -69 +69 @@
 ## e.g. # $pyHOME/drawing.py -o draw_river_pattern -S 'ORCHIDEEdef:ORCHIDEEdef:auto:cyl,l:None:-100.|-65.|-20.|15.:False:True:South|America|ORCHIDEE|routing|0.5|deg:pdf:0|10|1:True' -v Parana,Amazon
 ## e.g. # python ${pyHOME}/drawing.py -o draw_topofix_geogrid_boxes -f geo_em.d01.nc,geo_em.d02.nc -S 'None:RELAMPAGO!explicit!convection!configuration:pdf:cyl,i:d01$_{8k}$,d02$_{1.6k}$:0|10|1:auto:0.17|horizontal|8|10|-45.|0.8|0.15|#CCCCFF:auto:auto:auto:auto:auto:True'
+## e.g. # drawing.py -o draw_geogrid_landuse -f geo_em.d01.nc -S None:PICT!AMBA-NO!WRF!domains!land-use!d0!(15!km):pdf-png:cyl,l:usgs|0.9:auto:auto:auto:auto:auto:True
 
 #######
@@ -94 +95 @@
 # draw_cycle: Function to plot a variale with a circular cycle
 # draw_ensembles_time: Function to plot an ensembles of data following an axis-time
+# draw_geogrid_landuse: plotting land-use data from geo_em.d[nn].nc WPS files
 # draw_lines: Function to draw different lines at the same time from different files
 # draw_lines_time: Function to draw different lines at the same time from different files with times
@@ -140 +142 @@
   'draw_2D_shad_line_time', 'draw_bar', 'draw_bar_line', 'draw_bar_line_time',       \
   'draw_bar_time', 'draw_barbs', 'draw_basins', 'draw_cycle', 'draw_ensembles_time', \
-  'draw_2lines', 'draw_2lines_time',                                                 \
+  'draw_2lines', 'draw_2lines_time', 'draw_geogrid_landuse',                        \
   'draw_lines', 'draw_lines_time',                                                   \
   'draw_multi_2D_shad', 'draw_multi_SkewT', 'draw_multiWindRose',                    \
@@ -11479 +11481 @@
     return
 
+def draw_geogrid_landuse(ncfile, values):
+    """ plotting land-use data from geo_em.d[nn].nc WPS files
+      ncfile= geo_em.d[nn].nc file to use
+      values=[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[ColorLeg]:[mervals]:[parvals]:
+          [coastlvals]:[countrylvals]:[stslvals]:[close]
+        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] 
+          or None
+        title: title of the graph ('!' for spaces)
+        graphic_kind: kind of figure (jpg, pdf, png)
+        mapvalues: map characteristics [proj],[res]
+          see full documentation: http://matplotlib.org/basemap/
+          [proj]: projection
+            * 'cyl', cilindric
+            * 'lcc', lambert conformal
+          [res]: resolution:
+            * 'c', crude
+            * 'l', low
+            * 'i', intermediate
+            * 'h', high
+            * 'f', full
+        ColorLeg= [ijstart]|[lucolors] values for the legend of colors of the land-use
+          ([0.85, 'usgscolors'], default)
+          ijstart: initial position as fraction of figure on the oposite axis of 
+            orientation
+          lucolors: name of the dictionary with the equivalences between land-use and
+            color
+        mervals: [fontsize]|[color line]|[labels rotation]|[line width] values for 
+          the meridians  ([8,'#AAAAAA', 0, 0.5], auto)
+        parvals: [fontsize]|[color line]|[labels rotation]|[line width] values for 
+          the parallels ([8,'#AAAAAA', 0, 0.5], auto)
+        coastlvals: [line width|color line] values for the coastline (None for any, 
+          [0.25,'#161616'], auto)
+        countrylvals: [line width|color line] values for the countries (None for any,
+          [0.25, '#161616'], auto)
+        stslvals: [line width|color line] values for the states (None for any, None,
+          [0.25, '#080808'], auto) 
+        close: Whether figure should be finished or not
+    """
+    import matplotlib.pyplot as plt
+
+    fname = 'draw_geogrid_landuse'
+    lucolorsavail = ['modifmodis', 'ugs']
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_geogrid_landuse.__doc__
+        quit()
+
+    expectargs = '[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[ColorLeg]:' +        \
+      '[mervals]:[parvals]:[coastlvals]:[countrylvals]: [stslvals]:[close]'
+    drw.check_arguments(fname,values,expectargs,':')
+
+    lonlatLS = values.split(':')[0]
+    lonlatLv = lonlatLS.split(',')[0]
+
+    if lonlatLv == 'None':
+        lonlatL = None
+    else:
+        lonlatL = np.zeros((4), dtype=np.float)
+        lonlatL[0] = np.float(lonlatLS.split(',')[0])
+        lonlatL[1] = np.float(lonlatLS.split(',')[1])
+        lonlatL[2] = np.float(lonlatLS.split(',')[2])
+        lonlatL[3] = np.float(lonlatLS.split(',')[3])
+
+    grtit = values.split(':')[1].replace('!', ' ')
+    kindfig = values.split(':')[2]
+    mapvalues = values.split(':')[3]
+    ColorLeg = values.split(':')[4]
+    mervals = values.split(':')[5]
+    parvals = values.split(':')[6]
+    coastlvals = values.split(':')[7]
+    countrylvals = values.split(':')[8]
+    stslvals = values.split(':')[9]
+    close = gen.Str_Bool(values.split(':')[10])
+
+    objdomf = NetCDFFile(ncfile, 'r')
+
+    objhgt = objdomf.variables['HGT_M']
+    objlandsea = objdomf.variables['LANDMASK']
+    objlon0 = objdomf.variables['XLONG_M']
+    objlat0 = objdomf.variables['XLAT_M']
+    objlu = objdomf.variables['LU_INDEX']
+
+    lon = objlon0[0,:,:]
+    lat = objlat0[0,:,:]
+    topography = objhgt[0,:,:]
+    landsea = objlandsea[0,:,:]
+    lu = objlu[0,:,:]
+
+    lucolorsS = ColorLeg.split('|')[0]
+    ixleg = np.float(ColorLeg.split('|')[1])
+    if lucolorsS == 'usgs':
+        lucolors = drw.usgscolors
+    elif lucolorsS == 'modifmodis':
+        lucolors = drw.modigbpmodisnoahcolors
+    else:
+        print errormsg
+        print '  ' + fname + ": land-use colors dictinoary '" +lucolorsS +           \
+          "' not available !!"
+        print '    available ones:', lucolorsavail
+        quit(-1)
+
+    if mervals != 'auto':
+        Mervals = [int(mervals.split('|')[0]), str(mervals.split('|')[1]),           \
+          np.float(mervals.split('|')[2]), np.float(mervals.split('|')[3])]
+    else:
+        Mervals = [8,'#AAAAAA', 0, 0.5]
+
+    if parvals != 'auto':
+        Parvals = [int(parvals.split('|')[0]), str(parvals.split('|')[1]),           \
+          np.float(parvals.split('|')[2]), np.float(parvals.split('|')[3])]
+    else:
+        Parvals = [8,'#AAAAAA', 0, 0.5]
+
+    if coastlvals != 'auto':
+        Coastlvals = [np.float(coastlvals.split('|')[0]),                            \
+          str(coastlvals.split('|')[1])]
+    elif coastlvals == 'None':
+        Coastlvals = None
+    else:
+        Coastlvals = [0.25,'#161616']
+
+    if countrylvals != 'auto':
+        Countrylvals = [np.float(countrylvals.split('|')[0]),                        \
+          str(countrylvals.split('|')[1])]
+    elif countrylvals == 'None':
+        Countrylvals = None
+    else:
+        Countrylvals = [0.25,'#161616']
+
+    if stslvals != 'auto':
+        Stslvals = [np.float(stslvals.split('|')[0]), str(stslvals.split('|')[1])]
+    elif stslvals == 'None':
+        Stslvals = None
+    else:
+        Stslvals = [0.25,'#080808']
+
+    drw.plot_landuse(lu, lon, lat, lonlatL ,mapvalues, ixleg, lucolors, grtit,       \
+      kindfig, fname, Mervals, Parvals, Coastlvals, Countrylvals, Stslvals, close)
+
+    return
+
 #quit()
 
@@ -11590 +11734 @@
     elif oper == 'draw_ensembles_time':
         draw_ensembles_time(opts.ncfile, opts.values)
+    elif oper == 'draw_geogrid_landuse':
+        draw_geogrid_landuse(opts.ncfile, opts.values)
     elif oper == 'draw_Neighbourghood_evol':
         draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -219 +219 @@
 # plot_2Dshad_map: Function to plot a map of mask values on top of a map in shadow
 # plot_2D_shadow_time: Plotting a 2D field with one of the axes being time
+# plot_landuse: Function to plot a landuse map
 # plot_Neighbourghood_evol:Plotting neighbourghood evolution
 # plot_time_lag: Function to plot a time-lag figure (x, previous values; y, future values)
@@ -14950 +14951 @@
 
     return
+
+def plot_landuse(luv, olon, olat, lonlatLims, mapv, ixleg, landcols,figtitle, kfig,  \
+  fign, mervals=[8,'#AAAAAA',0,0.5], parvals=[8,'#AAAAAA',0,0.5],                    \
+  coastlvals=[0.25,'#AAAAAA'], countrylvals=[0.25,'#161616'], stslvals=None,         \
+  close=True):
+    """ Function to plot a landuse map
+      luv: land-use data
+      olon,olat: objects with longitudes and latitudes
+      lonlatLims: limits of the plot
+      mapv: map values
+      mervals = values for the meridians [fontsize, color line, labels rotation, 
+        line width] ([8,'#AAAAAA', 0, 0.5], default)
+      parvals = values for the parallels [fontsize, color line, labels rotation, 
+        line width] ([8,'#AAAAAA', 0, 0.5], default)
+      coastlvals = values for the coastline [line width, color line] None for any
+        ([0.25,'#161616'], default)
+      countrylvals = values for the countries [line width, color line] None for any
+        ([0.25, '#161616'], default)
+      ststlvals = values for the states [line width, color line] None for any (None,
+        default) 
+      ixleg: initial position of the legend of the colors
+      landcols: dictionary with the land colors
+      figtitle: title of the figure
+      kfig: kind of output graphical file
+      fign: name of the graphical output file
+      close: whether figure has to be close or not
+    """
+    fname = 'plot_landuse'
+
+    # https://matplotlib.org/tutorials/colors/colorbar_only.html
+    usgsc = []
+    Ncolors = len(landcols.keys())
+    labs = []
+    for iv in range(Ncolors):
+        vvals = landcols[iv+1]
+        labs.append(vvals[0])
+        usgsc.append(vvals[2])
+
+    bounds = range(Ncolors)
+    cmap = mpl.colors.ListedColormap(usgsc)
+    norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
+
+    if not mapv is None:
+        if len(olon[:].shape) == 3:
+            lon0 = olon[0,]
+            lat0 = olat[0,]
+        elif len(olon[:].shape) == 2:
+            lon0 = olon[:]
+            lat0 = olat[:]
+        elif len(olon[:].shape) == 1:
+            lon00 = olon[:]
+            lat00 = olat[:]
+            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
+            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
+
+            for iy in range(len(lat00)):
+                lon0[iy,:] = lon00
+            for ix in range(len(lon00)):
+                lat0[:,ix] = lat00
+
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+        dx = lon0.shape[1]
+        dy = lon0.shape[0]
+
+        if lonlatLims is not None:
+            print '  ' + fname + ': cutting the domain to plot !!!!'
+            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
+            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
+            nlon =  lonlatLims[0]
+            xlon =  lonlatLims[2]
+            nlat =  lonlatLims[1]
+            xlat =  lonlatLims[3]
+
+            if map_proj == 'lcc':
+                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
+                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
+        else:
+            nlon = np.min(lon0)
+            xlon = np.max(lon0)
+            nlat = np.min(lat0)
+            xlat = np.max(lat0)
+            lon2 = lon0[dy/2,dx/2]
+            lat2 = lat0[dy/2,dx/2]
+
+        plt.xlim(nlon, xlon)
+        plt.ylim(nlat, xlat)
+        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
+          xlon, ',', xlat
+
+        m, meridians, parallels, x, y = mapFigure(olon,olat,lonlatLims,mapv,5,5)
+
+    plt.pcolormesh(x,y,luv,cmap=cmap, vmin=1, vmax=Ncolors)
+    cbar = plt.colorbar(ticks=np.arange(Ncolors)+0.5)
+    cbar.ax.set_yticklabels(labs, fontsize=6)
+
+    if coastlvals is not None: 
+        m.drawcoastlines(linewidth=coastlvals[0], color=coastlvals[1])
+    if countrylvals is not None: 
+        m.drawcountries(linewidth=countrylvals[0], color=countrylvals[1])
+    if stslvals is not None: 
+        m.drawstates(linewidth=stslvals[0], color=stslvals[1])
+
+    m.drawmeridians(meridians, labels=[True,False,False,True], fontsize=mervals[0],  \
+      color=mervals[1], rotation=mervals[2], linewidth=mervals[3])
+    m.drawparallels(parallels, labels=[True,False,False,True], fontsize=parvals[0],  \
+      color=parvals[1], rotation=parvals[2], linewidth=parvals[3])
+
+    plt.title(figtitle)
+
+    output_kind(kfig, fign, close)
+    return
+