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

Last change on this file since 1142 was 1131, checked in by aslmd, 12 years ago

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