source: trunk/MESOSCALE/PLOT/PYTHON/mylib/myplot.py @ 194

Last change on this file since 194 was 194, checked in by aslmd, 14 years ago

MESOSCALE: corrected graphics script to display Ls. and added a figure in the first page of the user manual.

File size: 9.7 KB
Line 
1def latinterv (area):
2        if   area == "Europe": 
3                wlat = [20.,80.]
4                wlon = [-50.,50.]
5        elif area == "Central_America":
6                wlat = [-10.,40.]
7                wlon = [230.,300.]
8        elif area == "Africa":
9                wlat = [-20.,50.]
10                wlon = [-50.,50.]
11        elif area == "Whole":
12                wlat = [-90.,90.]
13                wlon = [-180.,180.]
14        elif area == "Southern_Hemisphere":
15                wlat = [-90.,60.]
16                wlon = [-180.,180.]
17        elif area == "Northern_Hemisphere":
18                wlat = [-60.,90.]
19                wlon = [-180.,180.]
20        elif area == "Tharsis":
21                wlat = [-30.,60.]
22                wlon = [-170.,-10.]
23        elif area == "Whole_No_High":
24                wlat = [-60.,60.]
25                wlon = [-180.,180.]
26        elif area == "Chryse":
27                wlat = [-60.,60.]
28                wlon = [-60.,60.]
29        elif area == "North_Pole":
30                wlat = [60.,90.]
31                wlon = [-180.,180.]
32        elif area == "Close_North_Pole":
33                wlat = [75.,90.]
34                wlon = [-180.,180.]
35        return wlon,wlat
36
37#def landers (map)
38#    map.plot(blue_calf_lon,blue_calf_lat, 'gs')
39#    return
40
41def api_onelevel (  path_to_input   = None, \
42                    input_name      = 'wrfout_d0?_????-??-??_??:00:00', \
43                    path_to_output  = None, \
44                    output_name     = 'output.nc', \
45                    process         = 'list', \
46                    fields          = 'tk,W,uvmet,HGT', \
47                    debug           = False, \
48                    bit64           = False, \
49                    oldvar          = True, \
50                    interp_method   = 4, \
51                    extrapolate     = 0, \
52                    unstagger_grid  = False, \
53                    onelevel        = 0.020  ):
54    import api
55    import numpy as np
56    if not path_to_input:   path_to_input  = './'
57    if not path_to_output:  path_to_output = path_to_input
58    api.api_main ( path_to_input, input_name, path_to_output, output_name, \
59                   process, fields, debug, bit64, oldvar, np.arange (299), \
60                   interp_method, extrapolate, unstagger_grid, onelevel )
61    return
62
63def getproj (nc):
64    map_proj = getattr(nc, 'MAP_PROJ')
65    cen_lat  = getattr(nc, 'CEN_LAT')
66    if map_proj == 2:
67        if cen_lat > 10.:   
68            proj="npstere"
69            print "NP stereographic polar domain" 
70        else:           
71            proj="spstere"
72            print "SP stereographic polar domain"
73    elif map_proj == 1: 
74        print "lambert projection domain" 
75        proj="lcc"
76    elif map_proj == 3: 
77        print "mercator projection"
78        proj="merc"
79    else:
80        proj="merc"
81    return proj   
82
83def ptitle (name):
84    from matplotlib.pyplot import title
85    title(name)
86    print name
87
88def simplinterv (lon2d,lat2d):
89    import numpy as np
90    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
91
92def wrfinterv (lon2d,lat2d):
93    nx = len(lon2d[0,:])-1
94    ny = len(lon2d[:,0])-1
95    return [[lon2d[0,0],lon2d[nx,ny]],[lat2d[0,0],lat2d[nx,ny]]]
96
97def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True):
98    import  matplotlib.pyplot as plt
99    res = int(res)
100    name = filename+"_"+str(res)+".png"
101    if folder != '':      name = folder+'/'+name
102    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
103    if disp:              display(name)         
104    return
105
106def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder=''):
107    makeplotpngres(filename,minres,     pad_inches_value=pad_inches_value,folder=folder)
108    makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False)
109    return
110
111def dumpbdy (field):
112    nx = len(field[0,:])-1
113    ny = len(field[:,0])-1
114    return field[5:ny-5,5:nx-5]
115
116def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
117    import numpy as np
118    if is1d:
119        lat = nc.variables[nlat][:]
120        lon = nc.variables[nlon][:]
121        [lon2d,lat2d] = np.meshgrid(lon,lat)
122    else:
123        lat = nc.variables[nlat][0,:,:]
124        lon = nc.variables[nlon][0,:,:]
125        [lon2d,lat2d] = [lon,lat]
126    return lon2d,lat2d
127
128def smooth (field, coeff):
129        ## actually blur_image could work with different coeff on x and y
130        if coeff > 1:   result = blur_image(field,int(coeff))
131        else:           result = field
132        return result
133
134def gauss_kern(size, sizey=None):
135        import numpy as np
136        ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth     
137        # Returns a normalized 2D gauss kernel array for convolutions
138        size = int(size)
139        if not sizey:
140                sizey = size
141        else:
142                sizey = int(sizey)
143        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
144        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
145        return g / g.sum()
146
147def blur_image(im, n, ny=None) :
148        from scipy.signal import convolve
149        ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
150        # blurs the image by convolving with a gaussian kernel of typical size n.
151        # The optional keyword argument ny allows for a different size in the y direction.
152        g = gauss_kern(n, sizey=ny)
153        improc = convolve(im, g, mode='same')
154        return improc
155
156def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
157    ## scale regle la reference du vecteur
158    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
159    import  matplotlib.pyplot               as plt
160    import  numpy                           as np
161    #posx = np.max(x) + np.std(x) / 3.  ## pb pour les domaines globaux ...
162    #posy = np.mean(y)
163    #posx = np.min(x)
164    #posy = np.max(x)
165    #posx = np.max(x) - np.std(x) / 10.
166    #posy = np.max(y) + np.std(y) / 10.
167    posx = np.min(x) - np.std(x) / 10.
168    posy = np.min(y) - np.std(y) / 10.
169    u = smooth(u,csmooth)
170    v = smooth(v,csmooth)
171    widthvec = 0.003 #0.005 #0.003
172    q = plt.quiver( x[::stride,::stride],\
173                    y[::stride,::stride],\
174                    u[::stride,::stride],\
175                    v[::stride,::stride],\
176                    angles='xy',color=color,\
177                    scale=factor,width=widthvec )
178    if color=='white':     kcolor='black'
179    elif color=='yellow':  kcolor=color
180    else:                  kcolor=color
181    if key: p = plt.quiverkey(q,posx,posy,scale,\
182                   str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
183    return 
184
185def display (name):
186    from os import system
187    system("display "+name+" > /dev/null 2> /dev/null &")
188    return name
189
190def findstep (wlon):
191    steplon = int((wlon[1]-wlon[0])/4.)  #3
192    step = 120.
193    while step > steplon and step > 15. :       step = step / 2.
194    if step <= 15.:
195        while step > steplon and step > 5.  :   step = step - 5.
196    if step <= 5.:
197        while step > steplon and step > 1.  :   step = step - 1.
198    if step <= 1.:
199        step = 1. 
200    return step
201
202def define_proj (char,wlon,wlat,back="."):
203    from    mpl_toolkits.basemap            import Basemap
204    import  numpy                           as np
205    import  matplotlib                      as mpl
206    meanlon = 0.5*(wlon[0]+wlon[1])
207    meanlat = 0.5*(wlat[0]+wlat[1])
208    if   wlat[0] >= 80.:   blat =  40. 
209    elif wlat[0] <= -80.:  blat = -40. 
210    else:                  blat = wlat[0]
211    h = 2000.
212    radius = 3397200
213    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
214                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
215    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
216    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat)
217    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
218                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
219    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
220    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.)
221    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
222    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
223                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
224    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
225    if char in ["cyl","lcc","merc"]:   step = findstep(wlon)
226    else:                              step = 10.
227    m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
228    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
229    if back == ".":      m.warpimage(marsmap(),scale=0.75)
230    elif back == None:   pass 
231    else:                m.warpimage(marsmap(back),scale=0.75)
232    return m
233
234def marsmap (whichone="vishires"):
235        whichlink =     { \
236                "vis":          "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
237                "vishires":     "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/MarsMap_2500x1250.jpg",\
238                "mola":         "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
239                "molabw":       "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/MarsElevation_2500x1250.jpg",\
240                        }
241        if whichone not in whichlink: 
242                print "marsmap: choice not defined... you'll get the default one... "
243                whichone = "vishires" 
244        return whichlink[whichone]
245
246def earthmap (whichone):
247        if   whichone == "contrast":    whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
248        elif whichone == "bw":          whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
249        elif whichone == "nice":        whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
250        return whichlink
251
Note: See TracBrowser for help on using the repository browser.