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

Last change on this file since 1062 was 1062, checked in by aslmd, 11 years ago

planetoplot_v2. added trans back showcb attributes to ppclass. simplified example scripts. added the web home page example.
File size: 35.8 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'] = "b,r,g,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
46cline = 0.8
47# (vectors)
48widthvec = 0.002
49reducevec = 30.
50# (colorbar)
51zeorientation="vertical"
52zefrac = 0.05
53# (save figures)
54pad_inches_value=0.25
55# (negative mode)
56def_negative = False
57###############################################
58
59### settings for 'negative-like' mode
60if def_negative:
61    mpl.rc('figure', facecolor='k', edgecolor='k')
62    mpl.rcParams['text.color'] = 'w'
63    mpl.rc('axes',labelcolor='w',facecolor='k',edgecolor='w')
64    mpl.rcParams['xtick.color'] = 'w'
65    mpl.rcParams['ytick.color'] = 'w'
66    mpl.rcParams['grid.color'] = 'w'
67    mpl.rc('savefig',facecolor='k',edgecolor='k')
68
69##########################
70# executed when imported #
71##########################
72###########################################
73# we load user-defined automatic settings #
74###########################################
75# initialize the warning variable about file not present...
76files_not_present = ""
77# ... and the whereset variable
78whereset = ppcompute.findset(whereset)
79
80# - variable settings
81# -------------------------------
82zefile = "set_var.txt"
83vf = {} ; vc = {} ; vl = {} ; vu = {}
84try: 
85    f = open(whereset+zefile, 'r')
86    for line in f:
87        if "#" in line: pass
88        else:
89            var, format, colorb, label, units = line.strip().split(';')
90            ind = var.strip() 
91            vf[ind] = format.strip()
92            vc[ind] = colorb.strip()
93            vl[ind] = label.strip()
94            vu[ind] = units.strip()
95    f.close()
96except IOError: 
97    files_not_present = files_not_present + zefile + " "
98
99## - file settings
100## -------------------------------
101#zefile = "set_file.txt"
102#prefix_t = {} ; suffix_t = {} ; name_t = {} ; proj_t = {} ; vecx_t = {} ; vecy_t = {}
103#try:
104#    f = open(whereset+zefile, 'r')
105#    for line in f:
106#        if "#" in line: pass
107#        else:
108#            prefix, suffix, name, proj, vecx, vecy = line.strip().split(';')
109#            ind = name.strip()
110#            prefix_t[ind] = prefix.strip()
111#            suffix_t[ind] = suffix.strip()
112#            #name_t[ind] = name.strip()
113#            proj_t[ind] = proj.strip()
114#            vecx_t[ind] = vecx.strip()
115#            vecy_t[ind] = vecy.strip()
116#    f.close()
117#except IOError:
118#    files_not_present = files_not_present + zefile + " "
119
120# - multiplot settings
121# -------------------------------
122zefile = "set_multiplot.txt"
123subv_t = {} ; subh_t = {} ; wspace_t = {} ; hspace_t = {} ; font_t = {}
124try:
125    f = open(whereset+zefile, 'r')
126    for line in f:
127        if "#" in line: pass
128        else:
129            num, subv, subh, wspace, hspace, font = line.strip().split(';')
130            ind = int(num.strip())
131            subv_t[ind] = int(subv.strip())
132            subh_t[ind] = int(subh.strip())
133            wspace_t[ind] = float(wspace.strip())
134            hspace_t[ind] = float(hspace.strip())
135            font_t[ind] = float(font.strip())
136    f.close()
137except IOError:
138    files_not_present = files_not_present + zefile + " "
139
140# - background settings
141# -------------------------------
142zefile = "set_back.txt"
143back = {}
144try:
145    f = open(whereset+zefile, 'r')
146    for line in f:
147        if "#" in line: pass
148        else:
149            name, link = line.strip().split(';')
150            ind = name.strip() 
151            back[ind] = link.strip()
152    f.close()
153except IOError:
154    files_not_present = files_not_present + zefile + " "
155
156# - area settings
157# -------------------------------
158zefile = "set_area.txt"
159area = {}
160try:
161    f = open(whereset+zefile, 'r')
162    for line in f:
163        if "#" in line: pass
164        else:
165            name, wlat1, wlat2, wlon1, wlon2 = line.strip().split(';')
166            area[name.strip()] = [[float(wlon1.strip()),float(wlon2.strip())],\
167                                  [float(wlat1.strip()),float(wlat2.strip())]]
168    f.close()
169except IOError:
170    files_not_present = files_not_present + zefile + " "
171# A note to the user about missing files
172if files_not_present != "":
173    print "warning: files "+files_not_present+" not in "+whereset+" ; those presets will be missing"
174
175# TBD: should change vector color with colormaps
176#"gist_heat":    "white",\
177#"hot":          "white",\
178#"gray":         "red",\
179
180####################
181# useful functions #
182####################
183
184# continuity with matplotlib
185def close():
186    mpl.close()
187
188# a function for predefined figure sizes
189def figuref(x=16,y=9):
190    fig = mpl.figure(figsize=(x,y))
191    return fig
192
193# a function to change color lines according to a color map (idea by A. Pottier)
194# ... two optional arguments: color maps + a number telling how much intervals are needed
195def rainbow(cb="jet",num=8):
196    ax = mpl.gca()
197    pal = mpl.cm.get_cmap(name=cb)
198    ax._get_lines.set_color_cycle([pal(i) for i in np.linspace(0,0.9,num)])
199
200# a function to define subplot
201# ... user can change settings in set_multiplot.txt read above
202# -------------------------------
203def definesubplot(numplot, fig):
204    try: 
205        mpl.rcParams['font.size'] = font_t[numplot]
206    except: 
207        mpl.rcParams['font.size'] = 18
208    try: 
209        fig.subplots_adjust(wspace = wspace_t[numplot], hspace = hspace_t[numplot])
210        subv, subh = subv_t[numplot],subh_t[numplot]
211    except: 
212        print "!! WARNING !! (ignore if superpose mode) no preset found from set_multiplot.txt, or this setting file was not found."
213        subv = 1 ; subh = numplot
214    return subv,subh
215
216# a function to calculate automatically bounds (or simply prescribe those)
217# -------------------------------
218def calculate_bounds(field,vmin=None,vmax=None):
219    # prescribed cases first
220    zevmin = vmin
221    zevmax = vmax
222    # computed cases
223    if zevmin is None or zevmax is None:
224       # select values
225       ind = np.where(field < 9e+35)
226       fieldcalc = field[ ind ] # field must be a numpy array
227       # calculate stdev and mean
228       dev = np.std(fieldcalc)*how_many_sigma
229       damean = ppcompute.mean(fieldcalc)
230       # fill min/max if needed
231       if vmin is None: zevmin = damean - dev
232       if vmax is None: zevmax = damean + dev
233       # special case: negative values with stddev while field is positive
234       if zevmin < 0. and ppcompute.min(fieldcalc) >= 0.: zevmin = 0.
235    # check that bounds are not too tight given the field
236    amin = ppcompute.min(field)
237    amax = ppcompute.max(field)
238    if np.abs(amin) < 1.e-15: 
239        cmin = 0.
240    else:
241        cmin = 100.*np.abs((amin - zevmin)/amin)
242    cmax = 100.*np.abs((amax - zevmax)/amax)
243    if cmin > 150. or cmax > 150.:
244        print "!! WARNING !! Bounds are a bit too tight. Might need to reconsider those."
245        print "!! WARNING !! --> actual",amin,amax,"adopted",zevmin,zevmax
246    return zevmin, zevmax   
247    #### treat vmin = vmax for continuity
248    #if vmin == vmax:  zevmin = damean - dev ; zevmax = damean + dev
249
250# a function to solve the problem with blank bounds !
251# -------------------------------
252def bounds(what_I_plot,zevmin,zevmax,miss=9e+35):
253    small_enough = 1.e-7
254    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1.-small_enough) ] = zevmin*(1.-small_enough)
255    else:          what_I_plot[ what_I_plot < zevmin*(1.+small_enough) ] = zevmin*(1.+small_enough)
256    what_I_plot[ what_I_plot > miss  ] = -miss
257    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1.-small_enough)
258    return what_I_plot
259
260# a function to change labels with modulo
261# ---------------------------------------
262def labelmodulo(ax,mod):
263    mpl.draw()
264    strtab = []
265    for tick in ax.get_xaxis().get_ticklabels():
266        onetick = tick.get_text()
267        if len(onetick) > 0:
268          num = float(onetick)
269          strtab.append(num % mod)
270    ax.get_xaxis().set_ticklabels(strtab)
271    return ax
272
273# a function to output an ascii file
274# ----------------------------------
275def writeascii (field=None,absc=None,name=None):
276    field = np.array(field)
277    absc = np.array(absc)
278    if name is None:
279        name = "prof"
280        for ttt in timelib.gmtime():
281            name = name + "_" + str(ttt)
282    if field is None or absc is None:
283        print "!! WARNING !! Not printing the file, incorrect field or absc."
284    else:
285        if field.ndim == 1:
286            myfile = open(name, 'w')
287            for ix in range(len(absc)):
288                myfile.write("%15.5e%15.5e\n" % ( absc[ix], field[ix] ))
289            myfile.close()
290        else:
291            print "!! WARNING !! Not printing the file, 2D fields not supported yet."
292    return
293
294# a generic function to show (GUI) or save a plot (PNG,EPS,PDF,...)
295# -------------------------------
296def save(mode="gui",filename="plot",folder="./",includedate=True,res=150,custom=False):
297    if mode != "nothing":
298      # a few settings
299      possiblesave = ['eps','ps','svg','png','jpg','pdf'] 
300      # now the main instructions
301      if mode == "gui": 
302          mpl.show()
303      elif mode in possiblesave:
304          ## name of plot
305          name = folder+'/'+filename
306          if includedate:
307              for ttt in timelib.gmtime():
308                  name = name + "_" + str(ttt)
309          name = name +"."+mode
310          ## save file
311          print "**** Saving file in "+mode+" format... Please wait."
312          if not custom:
313              # ... regular plots
314              mpl.savefig(name,dpi=res,pad_inches=pad_inches_value,bbox_inches='tight')
315          else:
316              # ... mapping mode, adapted space for labels etc...
317              mpl.savefig(name,dpi=res)
318      else:
319          print "!! ERROR !! File format not supported. Supported: ",possiblesave
320
321# a necessary addition for people used to matplotlib
322def show():
323    mpl.show()
324
325##################################
326# a generic class to make a plot #
327##################################
328class plot():
329
330    # print out a help string when help is invoked on the object
331    # -------------------------------
332    def __repr__(self):
333        whatprint = 'plot object. \"help(plot)\" for more information\n'
334        return whatprint
335
336    # default settings
337    # -- user can define settings by two methods.
338    # -- 1. yeah = plot2d(title="foo")
339    # -- 2. yeah = pp() ; yeah.title = "foo"
340    # -------------------------------
341    def __init__(self,\
342                 var=None,\
343                 f=None,\
344                 x=None,\
345                 xlabel="",\
346                 ylabel="",\
347                 div=20,\
348                 logx=False,\
349                 logy=False,\
350                 swap=False,\
351                 swaplab=True,\
352                 invert=False,\
353                 xcoeff=1.,\
354                 ycoeff=1.,\
355                 fmt=None,\
356                 colorbar="jet",\
357                 units="",\
358                 modx=None,\
359                 xmin=None,\
360                 ymin=None,\
361                 xmax=None,\
362                 ymax=None,\
363                 nxticks=10,\
364                 nyticks=10,\
365                 title=""):
366        ## what could be defined by the user
367        self.var = var
368        self.f = f
369        self.x = x
370        self.xlabel = xlabel
371        self.ylabel = ylabel
372        self.title = title
373        self.div = div
374        self.logx = logx
375        self.logy = logy
376        self.swap = swap
377        self.swaplab = swaplab # NB: swaplab only used if swap=True
378        self.invert = invert
379        self.xcoeff = xcoeff
380        self.ycoeff = ycoeff
381        self.fmt = fmt
382        self.units = units
383        self.colorbar = colorbar
384        self.modx = modx
385        self.xmin = xmin
386        self.ymin = ymin
387        self.xmax = xmax
388        self.ymax = ymax
389        self.nxticks = nxticks
390        self.nyticks = nyticks
391        ## other useful arguments
392        ## ... not used here in ppplot but need to be attached to plot object
393        self.axisbg = "white"
394        self.superpose = False
395
396    # check
397    # -------------------------------
398    def check(self):
399        if self.f is None: print "!! ERROR !! Please define a field to be plotted" ; exit()
400        else: self.f = np.array(self.f) # ensure this is a numpy array
401
402    # define_from_var
403    # ... this uses settings in set_var.txt
404    # -------------------------------
405    def define_from_var(self):
406        if self.var is not None:
407         if self.var.upper() in vl.keys():
408          self.title = vl[self.var.upper()]
409          self.units = vu[self.var.upper()]
410
411    # make
412    # this is generic to all plots
413    # -------------------------------
414    def make(self):
415        self.check()
416        # labels, title, etc...
417        mpl.xlabel(self.xlabel)
418        mpl.ylabel(self.ylabel)
419        if self.swap:
420         if self.swaplab:
421           mpl.xlabel(self.ylabel)
422           mpl.ylabel(self.xlabel)
423        mpl.title(self.title)
424        # if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
425        if type(self.f).__name__ in 'MaskedArray':
426            self.f[self.f.mask] = self.f.fill_value
427
428################################
429# a subclass to make a 1D plot #
430################################
431class plot1d(plot):
432
433    # print out a help string when help is invoked on the object
434    # -------------------------------
435    def __repr__(self):
436        whatprint = 'plot1d object. \"help(plot1d)\" for more information\n'
437        return whatprint
438
439    # default settings
440    # -- user can define settings by two methods.
441    # -- 1. yeah = plot1d(title="foo")
442    # -- 2. yeah = pp() ; yeah.title = "foo"
443    # -------------------------------
444    def __init__(self,\
445                 linestyle=None,\
446                 color=None,\
447                 marker='x',\
448                 legend=None):
449        ## get initialization from parent class
450        plot.__init__(self)
451        ## what could be defined by the user
452        self.linestyle = linestyle
453        self.color = color
454        self.marker = marker
455        self.legend = legend
456
457    # define_from_var
458    # ... this uses settings in set_var.txt
459    # -------------------------------
460    def define_from_var(self):
461        # get what is done in the parent class
462        plot.define_from_var(self)
463        # add specific stuff
464        if self.var is not None:
465         if self.var.upper() in vl.keys():
466          self.ylabel = vl[self.var.upper()] + " (" + vu[self.var.upper()] + ")"
467          self.title = ""
468          self.fmt = vf[self.var.upper()]
469
470    # make
471    # -------------------------------
472    def make(self):
473        # get what is done in the parent class
474        plot.make(self)
475        if self.fmt is None: self.fmt = '%.0f'
476        # add specific stuff
477        mpl.grid(color='grey')
478        if self.linestyle == "": self.linestyle = " " # to allow for no line at all with ""
479        # set dummy x axis if not defined
480        if self.x is None: 
481            self.x = np.array(range(self.f.shape[0]))
482            print "!! WARNING !! dummy coordinates on x axis"
483        else:
484            self.x = np.array(self.x) # ensure this is a numpy array
485        # swapping if requested
486        if self.swap:  x = self.f ; y = self.x
487        else:          x = self.x ; y = self.f
488        # coefficients on axis
489        x=x*self.xcoeff ; y=y*self.ycoeff
490        # check axis
491        if x.size != y.size:
492            print "!! ERROR !! x and y sizes don't match on 1D plot.", x.size, y.size
493            exit()
494        # make the 1D plot
495        # either request linestyle or let matplotlib decide
496        if self.linestyle is not None and self.color is not None:
497            mpl.plot(x,y,self.color+self.linestyle,marker=self.marker,label=self.legend)
498        elif self.color is not None:
499            mpl.plot(x,y,color=self.color,marker=self.marker,label=self.legend)
500        elif self.linestyle is not None:
501            mpl.plot(x,y,linestyle=self.linestyle,marker=self.marker,label=self.legend)
502        else:
503            mpl.plot(x,y,marker=self.marker,label=self.legend)
504        # AXES
505        ax = mpl.gca()
506        # make log axes and/or invert ordinate
507        # ... this must be after plot so that axis bounds are well-defined
508        # ... also inverting must be after making the thing logarithmic
509        if self.logx: mpl.xscale("log") # not mpl.semilogx() because excludes log on y
510        if self.logy: mpl.yscale("log") # not mpl.semilogy() because excludes log on x
511        if self.invert: ax.set_ylim(ax.get_ylim()[::-1])
512        if self.xmin is not None and self.xmax is not None:
513          if self.xmin > self.xmax:
514            ax.set_xlim(ax.get_xlim()[::-1])
515            self.xmin,self.xmax = self.xmax,self.xmin
516        # add a label for line(s)
517        if self.legend is not None:
518            if self.legend != "":
519                mpl.legend(loc="best",fancybox=True)
520        # format labels
521        if self.swap: ax.xaxis.set_major_formatter(FormatStrFormatter(self.fmt))
522        else: ax.yaxis.set_major_formatter(FormatStrFormatter(self.fmt))
523        # plot limits: ensure that no curve would be outside the window
524        # x-axis
525        x1, x2 = ax.get_xbound()
526        xmin = ppcompute.min(x)
527        xmax = ppcompute.max(x)
528        if xmin < x1: x1 = xmin
529        if xmax > x2: x2 = xmax
530        if self.xmin is not None: x1 = self.xmin
531        if self.xmax is not None: x2 = self.xmax
532        ax.set_xbound(lower=x1,upper=x2)
533        # y-axis
534        y1, y2 = ax.get_ybound()
535        ymin = ppcompute.min(y)
536        ymax = ppcompute.max(y)
537        if ymin < y1: y1 = ymin
538        if ymax > y2: y2 = ymax
539        if self.ymin is not None: y1 = self.ymin
540        if self.ymax is not None: y2 = self.ymax
541        ax.set_ybound(lower=y1,upper=y2)
542        ## set with .div the number of ticks. (is it better than automatic?)
543        if not self.logx:
544            ax.xaxis.set_major_locator(MaxNLocator(self.div/2)) #TBD: with nxticks?
545        else:
546            print "!! WARNING. in logx mode, ticks are set automatically."
547        ## specific modulo labels
548        if self.modx is not None:
549            ax = labelmodulo(ax,self.modx)
550
551    # makeshow = make + show
552    # -------------------------------
553    def makeshow(self):
554        self.make()
555        mpl.show()
556
557################################
558# a subclass to make a 2D plot #
559################################
560class plot2d(plot):
561
562    # print out a help string when help is invoked on the object
563    # -------------------------------
564    def __repr__(self):
565        whatprint = 'plot2d object. \"help(plot2d)\" for more information\n'
566        return whatprint
567
568    # default settings
569    # -- user can define settings by two methods.
570    # -- 1. yeah = plot2d(title="foo")
571    # -- 2. yeah = pp() ; yeah.title = "foo"
572    # -------------------------------
573    def __init__(self,\
574                 y=None,\
575                 mapmode=False,\
576                 proj="cyl",\
577                 back=None,\
578                 trans=1.0,\
579                 vx=None,\
580                 vy=None,\
581                 c=None,\
582                 blon=None,\
583                 blat=None,\
584                 area=None,\
585                 vmin=None,\
586                 vmax=None,\
587                 showcb=True,\
588                 wscale=None,\
589                 svx=1,\
590                 svy=1,\
591                 leftcorrect=False,\
592                 colorvec="black"):
593        ## get initialization from parent class
594        plot.__init__(self)
595        ## what could be defined by the user
596        self.y = y
597        self.mapmode = mapmode
598        self.proj = proj
599        self.back = back
600        self.trans = trans
601        self.vx = vx
602        self.vy = vy
603        self.colorvec = colorvec
604        self.c = c
605        self.blon = blon ; self.blat = blat
606        self.area = area
607        self.vmin = vmin ; self.vmax = vmax
608        self.showcb = showcb
609        self.wscale = wscale
610        self.svx = svx
611        self.svy = svy
612        self.leftcorrect = leftcorrect
613
614    # define_from_var
615    # ... this uses settings in set_var.txt
616    # -------------------------------
617    def define_from_var(self):
618        # get what is done in the parent class
619        plot.define_from_var(self)
620        # add specific stuff
621        if self.var is not None:
622         if self.var.upper() in vl.keys():
623          self.colorbar = vc[self.var.upper()]
624          self.fmt = vf[self.var.upper()]
625
626    # make
627    # -------------------------------
628    def make(self):
629        # get what is done in the parent class...
630        plot.make(self)
631        if self.fmt is None: self.fmt = "%.1e"
632        # ... then add specific stuff
633        ############################################################################################
634        ### PRE-SETTINGS
635        ############################################################################################
636        # set dummy xy axis if not defined
637        if self.x is None: 
638            self.x = np.array(range(self.f.shape[0]))
639            self.mapmode = False
640            print "!! WARNING !! dummy coordinates on x axis"
641        if self.y is None: 
642            self.y = np.array(range(self.f.shape[1]))
643            self.mapmode = False
644            print "!! WARNING !! dummy coordinates on y axis"
645        # transposing if necessary
646        shape = self.f.shape
647        if shape[0] != shape[1]:
648         if len(self.x) == shape[0] and len(self.y) == shape[1]:
649            print "!! WARNING !! Transposing axes"
650            self.f = np.transpose(self.f)
651        # bound field
652        zevmin, zevmax = calculate_bounds(self.f,vmin=self.vmin,vmax=self.vmax)
653        what_I_plot = bounds(self.f,zevmin,zevmax)
654        # define contour field levels. define color palette
655        ticks = self.div + 1
656        zelevels = np.linspace(zevmin,zevmax,ticks)
657        palette = get_cmap(name=self.colorbar)
658        # do the same thing for possible contourline entries
659        if self.c is not None:
660            # if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
661            if type(self.c).__name__ in 'MaskedArray':
662               self.c[self.c.mask] = self.c.fill_value
663            zevminc, zevmaxc = calculate_bounds(self.c)
664            what_I_contour = bounds(self.c,zevminc,zevmaxc)
665            ticks = self.div + 1
666            zelevelsc = np.linspace(zevminc,zevmaxc,ticks)
667        ############################################################################################
668        ### MAIN PLOT
669        ### NB: contour lines are done before contour shades otherwise colorar error
670        ############################################################################################
671        if not self.mapmode:
672            ## A SIMPLE 2D PLOT
673            ###################
674            # swapping if requested
675            if self.swap:  x = self.y ; y = self.x
676            else:          x = self.x ; y = self.y
677            # coefficients on axis
678            x=x*self.xcoeff ; y=y*self.ycoeff
679            # make shaded and line contours
680            if self.c is not None: 
681                objC = mpl.contour(x, y, what_I_contour, \
682                            zelevelsc, colors = ccol, linewidths = cline)
683                #mpl.clabel(objC, inline=1, fontsize="small",\
684                #             inline_spacing=1,fmt="%.0f")
685            mpl.contourf(x, y, \
686                         self.f, \
687                         zelevels, cmap=palette)
688            #mpl.pcolor(x,y,\
689            #             self.f, \
690            #             cmap=palette)
691            # make log axes and/or invert ordinate
692            ax = mpl.gca()
693            if self.logx: mpl.semilogx()
694            if self.logy: mpl.semilogy()
695            if self.invert: ax.set_ylim(ax.get_ylim()[::-1])
696            if self.xmin is not None and self.xmax is not None:
697              if self.xmin > self.xmax: 
698                ax.set_xlim(ax.get_xlim()[::-1])
699                self.xmin,self.xmax = self.xmax,self.xmin
700            if self.xmin is not None: ax.set_xbound(lower=self.xmin)
701            if self.xmax is not None: ax.set_xbound(upper=self.xmax)
702            if self.ymin is not None: ax.set_ybound(lower=self.ymin)
703            if self.ymax is not None: ax.set_ybound(upper=self.ymax)
704            # set the number of ticks
705            if not self.logx:
706                ax.xaxis.set_major_locator(MaxNLocator(self.nxticks))
707            else:
708                print "!! WARNING. in logx mode, ticks are set automatically."
709            ## specific modulo labels
710            if self.modx is not None:
711                ax = labelmodulo(ax,self.modx)
712        else:
713            ## A 2D MAP USING PROJECTIONS (basemap)
714            #######################################
715            mpl.xlabel("") ; mpl.ylabel("")
716            # additional security in case self.proj is None here
717            # ... we set cylindrical projection (the simplest one)
718            if self.proj is None: self.proj = "cyl"
719            # get lon and lat in 2D version.
720            # (but first ensure we do have 2D coordinates)
721            if self.x.ndim == 1:        [self.x,self.y] = np.meshgrid(self.x,self.y)
722            elif self.x.ndim > 2:    print "!! ERROR !! lon and lat arrays must be 1D or 2D"
723            # get lat lon intervals and associated settings
724            wlon = [np.min(self.x),np.max(self.x)]
725            wlat = [np.min(self.y),np.max(self.y)]
726            # -- area presets are in set_area.txt
727            if self.area is not None:
728             if self.area in area.keys():
729                wlon, wlat = area[self.area]
730            # -- settings for meridians and parallels
731            steplon = int(abs(wlon[1]-wlon[0])/6.)
732            steplat = int(abs(wlat[1]-wlat[0])/3.)
733            #mertab = np.r_[wlon[0]:wlon[1]:steplon] ; merlab = [0,0,0,1]
734            #partab = np.r_[wlat[0]:wlat[1]:steplat] ; parlab = [1,0,0,0]
735            mertab = np.r_[-180.:180.:steplon] ; merlab = [0,0,0,1] #-360:360.
736            partab = np.r_[-90.:90.:steplat] ; parlab = [1,0,0,0]
737            format = '%.1f'
738            # -- center of domain and bounding lats
739            lon_0 = 0.5*(wlon[0]+wlon[1])
740            lat_0 = 0.5*(wlat[0]+wlat[1])
741            # some tests, bug fixes, and good-looking settings
742            # ... cyl is good for global and regional
743            if self.proj == "cyl":
744                format = '%.0f'
745            # ... global projections
746            elif self.proj in ["ortho","moll","robin"]:
747                wlat[0] = None ; wlat[1] = None ; wlon[0] = None ; wlon[1] = None
748                steplon = 30. ; steplat = 30.
749                if self.proj in ["moll"]: steplon = 60.
750                if self.proj in ["robin"]: steplon = 90.
751                mertab = np.r_[-360.:360.:steplon]
752                partab = np.r_[-90.:90.:steplat]
753                if self.proj == "ortho": 
754                    merlab = [0,0,0,0] ; parlab = [0,0,0,0]
755                    # in ortho projection, blon and blat can be used to set map center
756                    if self.blon is not None: lon_0 = self.blon
757                    if self.blat is not None: lat_0 = self.blat
758                elif self.proj == "moll":
759                    merlab = [0,0,0,0]
760                format = '%.0f'
761            # ... regional projections
762            elif self.proj in ["lcc","laea","merc"]:
763                if self.proj in ["lcc","laea"] and wlat[0] == -wlat[1]: 
764                    print "!! ERROR !! with Lambert lat1 must be different than lat2" ; exit()
765                if wlat[0] < -80. and wlat[1] > 80.:
766                    print "!! ERROR !! set an area (not global)" ; exit()
767                format = '%.0f'
768            elif self.proj in ["npstere","spstere"]:
769                # in polar projections, blat gives the bounding lat
770                # if not set, set something reasonable
771                if self.blat is None:   self.blat = 60.
772                # help the user who forgets self.blat would better be negative in spstere
773                # (this actually serves for the default setting just above)
774                if self.proj == "spstere" and self.blat > 0: self.blat = -self.blat
775                # labels
776                mertab = np.r_[-360.:360.:15.]
777                partab = np.r_[-90.:90.:5.]
778            # ... unsupported projections
779            else:
780                print "!! ERROR !! unsupported projection. supported: "+\
781                      "cyl, npstere, spstere, ortho, moll, robin, lcc, laea, merc"
782            # finally define projection
783            m = Basemap(projection=self.proj,\
784                        lat_0=lat_0,lon_0=lon_0,\
785                        boundinglat=self.blat,\
786                        llcrnrlat=wlat[0],urcrnrlat=wlat[1],\
787                        llcrnrlon=wlon[0],urcrnrlon=wlon[1])
788            # in some case need to translated to the left for colorbar + labels
789            # TBD: break stuff. a better solution should be found.
790            if self.leftcorrect:
791                ax = mpl.gca()
792                pos = ax.get_position().bounds
793                newpos = [0.,pos[1],pos[2],pos[3]]
794                ax.set_position(newpos)
795            # draw meridians and parallels
796            ft = int(mpl.rcParams['font.size']*3./4.)
797            zelatmax = 85.
798            m.drawmeridians(mertab,labels=merlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
799            m.drawparallels(partab,labels=parlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
800            # define background (see set_back.txt)
801            if self.back is not None:
802              if self.back in back.keys():
803                 print "**** info: loading a background, please wait.",self.back
804                 if self.back not in ["coast","sea"]:
805                    try: m.warpimage(back[self.back],scale=0.75)
806                    except: print "!! ERROR !! no background image could be loaded. probably not connected to the internet?"
807                 elif self.back == "coast":
808                    m.drawcoastlines()
809                 elif self.back == "sea":
810                    m.drawlsmask(land_color='white',ocean_color='aqua')
811              else:
812                 print "!! ERROR !! requested background not defined. change name or fill in set_back.txt" ; exit()
813            # define x and y given the projection
814            x, y = m(self.x, self.y)
815            # contour field. first line contour then shaded contour.
816            if self.c is not None: 
817                objC2 = m.contour(x, y, what_I_contour, \
818                            zelevelsc, colors = ccol, linewidths = cline)
819                #mpl.clabel(objC2, inline=1, fontsize=10,manual=True,fmt='-%2.0f$^{\circ}$C',colors='r')
820            m.contourf(x, y, what_I_plot, zelevels, cmap = palette, alpha = self.trans, antialiased=True)
821        ############################################################################################
822        ### COLORBAR
823        ############################################################################################
824        if self.trans > 0. and self.showcb:
825            ## draw colorbar. settings are different with projections. or if not mapmode.
826            if not self.mapmode: orientation=zeorientation ; frac = 0.075 ; pad = 0.03 ; lu = 0.5
827            elif self.proj in ['moll']: orientation="horizontal" ; frac = 0.08 ; pad = 0.03 ; lu = 1.0
828            elif self.proj in ['cyl']: orientation="vertical" ; frac = 0.023 ; pad = 0.03 ; lu = 0.5
829            else: orientation = zeorientation ; frac = zefrac ; pad = 0.03 ; lu = 0.5
830            zelevpal = np.linspace(zevmin,zevmax,num=min([ticks/2+1,21]))
831            zecb = mpl.colorbar(fraction=frac,pad=pad,\
832                                format=self.fmt,orientation=orientation,\
833                                ticks=zelevpal,\
834                                extend='neither',spacing='proportional')
835            if zeorientation == "horizontal": zecb.ax.set_xlabel(self.title) ; self.title = ""
836            # colorbar title --> units
837            if self.units not in ["dimless",""]:
838                zecb.ax.set_title("["+self.units+"]",fontsize=3.*mpl.rcParams['font.size']/4.,x=lu,y=1.025)
839
840        ############################################################################################
841        ### VECTORS. must be after the colorbar. we could also leave possibility for streamlines.
842        ############################################################################################
843        ### not expecting NaN in self.vx and self.vy. masked arrays is just enough.
844        if self.vx is not None and self.vy is not None: 
845                # vectors on map projection or simple 2D mapping
846                if self.mapmode: 
847                   [vecx,vecy] = m.rotate_vector(self.vx,self.vy,self.x,self.y) # for metwinds only ?
848                else:
849                   vecx,vecy = self.vx,self.vy
850                   if x.ndim < 2 and y.ndim < 2: x,y = np.meshgrid(x,y)
851                # reference vector is scaled
852                if self.wscale is None: 
853                    self.wscale = ppcompute.mean(np.sqrt(self.vx*self.vx+self.vy*self.vy))
854                # make vector field
855                if self.mapmode: 
856                    q = m.quiver( x[::self.svy,::self.svx],y[::self.svy,::self.svx],\
857                                  vecx[::self.svy,::self.svx],vecy[::self.svy,::self.svx],\
858                                  angles='xy',color=self.colorvec,pivot='tail',\
859                                  scale=self.wscale*reducevec,width=widthvec )
860                else:
861                    q = mpl.quiver( x[::self.svy,::self.svx],y[::self.svy,::self.svx],\
862                                    vecx[::self.svy,::self.svx],vecy[::self.svy,::self.svx],\
863                                    angles='xy',color=self.colorvec,pivot='tail',\
864                                    scale=self.wscale*reducevec,width=widthvec )
865                # make vector key.
866                #keyh = 1.025 ; keyv = 1.05 # upper right corner over colorbar
867                keyh = 0.97 ; keyv = 1.05
868                #keyh = -0.03 ; keyv = 1.08 # upper left corner
869                p = mpl.quiverkey(q,keyh,keyv,\
870                                  self.wscale,str(int(self.wscale)),\
871                                  fontproperties={'size': 'small'},\
872                                  color='black',labelpos='S',labelsep = 0.07)
873        ############################################################################################
874        ### TEXT. ANYWHERE. add_text.txt should be present with lines x ; y ; text ; color
875        ############################################################################################
876        try:
877            f = open("add_text.txt", 'r')
878            for line in f:
879              if "#" in line: pass
880              else:
881                  userx, usery, usert, userc = line.strip().split(';')
882                  userc = userc.strip()
883                  usert = usert.strip()
884                  userx = float(userx.strip())
885                  usery = float(usery.strip())
886                  if self.mapmode: userx,usery = m(userx,usery)
887                  mpl.text(userx,usery,usert,\
888                           color = userc,\
889                           horizontalalignment='center',\
890                           verticalalignment='center')
891            f.close()
892        except IOError:
893            pass
894
895    # makeshow = make + show
896    # -------------------------------
897    def makeshow(self):
898        self.make()
899        mpl.show()
Note: See TracBrowser for help on using the repository browser.