Changeset 2014 in lmdz_wrf

Timestamp:

Aug 2, 2018, 5:25:56 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `draw_stations_map': Function to plot a map with the stations

Location:

trunk/tools

Files:

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

trunk/tools/drawing.py

@@ -64 +64 @@
 ## e.g. # drawing.py -o draw_2D_shad_contdisc -f 'simout_sfcdiags.nc;tdas;XLONG;XLAT;Time|0,time|0,west_east|-1,south_north|-1@all_single-stations.nc;tdas;stslon;stslat;time|0,lon|-1,lat|-1' -S 'tdas:west_east,south_north:auto:YlGnBu,auto,auto:286.,298.:auto:obs!&!sim!tdas!on!1995/01/01!00!UTC:png:None:cyl,f:sponge,0.1,0.1:yes'
 ## e.g. # drawing.py -o draw_multiWindRose -f '/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/sounding_uava_10868.nc@pres|-1;time|1@ua,va#/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/simout_vars_sndpt_10868_38lev.nc@bottom_top|-1;time|0@ua,va#/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/simout_vars_sndpt_10868_50lev.nc@bottom_top|-1;time|0@ua,va#/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/simout_vars_sndpt_10868_50lev_assigned.nc@bottom_top|-1;time|0@ua,va#/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/simout_vars_sndpt_10868_80lev.nc@bottom_top|-1;time|0@ua,va#/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/simout_vars_sndpt_10868_120lev.nc@bottom_top|-1;time|0@ua,va#/media/lluis/ExtDiskC_ext3/DATA/estudios/FPS_Alps/additional/IOP/select/out/simout_vars_sndpt_10868_NOaerosol.nc@bottom_top|-1;time|0@ua,va' -S 'linepoint;multicol;pres;-;auto;spectral;auto:obs,38lev,50lev,50leva,80lev,120lev:WindRose!obs!sim!comparison!on!2012/10/23!00!UTC:png:auto:True'
+## e.g. # drawing.py -o draw_stations_map -f stations.inf -S 'snd|x|#FFAAAA|2|8@sfc|x|#AAFFAA|2|8:cyl,l:1.,20.,35.,49.:FPS!Alps!stations:pdf:true'
 
 #######
@@ -98 +99 @@
 # draw_ptZvals: Function to plot a given list of points by their Z value and a colorbar
 # draw_SkewT: creation of a SkewT-logP diagram using matplotlib's API example  
+# draw_stations_map: Function to plot a map with the stations
 # draw_time_lag: Function to plot a time-lag figure with multiple sources (x, previous values; y, future values)
 # draw_timeSeries: Function to draw a time-series
@@ -131 +133 @@
   'draw_Neighbourghood_evol',                                                        \
   'draw_points', 'draw_points_lonlat',                                               \
-  'draw_ptZvals', 'draw_river_desc', 'draw_SkewT', 'draw_subbasin', 'draw_Taylor',   \
+  'draw_ptZvals', 'draw_river_desc', 'draw_SkewT', 'draw_subbasin',                  \
+  'draw_stations_map', 'draw_Taylor',                                                \
   'draw_time_lag', 'draw_timeSeries', 'draw_topo_geogrid',                           \
   'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
@@ -10165 +10168 @@
 #draw_multiWindRose(filenames, figvalues)
 
+def draw_stations_map(filename, values):
+    """ Function to plot a map with the stations
+      filename= name of ASCII file with the values of the stations as: ('#' comments)
+        [sttype],[label],[lon],[lat],[height]
+          [sttype]: type of the station
+          [label]: label of the station (LaTeX format, 'None' for no label)
+          [lon]: longitude of the station
+          [lat]: latitude of the station
+          [height]: height of the station
+      values=[stypevals]:[mapvals]:[SWNEbox]:[imgtit]:[imgkind]:[close]
+        stypevals='@' separated list of [sttype]|[marker]|[color]|[size]|[fontsize]
+          [sttype]: label of the type of the marker
+          [marker]: type of the marker
+          [color]: color of the marker
+          [size]: size of the marker
+        mapvals= 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
+        SWNEbox: ',' list of the vertexs of the map [SWlon, NElon, SWlat, NElat]
+        imgtit: title of the image ('!' for spaces)
+        imgkind: kind of file output of the image (ps, pns, pdf, ...)
+        close: whether figure should be closed or not
+    """
+    fname = 'draw_stations_map'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_stations_map.__doc__
+        quit()
+
+    expectargs = '[stypevals]:[mapvals]:[SWNEbox]:[imgtit]:[imgkind]:[close]'
+    drw.check_arguments(fname,values,expectargs,':')
+
+    stypevals = values.split(':')[0].split('@')
+    mapvals = values.split(':')[1]
+    SWNEbox = gen.str_list_k(values.split(':')[2], ',', 'F')
+    imgtit = values.split(':')[3].replace('!', ' ')
+    imgkind = values.split(':')[4]
+    close = gen.Str_Bool(values.split(':')[5])
+
+    if not os.path.isfile(filename):
+        print errormsg
+        print '  ' + fname + ": file '" + filename + "' does not exist !!"
+        quit(-1)
+
+    # Getting stations values
+    styypes = {}
+    sttyps = []
+    for stvs in stypevals:
+        stv = stvs.split('|')
+        styp = stv[0]
+        smk = stv[1]
+        scol = stv[2]
+        ssz = np.float(stv[3])
+        fsz = stv[4]
+        styypes[styp] = [smk, scol, ssz, fsz]
+        sttyps.append(styp)
+
+    ofile = open(filename, 'r')
+    dicstations = {}
+    for line in ofile:
+        if line[0:1] != '#' and len(line) > 1: 
+            linev = line.replace('\n','').replace('\r','').split(',')
+            ltyp = linev[0]
+            lid = linev[1]
+            llon = np.float(linev[2])
+            llat = np.float(linev[3])
+            lhgt = np.float(linev[4])
+
+            if not gen.searchInlist(sttyps, ltyp):
+                print errormsg
+                print '  ' + fname + ": station type '" + ltyp + "' not ready !!"
+                print '    available ones:', sttyps
+                quit(-1)
+            
+            typsv = styypes[ltyp]
+            dicstations[lid] = [llon, llat, lhgt, lid, typsv[0], typsv[1], typsv[2], \
+              typsv[3]]
+
+    ofile.close()
+
+    Nd = 100
+    nlon = SWNEbox[0]
+    xlon = SWNEbox[1]
+    nlat = SWNEbox[2]
+    xlat = SWNEbox[3]
+
+    ddlon = (xlon - nlon)/Nd
+    ddlat = (xlat - nlat)/Nd
+
+    lon1D = np.arange(nlon, xlon+ddlon, ddlon)
+    lat1D = np.arange(nlat, xlat+ddlat, ddlat)
+   
+    drw.plot_stations_map(dicstations, lon1D, lat1D, mapvals, SWNEbox, ddlon, imgtit,\
+      imgkind, close)
+
+    return
+
 #quit()
 
@@ -10299 +10408 @@
     elif oper == 'draw_SkewT':
         draw_SkewT(opts.ncfile, opts.values, opts.varname)
+    elif oper == 'draw_stations_map':
+        draw_stations_map(opts.ncfile, opts.values)
     elif oper == 'draw_Taylor':
         draw_Taylor(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -187 +187 @@
 # plot_multiWindRose: Function to plot multi wind roses (from where the dinw blows)
 # plot_SkewT: Function to plot a SkewT-logP diagram using matplotlib's API example 
+# plot_stations_map: Function to plot a map with a series of stations
 # plot_Taylor: Function to draw a Taylor diagram (Taylor 2001)
 # plot_TimeEnsembles: Function to plot a of an Ensemble of values with a time-axis
@@ -12630 +12631 @@
 
     return
+
+def plot_stations_map(dicstations, dimxv, dimyv, mapv, lonlatbox, ddlon, figtitle,   \
+  kfig, close):
+    """ Function to plot a map with a series of stations
+      dicstations: dictionary with the values of the stations as:
+        dicstations[refst]= [[lon], [lat], [height], [label], [marker], [color], 
+            [size], [fontsize]]
+          [lon]: longitude of the station
+          [lat]: latitude of the station
+          [height]: height of the station
+          [label]: label in the figure ('None' for no label)
+          [marker]: type of marker
+          [color]: coor of the marker
+          [size]: size of the marker   
+          [fontsize]: size of the label font
+      dimxv: values along x-axis
+      dimyv: values along y-axis
+      mapv= 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
+      lonlatbox: list of the vertexs of the map [SWlon, NElon, SWlat, NElat]
+      ddlon: diferential of longitude
+      figtitle: title of the figure
+      kfig: kind of figure file (ps, png, pdf, ...) 
+      close: whether figure should be closed or not
+    """
+    fname = 'plot_stations_map'
+
+    nlon = lonlatbox[0]
+    xlon = lonlatbox[1]
+    nlat = lonlatbox[2]
+    xlat = lonlatbox[3]
+
+    plt.rc('text', usetex=True)
+
+    if not mapv is None:
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        lon, lat = np.meshgrid(dimxv, dimyv)
+
+        lon2 = (nlon + xlon)/2.
+        lat2 = (nlat + xlat)/2.
+
+        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
+          xlon, ',', xlat
+
+        if map_proj == 'cyl':
+            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
+              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        elif map_proj == 'lcc':
+            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
+              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        else:
+            print errormsg
+            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
+            print '    available: cyl, lcc'
+            quit(-1)
+ 
+        x,y = m(lon,lat)
+
+        m.drawcoastlines()
+
+        meridians = gen.pretty_int(nlon,xlon,5)
+        m.drawmeridians(meridians,labels=[True,False,False,True])
+
+        parallels = gen.pretty_int(nlat,xlat,5)
+        m.drawparallels(parallels,labels=[False,True,True,False])
+
+    else:
+        # No following data values
+        x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv))
+        y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv))
+
+    for ist in dicstations.keys():
+        stv = dicstations[ist]
+        stlon = stv[0]
+        stlat = stv[1]
+        stlab = stv[3]
+        stmrk = stv[4]
+        stcol = stv[5]
+        stsiz = stv[6]
+        fontsiz = stv[7]
+
+        plt.plot(stlon, stlat, stmrk, color=stcol, markersize=stsiz)
+        if stlab != 'None':
+            plt.annotate(gen.latex_text(stlab), xy=(stlon+ddlon,stlat),         \
+              xycoords='data', color=stcol, fontsize=fontsiz)
+
+    plt.title(gen.latex_text(figtitle))
+
+    figname = 'stations_map'
+    output_kind(kfig, figname, close)
+
+    return
+