source: trunk/UTIL/PYTHON/planetoplot_v2/ppplot.py @ 938

Last change on this file since 938 was 936, checked in by aslmd, 12 years ago

UTIL PYTHON planetoplot_v2. fixed bug with minimum close to 0. number of ticks in 1D plots can be set with self.div. better labelling (actual values). better handling of time in case this is Ls.
File size: 28.3 KB
Line 
1###############################################
2## PLANETOPLOT                               ##
3## --> PPPLOT                                ##
4###############################################
5## Author: Aymeric Spiga. 02-03/2013         ##
6###############################################
7# python built-in librairies
8import time as timelib
9# added librairies
10import numpy as np
11import matplotlib.pyplot as mpl
12from matplotlib.cm import get_cmap
13from mpl_toolkits.basemap import Basemap
14from matplotlib.ticker import FormatStrFormatter,MaxNLocator
15# personal librairies
16import ppcompute
17###############################################
18
19#################################
20# global variables and settings #
21#################################
22
23# matplotlib settings
24# http://matplotlib.org/users/customizing.html
25# -------------------------------
26mpl.rcParams['font.family'] = "serif"
27mpl.rcParams['axes.color_cycle'] = "r,g,b,k"
28mpl.rcParams['contour.negative_linestyle'] = "dashed" # ou "solid"
29mpl.rcParams['verbose.level'] = "silent"
30mpl.rcParams['lines.linewidth'] = 1.5
31mpl.rcParams['lines.markersize'] = 10
32mpl.rcParams['xtick.major.pad'] = 10
33mpl.rcParams['ytick.major.pad'] = 10
34
35# global variables
36# -------------------------------
37# - where settings files are located
38#   (None means planetoplot_v2 in PYTHONPATH)
39whereset = None
40# - some good default settings.
41# (bounds)
42how_many_sigma = 3.0 
43# (contours)
44ccol = 'black'
45cline = 0.55
46# (vectors)
47widthvec = 0.002
48reducevec = 30.
49# (colorbar)
50zeorientation="vertical"
51zefrac = 0.05
52# (save figures)
53pad_inches_value=0.25
54# (negative mode)
55def_negative = False
56###############################################
57
58### settings for 'negative-like' mode
59if def_negative:
60    mpl.rc('figure', facecolor='k', edgecolor='k')
61    mpl.rcParams['text.color'] = 'w'
62    mpl.rc('axes',labelcolor='w',facecolor='k',edgecolor='w')
63    mpl.rcParams['xtick.color'] = 'w'
64    mpl.rcParams['ytick.color'] = 'w'
65    mpl.rcParams['grid.color'] = 'w'
66    mpl.rc('savefig',facecolor='k',edgecolor='k')
67
68##########################
69# executed when imported #
70##########################
71###########################################
72# we load user-defined automatic settings #
73###########################################
74# initialize the warning variable about file not present...
75files_not_present = ""
76# ... and the whereset variable
77whereset = ppcompute.findset(whereset)
78
79# - variable settings
80# -------------------------------
81zefile = "set_var.txt"
82vf = {} ; vc = {} ; vl = {} ; vu = {}
83try: 
84    f = open(whereset+zefile, 'r')
85    for line in f:
86        if "#" in line: pass
87        else:
88            var, format, colorb, label, units = line.strip().split(';')
89            ind = var.strip() 
90            vf[ind] = format.strip()
91            vc[ind] = colorb.strip()
92            vl[ind] = label.strip()
93            vu[ind] = units.strip()
94    f.close()
95except IOError: 
96    files_not_present = files_not_present + zefile + " "
97
98## - file settings
99## -------------------------------
100#zefile = "set_file.txt"
101#prefix_t = {} ; suffix_t = {} ; name_t = {} ; proj_t = {} ; vecx_t = {} ; vecy_t = {}
102#try:
103#    f = open(whereset+zefile, 'r')
104#    for line in f:
105#        if "#" in line: pass
106#        else:
107#            prefix, suffix, name, proj, vecx, vecy = line.strip().split(';')
108#            ind = name.strip()
109#            prefix_t[ind] = prefix.strip()
110#            suffix_t[ind] = suffix.strip()
111#            #name_t[ind] = name.strip()
112#            proj_t[ind] = proj.strip()
113#            vecx_t[ind] = vecx.strip()
114#            vecy_t[ind] = vecy.strip()
115#    f.close()
116#except IOError:
117#    files_not_present = files_not_present + zefile + " "
118
119# - multiplot settings
120# -------------------------------
121zefile = "set_multiplot.txt"
122subv_t = {} ; subh_t = {} ; wspace_t = {} ; hspace_t = {} ; font_t = {}
123try:
124    f = open(whereset+zefile, 'r')
125    for line in f:
126        if "#" in line: pass
127        else:
128            num, subv, subh, wspace, hspace, font = line.strip().split(';')
129            ind = int(num.strip())
130            subv_t[ind] = int(subv.strip())
131            subh_t[ind] = int(subh.strip())
132            wspace_t[ind] = float(wspace.strip())
133            hspace_t[ind] = float(hspace.strip())
134            font_t[ind] = float(font.strip())
135    f.close()
136except IOError:
137    files_not_present = files_not_present + zefile + " "
138
139# - background settings
140# -------------------------------
141zefile = "set_back.txt"
142back = {}
143try:
144    f = open(whereset+zefile, 'r')
145    for line in f:
146        if "#" in line: pass
147        else:
148            name, link = line.strip().split(';')
149            ind = name.strip() 
150            back[ind] = link.strip()
151    f.close()
152except IOError:
153    files_not_present = files_not_present + zefile + " "
154
155# - area settings
156# -------------------------------
157zefile = "set_area.txt"
158area = {}
159try:
160    f = open(whereset+zefile, 'r')
161    for line in f:
162        if "#" in line: pass
163        else:
164            name, wlat1, wlat2, wlon1, wlon2 = line.strip().split(';')
165            area[name.strip()] = [[float(wlon1.strip()),float(wlon2.strip())],\
166                                  [float(wlat1.strip()),float(wlat2.strip())]]
167    f.close()
168except IOError:
169    files_not_present = files_not_present + zefile + " "
170# A note to the user about missing files
171if files_not_present != "":
172    print "warning: files "+files_not_present+" not in "+whereset+" ; those presets will be missing"
173
174# TBD: should change vector color with colormaps
175#"gist_heat":    "white",\
176#"hot":          "white",\
177#"gray":         "red",\
178
179####################
180# useful functions #
181####################
182
183# a function to define subplot
184# ... user can change settings in set_multiplot.txt read above
185# -------------------------------
186def definesubplot(numplot, fig):
187    try: 
188        mpl.rcParams['font.size'] = font_t[numplot]
189    except: 
190        mpl.rcParams['font.size'] = 18
191    try: 
192        fig.subplots_adjust(wspace = wspace_t[numplot], hspace = hspace_t[numplot])
193        subv, subh = subv_t[numplot],subh_t[numplot]
194    except: 
195        print "!! WARNING !! no preset found from set_multiplot.txt, or this setting file was not found."
196        subv = 1 ; subh = numplot
197    return subv,subh
198
199# a function to calculate automatically bounds (or simply prescribe those)
200# -------------------------------
201def calculate_bounds(field,vmin=None,vmax=None):
202    # prescribed cases first
203    zevmin = vmin
204    zevmax = vmax
205    # computed cases
206    if zevmin is None or zevmax is None:
207       # select values
208       ind = np.where(field < 9e+35)
209       fieldcalc = field[ ind ] # field must be a numpy array
210       # calculate stdev and mean
211       dev = np.std(fieldcalc)*how_many_sigma
212       damean = ppcompute.mean(fieldcalc)
213       # fill min/max if needed
214       if vmin is None: zevmin = damean - dev
215       if vmax is None: zevmax = damean + dev
216       # special case: negative values with stddev while field is positive
217       if zevmin < 0. and ppcompute.min(fieldcalc) >= 0.: zevmin = 0.
218    # check that bounds are not too tight given the field
219    amin = ppcompute.min(field)
220    amax = ppcompute.max(field)
221    if np.abs(amin) < 1.e-15: 
222        zevmin = 0.
223        cmin = 0.
224    else:
225        cmin = 100.*np.abs((amin - zevmin)/amin)
226    cmax = 100.*np.abs((amax - zevmax)/amax)
227    if cmin > 150. or cmax > 150.:
228        print "!! WARNING !! Bounds are a bit too tight. Might need to reconsider those."
229        print "!! WARNING !! --> actual",amin,amax,"adopted",zevmin,zevmax
230    return zevmin, zevmax   
231    #### treat vmin = vmax for continuity
232    #if vmin == vmax:  zevmin = damean - dev ; zevmax = damean + dev
233
234# a function to solve the problem with blank bounds !
235# -------------------------------
236def bounds(what_I_plot,zevmin,zevmax,miss=9e+35):
237    small_enough = 1.e-7
238    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1.-small_enough) ] = zevmin*(1.-small_enough)
239    else:          what_I_plot[ what_I_plot < zevmin*(1.+small_enough) ] = zevmin*(1.+small_enough)
240    what_I_plot[ what_I_plot > miss  ] = -miss
241    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1.-small_enough)
242    return what_I_plot
243
244# a generic function to show (GUI) or save a plot (PNG,EPS,PDF,...)
245# -------------------------------
246def save(mode="gui",filename="plot",folder="./",includedate=True,res=150,custom=False):
247    if mode != "nothing":
248      # a few settings
249      possiblesave = ['eps','ps','svg','png','jpg','pdf'] 
250      # now the main instructions
251      if mode == "gui": 
252          mpl.show()
253      elif mode in possiblesave:
254          ## name of plot
255          name = folder+'/'+filename
256          if includedate:
257              for ttt in timelib.gmtime():
258                  name = name + "_" + str(ttt)
259          name = name +"."+mode
260          ## save file
261          print "**** Saving file in "+mode+" format... Please wait."
262          if not custom:
263              # ... regular plots
264              mpl.savefig(name,dpi=res,pad_inches=pad_inches_value,bbox_inches='tight')
265          else:
266              # ... mapping mode, adapted space for labels etc...
267              mpl.savefig(name,dpi=res)
268      else:
269          print "!! ERROR !! File format not supported. Supported: ",possiblesave
270
271##################################
272# a generic class to make a plot #
273##################################
274class plot():
275
276    # print out a help string when help is invoked on the object
277    # -------------------------------
278    def __repr__(self):
279        whatprint = 'plot object. \"help(plot)\" for more information\n'
280        return whatprint
281
282    # default settings
283    # -- user can define settings by two methods.
284    # -- 1. yeah = plot2d(title="foo")
285    # -- 2. yeah = pp() ; yeah.title = "foo"
286    # -------------------------------
287    def __init__(self,\
288                 var=None,\
289                 field=None,\
290                 absc=None,\
291                 xlabel="",\
292                 ylabel="",\
293                 div=20,\
294                 logx=False,\
295                 logy=False,\
296                 swap=False,\
297                 swaplab=True,\
298                 invert=False,\
299                 xcoeff=1.,\
300                 ycoeff=1.,\
301                 fmt=None,\
302                 units="",\
303                 title=""):
304        ## what could be defined by the user
305        self.var = var
306        self.field = field
307        self.absc = absc
308        self.xlabel = xlabel
309        self.ylabel = ylabel
310        self.title = title
311        self.div = div
312        self.logx = logx
313        self.logy = logy
314        self.swap = swap
315        self.swaplab = swaplab # NB: swaplab only used if swap=True
316        self.invert = invert
317        self.xcoeff = xcoeff
318        self.ycoeff = ycoeff
319        self.fmt = fmt
320        self.units = units
321        ## other useful arguments
322        ## ... not used here in ppplot but need to be attached to plot object
323        self.axisbg = "white"
324        self.superpose = False
325
326    # check
327    # -------------------------------
328    def check(self):
329        if self.field is None: print "!! ERROR !! Please define a field to be plotted" ; exit()
330
331    # define_from_var
332    # ... this uses settings in set_var.txt
333    # -------------------------------
334    def define_from_var(self):
335        if self.var is not None:
336         if self.var.upper() in vl.keys():
337          self.title = vl[self.var.upper()]
338          self.units = vu[self.var.upper()]
339
340    # make
341    # this is generic to all plots
342    # -------------------------------
343    def make(self):
344        self.check()
345        # labels, title, etc...
346        mpl.xlabel(self.xlabel)
347        mpl.ylabel(self.ylabel)
348        if self.swap:
349         if self.swaplab:
350           mpl.xlabel(self.ylabel)
351           mpl.ylabel(self.xlabel)
352        mpl.title(self.title)
353        # if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
354        if type(self.field).__name__ in 'MaskedArray':
355            self.field[self.field.mask] = self.field.fill_value
356
357################################
358# a subclass to make a 1D plot #
359################################
360class plot1d(plot):
361
362    # print out a help string when help is invoked on the object
363    # -------------------------------
364    def __repr__(self):
365        whatprint = 'plot1d object. \"help(plot1d)\" for more information\n'
366        return whatprint
367
368    # default settings
369    # -- user can define settings by two methods.
370    # -- 1. yeah = plot1d(title="foo")
371    # -- 2. yeah = pp() ; yeah.title = "foo"
372    # -------------------------------
373    def __init__(self,\
374                 lstyle='-',\
375                 color='b',\
376                 marker='x',\
377                 label=None):
378        ## get initialization from parent class
379        plot.__init__(self)
380        ## what could be defined by the user
381        self.lstyle = lstyle
382        self.color = color
383        self.marker = marker
384        self.label = label
385
386    # define_from_var
387    # ... this uses settings in set_var.txt
388    # -------------------------------
389    def define_from_var(self):
390        # get what is done in the parent class
391        plot.define_from_var(self)
392        # add specific stuff
393        if self.var is not None:
394         if self.var.upper() in vl.keys():
395          self.ylabel = vl[self.var.upper()]
396          self.title = ""
397          self.fmt = vf[self.var.upper()]
398
399    # make
400    # -------------------------------
401    def make(self):
402        # get what is done in the parent class
403        plot.make(self)
404        if self.fmt is None: self.fmt = '%.0f'
405        # add specific stuff
406        mpl.grid(color='grey')
407        if self.lstyle == "": self.lstyle = " " # to allow for no line at all with ""
408        # swapping if requested
409        if self.swap:  x = self.field ; y = self.absc
410        else:          x = self.absc ; y = self.field
411        # coefficients on axis
412        x=x*self.xcoeff ; y=y*self.ycoeff
413        # check axis
414        if x.size != y.size:
415            print "!! ERROR !! x and y sizes don't match on 1D plot.", x.size, y.size
416            exit()
417        # make the 1D plot
418        # either request linestyle or let matplotlib decide
419        if self.lstyle is not None and self.color is not None:
420            mpl.plot(x,y,self.color+self.lstyle,marker=self.marker,label=self.label)
421        else:
422            mpl.plot(x,y,marker=self.marker,label=self.label)
423        # make log axes and/or invert ordinate
424        # ... this must be after plot so that axis bounds are well-defined
425        # ... also inverting must be after making the thing logarithmic
426        if self.logx: mpl.semilogx()
427        if self.logy: mpl.semilogy()
428        if self.invert: ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
429        # add a label for line(s)
430        if self.label is not None:
431            if self.label != "":
432                mpl.legend(loc="best",fancybox=True)
433        # AXES
434        ax = mpl.gca()
435        # format labels
436        if self.swap: ax.xaxis.set_major_formatter(FormatStrFormatter(self.fmt))
437        else: ax.yaxis.set_major_formatter(FormatStrFormatter(self.fmt))
438        # plot limits: ensure that no curve would be outside the window
439        # x-axis
440        x1, x2 = ax.get_xbound()
441        xmin = ppcompute.min(x)
442        xmax = ppcompute.max(x)
443        if xmin < x1: x1 = xmin
444        if xmax > x2: x2 = xmax
445        ax.set_xbound(lower=x1,upper=x2)
446        # y-axis
447        y1, y2 = ax.get_ybound()
448        ymin = ppcompute.min(y)
449        ymax = ppcompute.max(y)
450        if ymin < y1: y1 = ymin
451        if ymax > y2: y2 = ymax
452        ax.set_ybound(lower=y1,upper=y2)
453        ## set with .div the number of ticks. (is it better than automatic?)
454        ax.xaxis.set_major_locator(MaxNLocator(self.div/2))
455
456################################
457# a subclass to make a 2D plot #
458################################
459class plot2d(plot):
460
461    # print out a help string when help is invoked on the object
462    # -------------------------------
463    def __repr__(self):
464        whatprint = 'plot2d object. \"help(plot2d)\" for more information\n'
465        return whatprint
466
467    # default settings
468    # -- user can define settings by two methods.
469    # -- 1. yeah = plot2d(title="foo")
470    # -- 2. yeah = pp() ; yeah.title = "foo"
471    # -------------------------------
472    def __init__(self,\
473                 ordi=None,\
474                 mapmode=False,\
475                 proj="cyl",\
476                 back=None,\
477                 colorb="jet",\
478                 trans=1.0,\
479                 addvecx=None,\
480                 addvecy=None,\
481                 addcontour=None,\
482                 blon=None,\
483                 blat=None,\
484                 area=None,\
485                 vmin=None,\
486                 vmax=None,\
487                 colorvec="black"):
488        ## get initialization from parent class
489        plot.__init__(self)
490        ## what could be defined by the user
491        self.ordi = ordi
492        self.mapmode = mapmode
493        self.proj = proj
494        self.back = back
495        self.colorb = colorb
496        self.trans = trans
497        self.addvecx = addvecx
498        self.addvecy = addvecy
499        self.colorvec = colorvec
500        self.addcontour = addcontour
501        self.blon = blon ; self.blat = blat
502        self.area = area
503        self.vmin = vmin ; self.vmax = vmax
504
505    # define_from_var
506    # ... this uses settings in set_var.txt
507    # -------------------------------
508    def define_from_var(self):
509        # get what is done in the parent class
510        plot.define_from_var(self)
511        # add specific stuff
512        if self.var is not None:
513         if self.var.upper() in vl.keys():
514          self.colorb = vc[self.var.upper()]
515          self.fmt = vf[self.var.upper()]
516
517    # make
518    # -------------------------------
519    def make(self):
520        # get what is done in the parent class...
521        plot.make(self)
522        if self.fmt is None: self.fmt = "%.2e"
523        # ... then add specific stuff
524        ############################################################################################
525        ### PRE-SETTINGS
526        ############################################################################################
527        # transposing if necessary
528        shape = self.field.shape
529        if shape[0] != shape[1]:
530         if len(self.absc) == shape[0] and len(self.ordi) == shape[1]:
531            print "!! WARNING !! Transposing axes"
532            self.field = np.transpose(self.field)
533        # bound field
534        zevmin, zevmax = calculate_bounds(self.field,vmin=self.vmin,vmax=self.vmax)
535        what_I_plot = bounds(self.field,zevmin,zevmax)
536        # define contour field levels. define color palette
537        ticks = self.div + 1
538        zelevels = np.linspace(zevmin,zevmax,ticks)
539        palette = get_cmap(name=self.colorb)
540        # do the same thing for possible contourline entries
541        if self.addcontour is not None:
542            # if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
543            if type(self.addcontour).__name__ in 'MaskedArray':
544               self.addcontour[self.addcontour.mask] = self.addcontour.fill_value
545            zevminc, zevmaxc = calculate_bounds(self.addcontour)
546            what_I_contour = bounds(self.addcontour,zevminc,zevmaxc)
547            ticks = self.div + 1
548            zelevelsc = np.linspace(zevminc,zevmaxc,ticks)
549        ############################################################################################
550        ### MAIN PLOT
551        ### NB: contour lines are done before contour shades otherwise colorar error
552        ############################################################################################
553        if not self.mapmode:
554            ## A SIMPLE 2D PLOT
555            ###################
556            # swapping if requested
557            if self.swap:  x = self.ordi ; y = self.absc
558            else:          x = self.absc ; y = self.ordi
559            # coefficients on axis
560            x=x*self.xcoeff ; y=y*self.ycoeff
561            # make shaded and line contours
562            if self.addcontour is not None: 
563                mpl.contour(x, y, what_I_contour, \
564                            zelevelsc, colors = ccol, linewidths = cline)
565            mpl.contourf(x, y, \
566                         self.field, \
567                         zelevels, cmap=palette)
568            # make log axes and/or invert ordinate
569            if self.logx: mpl.semilogx()
570            if self.logy: mpl.semilogy()
571            if self.invert: ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
572        else:
573            ## A 2D MAP USING PROJECTIONS (basemap)
574            #######################################
575            mpl.xlabel("") ; mpl.ylabel("")
576            # additional security in case self.proj is None here
577            # ... we set cylindrical projection (the simplest one)
578            if self.proj is None: self.proj = "cyl"
579            # get lon and lat in 2D version.
580            # (but first ensure we do have 2D coordinates)
581            if self.absc.ndim == 1:     [self.absc,self.ordi] = np.meshgrid(self.absc,self.ordi)
582            elif self.absc.ndim > 2:    print "!! ERROR !! lon and lat arrays must be 1D or 2D"
583            # get lat lon intervals and associated settings
584            wlon = [np.min(self.absc),np.max(self.absc)]
585            wlat = [np.min(self.ordi),np.max(self.ordi)]
586            # -- area presets are in set_area.txt
587            if self.area is not None:
588             if self.area in area.keys():
589                wlon, wlat = area[self.area]
590            # -- settings for meridians and parallels
591            steplon = int(abs(wlon[1]-wlon[0])/6.)
592            steplat = int(abs(wlat[1]-wlat[0])/3.)
593            #mertab = np.r_[wlon[0]:wlon[1]:steplon] ; merlab = [0,0,0,1]
594            #partab = np.r_[wlat[0]:wlat[1]:steplat] ; parlab = [1,0,0,0]
595            mertab = np.r_[-360.:360.:steplon] ; merlab = [0,0,0,1]
596            partab = np.r_[-90.:90.:steplat] ; parlab = [1,0,0,0]
597            format = '%.1f'
598            # -- center of domain and bounding lats
599            lon_0 = 0.5*(wlon[0]+wlon[1])
600            lat_0 = 0.5*(wlat[0]+wlat[1])
601            # some tests, bug fixes, and good-looking settings
602            # ... cyl is good for global and regional
603            if self.proj == "cyl":
604                format = '%.0f'
605            # ... global projections
606            elif self.proj in ["ortho","moll","robin"]:
607                wlat[0] = None ; wlat[1] = None ; wlon[0] = None ; wlon[1] = None
608                steplon = 30. ; steplat = 30.
609                if self.proj in ["robin","moll"]: steplon = 60.
610                mertab = np.r_[-360.:360.:steplon]
611                partab = np.r_[-90.:90.:steplat]
612                if self.proj == "ortho": 
613                    merlab = [0,0,0,0] ; parlab = [0,0,0,0]
614                    # in ortho projection, blon and blat can be used to set map center
615                    if self.blon is not None: lon_0 = self.blon
616                    if self.blat is not None: lat_0 = self.blat
617                elif self.proj == "moll":
618                    merlab = [0,0,0,0]
619                format = '%.0f'
620            # ... regional projections
621            elif self.proj in ["lcc","laea","merc"]:
622                if self.proj in ["lcc","laea"] and wlat[0] == -wlat[1]: 
623                    print "!! ERROR !! with Lambert lat1 must be different than lat2" ; exit()
624                if wlat[0] < -80. and wlat[1] > 80.:
625                    print "!! ERROR !! set an area (not global)" ; exit()
626                format = '%.0f'
627            elif self.proj in ["npstere","spstere"]:
628                # in polar projections, blat gives the bounding lat
629                # if not set, set something reasonable
630                if self.blat is None:   self.blat = 60.
631                # help the user who forgets self.blat would better be negative in spstere
632                # (this actually serves for the default setting just above)
633                if self.proj == "spstere" and self.blat > 0: self.blat = -self.blat
634                # labels
635                mertab = np.r_[-360.:360.:15.]
636                partab = np.r_[-90.:90.:5.]
637            # ... unsupported projections
638            else:
639                print "!! ERROR !! unsupported projection. supported: "+\
640                      "cyl, npstere, spstere, ortho, moll, robin, lcc, laea, merc"
641            # finally define projection
642            m = Basemap(projection=self.proj,\
643                        lat_0=lat_0,lon_0=lon_0,\
644                        boundinglat=self.blat,\
645                        llcrnrlat=wlat[0],urcrnrlat=wlat[1],\
646                        llcrnrlon=wlon[0],urcrnrlon=wlon[1])
647            # draw meridians and parallels
648            ft = int(mpl.rcParams['font.size']*3./4.)
649            zelatmax = 85.
650            m.drawmeridians(mertab,labels=merlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
651            m.drawparallels(partab,labels=parlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
652            # define background (see set_back.txt)
653            if self.back is not None:
654              if self.back in back.keys():
655                 print "**** info: loading a background, please wait.",self.back
656                 if self.back not in ["coast","sea"]:   m.warpimage(back[self.back],scale=0.75)
657                 elif self.back == "coast":             m.drawcoastlines()
658                 elif self.back == "sea":               m.drawlsmask(land_color='white',ocean_color='aqua')
659              else:
660                 print "!! ERROR !! requested background not defined. change name or fill in set_back.txt" ; exit()
661            # define x and y given the projection
662            x, y = m(self.absc, self.ordi)
663            # contour field. first line contour then shaded contour.
664            if self.addcontour is not None: 
665                m.contour(x, y, what_I_contour, \
666                            zelevelsc, colors = ccol, linewidths = cline)
667            m.contourf(x, y, what_I_plot, zelevels, cmap = palette, alpha = self.trans)
668        ############################################################################################
669        ### COLORBAR
670        ############################################################################################
671        if self.trans > 0.:
672            ## draw colorbar. settings are different with projections. or if not mapmode.
673            if not self.mapmode: orientation=zeorientation ; frac = 0.075 ; pad = 0.03 ; lu = 0.5
674            elif self.proj in ['moll']: orientation="horizontal" ; frac = 0.075 ; pad = 0.03 ; lu = 1.0
675            elif self.proj in ['cyl']: orientation="vertical" ; frac = 0.023 ; pad = 0.03 ; lu = 0.5
676            else: orientation = zeorientation ; frac = zefrac ; pad = 0.03 ; lu = 0.5
677            zelevpal = np.linspace(zevmin,zevmax,num=min([ticks/2+1,21]))
678            zecb = mpl.colorbar(fraction=frac,pad=pad,\
679                                format=self.fmt,orientation=orientation,\
680                                ticks=zelevpal,\
681                                extend='neither',spacing='proportional')
682            if zeorientation == "horizontal": zecb.ax.set_xlabel(self.title) ; self.title = ""
683            # colorbar title --> units
684            if self.units not in ["dimless",""]:
685                zecb.ax.set_title("["+self.units+"]",fontsize=3.*mpl.rcParams['font.size']/4.,x=lu,y=1.025)
686
687        ############################################################################################
688        ### VECTORS. must be after the colorbar. we could also leave possibility for streamlines.
689        ############################################################################################
690        ### not expecting NaN in self.addvecx and self.addvecy. masked arrays is just enough.
691        stride = 3
692        if self.addvecx is not None and self.addvecy is not None and self.mapmode:
693                ### for metwinds only ???
694                [vecx,vecy] = m.rotate_vector(self.addvecx,self.addvecy,self.absc,self.ordi) 
695                #vecx,vecy = self.addvecx,self.addvecy
696                # reference vector is scaled
697                zescale = ppcompute.mean(np.sqrt(self.addvecx*self.addvecx+self.addvecy*self.addvecy))
698                # make vector field
699                q = m.quiver( x[::stride,::stride],y[::stride,::stride],\
700                              vecx[::stride,::stride],vecy[::stride,::stride],\
701                              angles='xy',color=self.colorvec,pivot='middle',\
702                              scale=zescale*reducevec,width=widthvec )
703                # make vector key.
704                #keyh = 1.025 ; keyv = 1.05 # upper right corner over colorbar
705                keyh = 0.98 ; keyv = 1.06
706                #keyh = -0.03 ; keyv = 1.08 # upper left corner
707                p = mpl.quiverkey(q,keyh,keyv,\
708                                  zescale,str(int(zescale)),\
709                                  color='black',labelpos='S',labelsep = 0.07)
Note: See TracBrowser for help on using the repository browser.