source: trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py @ 226

def latinterv (area):
        if   area == "Europe": 
                wlat = [20.,80.]
                wlon = [-50.,50.]
        elif area == "Central_America":
                wlat = [-10.,40.]
                wlon = [230.,300.]
        elif area == "Africa":
                wlat = [-20.,50.]
                wlon = [-50.,50.]
        elif area == "Whole":
                wlat = [-90.,90.]
                wlon = [-180.,180.]
        elif area == "Southern_Hemisphere":
                wlat = [-90.,60.]
                wlon = [-180.,180.]
        elif area == "Northern_Hemisphere":
                wlat = [-60.,90.]
                wlon = [-180.,180.]
        elif area == "Tharsis":
                wlat = [-30.,60.]
                wlon = [-170.,-10.]
        elif area == "Whole_No_High":
                wlat = [-60.,60.]
                wlon = [-180.,180.]
        elif area == "Chryse":
                wlat = [-60.,60.]
                wlon = [-60.,60.]
        elif area == "North_Pole":
                wlat = [60.,90.]
                wlon = [-180.,180.]
        elif area == "Close_North_Pole":
                wlat = [75.,90.]
                wlon = [-180.,180.]
        return wlon,wlat

#def landers (map) 
#    map.plot(blue_calf_lon,blue_calf_lat, 'gs')
#    return

def getlschar ( namefile ):
    #### strangely enough this does not work for api or ncrcat results!
    from netCDF4 import Dataset
    from timestuff import sol2ls
    nc  = Dataset(namefile)
    if namefile[0]+namefile[1]+namefile[2] != "geo" \
       and 'Times' in nc.variables \
       and 'vert' not in nc.variables:
        zetime = nc.variables['Times'][0]
        zetimestart = getattr(nc, 'START_DATE')
        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
        if zeday < 0:    lschar=""  ## might have crossed a month... fix soon
        else:            lschar="_Ls"+str( int( sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) )
        ###
        zetime2 = nc.variables['Times'][1]
        one  = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
        next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. 
        zehour    = one
        zehourin  = abs ( next - one )
    else:
        lschar=""
        zehour = 0
        zehourin = 1  
    return lschar, zehour, zehourin

def getprefix (nc):
    prefix = 'LMD_MMM_'
    prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
    prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
    return prefix

def getproj (nc):
    map_proj = getattr(nc, 'MAP_PROJ')
    cen_lat  = getattr(nc, 'CEN_LAT')
    if map_proj == 2:
        if cen_lat > 10.:    
            proj="npstere"
            print "NP stereographic polar domain" 
        else:            
            proj="spstere"
            print "SP stereographic polar domain"
    elif map_proj == 1: 
        print "lambert projection domain" 
        proj="lcc"
    elif map_proj == 3: 
        print "mercator projection"
        proj="merc"
    else:
        proj="merc"
    return proj    

def ptitle (name):
    from matplotlib.pyplot import title
    title(name)
    print name

def simplinterv (lon2d,lat2d):
    import numpy as np
    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]

def wrfinterv (lon2d,lat2d):
    nx = len(lon2d[0,:])-1
    ny = len(lon2d[:,0])-1
    lon1 = lon2d[0,0] 
    lon2 = lon2d[nx,ny] 
    lat1 = lat2d[0,0] 
    lat2 = lat2d[nx,ny] 
    wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
    if lon1 < lon2:  wlon = [lon1, lon2 + wider]  ## a tester en normal
    else:            wlon = [lon2, lon1 + wider]
    if lat1 < lat2:  wlat = [lat1, lat2]
    else:            wlat = [lat2, lat1]
    return [wlon,wlat]

def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True):
    import  matplotlib.pyplot as plt
    res = int(res)
    name = filename+"_"+str(res)+".png"
    if folder != '':      name = folder+'/'+name
    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
    if disp:              display(name)         
    return

def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder=''):
    makeplotpngres(filename,minres,     pad_inches_value=pad_inches_value,folder=folder)
    makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False)
    return

def dumpbdy (field):
    nx = len(field[0,:])-1
    ny = len(field[:,0])-1
    return field[5:ny-5,5:nx-5]

def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
    import numpy as np
    if is1d:
        lat = nc.variables[nlat][:]
        lon = nc.variables[nlon][:]
        [lon2d,lat2d] = np.meshgrid(lon,lat)
    else:
        lat = nc.variables[nlat][0,:,:]
        lon = nc.variables[nlon][0,:,:]
        [lon2d,lat2d] = [lon,lat]
    return lon2d,lat2d

def smooth (field, coeff):
        ## actually blur_image could work with different coeff on x and y
        if coeff > 1:   result = blur_image(field,int(coeff))
        else:           result = field
        return result

