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 '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 api_onelevel ( path_to_input = None, \ input_name = 'wrfout_d0?_????-??-??_??:00:00', \ path_to_output = None, \ output_name = 'output.nc', \ process = 'list', \ fields = 'tk,W,uvmet,HGT', \ debug = False, \ bit64 = False, \ oldvar = True, \ interp_method = 4, \ extrapolate = 0, \ unstagger_grid = False, \ onelevel = 0.020 ): import api import numpy as np if not path_to_input: path_to_input = './' if not path_to_output: path_to_output = path_to_input api.api_main ( path_to_input, input_name, path_to_output, output_name, \ process, fields, debug, bit64, oldvar, np.arange (299), \ interp_method, extrapolate, unstagger_grid, onelevel ) return 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 return [[lon2d[0,0],lon2d[nx,ny]],[lat2d[0,0],lat2d[nx,ny]]] 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 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.max(x) + np.std(x) / 3. ## pb pour les domaines globaux ... #posy = np.mean(y) #posx = np.min(x) #posy = np.max(x) #posx = np.max(x) - np.std(x) / 10. #posy = np.max(y) + np.std(y) / 10. 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=='white': kcolor='black' elif color=='yellow': kcolor=color 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 = 2000. 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 == "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"]: step = findstep(wlon) else: step = 10. 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 marsmap (whichone="vishires"): whichlink = { \ "vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ "vishires": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/MarsMap_2500x1250.jpg",\ "mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ "molabw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/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