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

Last change on this file since 1031 was 1029, checked in by aslmd, 12 years ago

UTIL PYTHON planetoplot_v2

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