def gauss_kern(size, sizey=None):
        import numpy as np
        ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth     
        # Returns a normalized 2D gauss kernel array for convolutions
        size = int(size)
        if not sizey:
                sizey = size
        else:
                sizey = int(sizey)
        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
        return g / g.sum()

def blur_image(im, n, ny=None) :
        from scipy.signal import convolve
        ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
        # blurs the image by convolving with a gaussian kernel of typical size n. 
        # The optional keyword argument ny allows for a different size in the y direction.
        g = gauss_kern(n, sizey=ny)
        improc = convolve(im, g, mode='same')
        return improc

def getwinds (nc,charu='Um',charv='Vm'):
    import numpy as np
    u = nc.variables[charu]
    v = nc.variables[charv]
    if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1]
    if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :]
                     ### ou alors prendre les coordonnees speciales
    return u,v

def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
    ## scale regle la reference du vecteur
    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
    import  matplotlib.pyplot               as plt
    import  numpy                           as np
    posx = np.min(x) - np.std(x) / 10.
    posy = np.min(y) - np.std(y) / 10.
    u = smooth(u,csmooth)
    v = smooth(v,csmooth)
    widthvec = 0.003 #0.005 #0.003
    q = plt.quiver( x[::stride,::stride],\
                    y[::stride,::stride],\
                    u[::stride,::stride],\
                    v[::stride,::stride],\
                    angles='xy',color=color,\
                    scale=factor,width=widthvec )
    if color in ['white','yellow']:     kcolor='black'
    else:                               kcolor=color
    if key: p = plt.quiverkey(q,posx,posy,scale,\
                   str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
    return 

def display (name):
    from os import system
    system("display "+name+" > /dev/null 2> /dev/null &")
    return name

def findstep (wlon):
    steplon = int((wlon[1]-wlon[0])/4.)  #3
    step = 120.
    while step > steplon and step > 15. :       step = step / 2.
    if step <= 15.:
        while step > steplon and step > 5.  :   step = step - 5.
    if step <= 5.:
        while step > steplon and step > 1.  :   step = step - 1.
    if step <= 1.:
        step = 1. 
    return step

def define_proj (char,wlon,wlat,back="."):
    from    mpl_toolkits.basemap            import Basemap
    import  numpy                           as np
    import  matplotlib                      as mpl
    meanlon = 0.5*(wlon[0]+wlon[1])
    meanlat = 0.5*(wlat[0]+wlat[1])
    if   wlat[0] >= 80.:   blat =  40. 
    elif wlat[0] <= -80.:  blat = -40. 
    else:                  blat = wlat[0]
    h = 50.  ## en km
    radius = 3397200.
    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat)
    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.)
    elif char == "nplaea":  m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
    elif char == "laea":    m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
    if char in ["cyl","lcc","merc","nsper","laea"]:   step = findstep(wlon)
    else:                                             step = 10.
    print step
    m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
    if back == ".":      m.warpimage(marsmap(),scale=0.75)
    elif back == None:   pass 
    else:                m.warpimage(marsmap(back),scale=0.75)
    return m

def fmtvar (whichvar="def"):
    fmtvar    =     { \
             "tk":           "%.0f",\
             "tpot":         "%.0f",\
             "def":          "%.1e",\
             "PTOT":         "%.0f",\
             "HGT":          "%.1e",\
             "USTM":         "%.2f",\
             "HFX":          "%.0f",\
                    }
    if whichvar not in fmtvar:
        whichvar = "def"
    return fmtvar[whichvar]

def defcolorb (whichone="def"):
    whichcolorb =    { \
             "def":          "spectral",\
             "HGT":          "spectral",\
             "tk":           "gist_heat",\
             "QH2O":         "PuBu",\
             "USTM":         "YlOrRd",\
             "HFX":          "RdYlBu",\
                     }
    if whichone not in whichcolorb:
        whichone = "def"
    return whichcolorb[whichone]

def definecolorvec (whichone="def"):
        whichcolor =    { \
                "def":          "black",\
                "vis":          "yellow",\
                "vishires":     "yellow",\
                "molabw":       "yellow",\
                "mola":         "black",\
                "gist_heat":    "white",\
                "hot":          "tk",\
                "gist_rainbow": "black",\
                "spectral":     "black",\
                "gray":         "red",\
                "PuBu":         "black",\
                        }
        if whichone not in whichcolor:
                whichone = "def"
        return whichcolor[whichone]

def marsmap (whichone="vishires"):
        whichlink =     { \
                "vis":          "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
                "vishires":     "http://dl.dropbox.com/u/11078310/MarsMap_2500x1250.jpg",\
                "mola":         "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
                "molabw":       "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
                        }
        if whichone not in whichlink: 
                print "marsmap: choice not defined... you'll get the default one... "
                whichone = "vishires"  
        return whichlink[whichone]

def earthmap (whichone):
        if   whichone == "contrast":    whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
        elif whichone == "bw":          whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
        elif whichone == "nice":        whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
        return whichlink