| 1 | def 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 | |
|---|
| 41 | def getlschar ( namefile ): |
|---|
| 42 | #### strangely enough this does not work for api or ncrcat results! |
|---|
| 43 | from netCDF4 import Dataset |
|---|
| 44 | from timestuff import sol2ls |
|---|
| 45 | nc = Dataset(namefile) |
|---|
| 46 | if 'Times' in nc.variables and 'vert' not in nc.variables: |
|---|
| 47 | zetime = nc.variables['Times'][0] |
|---|
| 48 | zetimestart = getattr(nc, 'START_DATE') |
|---|
| 49 | zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) |
|---|
| 50 | if zeday < 0: lschar="" ## might have crossed a month... fix soon |
|---|
| 51 | else: lschar="_Ls"+str( int( sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) ) |
|---|
| 52 | ### |
|---|
| 53 | zetime2 = nc.variables['Times'][1] |
|---|
| 54 | one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37. |
|---|
| 55 | next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. |
|---|
| 56 | zehour = one |
|---|
| 57 | zehourin = abs ( next - one ) |
|---|
| 58 | else: |
|---|
| 59 | lschar="" |
|---|
| 60 | zehour = 0 |
|---|
| 61 | zehourin = 1 |
|---|
| 62 | return lschar, zehour, zehourin |
|---|
| 63 | |
|---|
| 64 | def getprefix (nc): |
|---|
| 65 | prefix = 'LMD_MMM_' |
|---|
| 66 | prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_' |
|---|
| 67 | prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_' |
|---|
| 68 | return prefix |
|---|
| 69 | |
|---|
| 70 | def api_onelevel ( path_to_input = None, \ |
|---|
| 71 | input_name = 'wrfout_d0?_????-??-??_??:00:00', \ |
|---|
| 72 | path_to_output = None, \ |
|---|
| 73 | output_name = 'output.nc', \ |
|---|
| 74 | process = 'list', \ |
|---|
| 75 | fields = 'tk,W,uvmet,HGT', \ |
|---|
| 76 | debug = False, \ |
|---|
| 77 | bit64 = False, \ |
|---|
| 78 | oldvar = True, \ |
|---|
| 79 | interp_method = 4, \ |
|---|
| 80 | extrapolate = 0, \ |
|---|
| 81 | unstagger_grid = False, \ |
|---|
| 82 | onelevel = 0.020 ): |
|---|
| 83 | import api |
|---|
| 84 | import numpy as np |
|---|
| 85 | if not path_to_input: path_to_input = './' |
|---|
| 86 | if not path_to_output: path_to_output = path_to_input |
|---|
| 87 | api.api_main ( path_to_input, input_name, path_to_output, output_name, \ |
|---|
| 88 | process, fields, debug, bit64, oldvar, np.arange (299), \ |
|---|
| 89 | interp_method, extrapolate, unstagger_grid, onelevel ) |
|---|
| 90 | return |
|---|
| 91 | |
|---|
| 92 | def getproj (nc): |
|---|
| 93 | map_proj = getattr(nc, 'MAP_PROJ') |
|---|
| 94 | cen_lat = getattr(nc, 'CEN_LAT') |
|---|
| 95 | if map_proj == 2: |
|---|
| 96 | if cen_lat > 10.: |
|---|
| 97 | proj="npstere" |
|---|
| 98 | print "NP stereographic polar domain" |
|---|
| 99 | else: |
|---|
| 100 | proj="spstere" |
|---|
| 101 | print "SP stereographic polar domain" |
|---|
| 102 | elif map_proj == 1: |
|---|
| 103 | print "lambert projection domain" |
|---|
| 104 | proj="lcc" |
|---|
| 105 | elif map_proj == 3: |
|---|
| 106 | print "mercator projection" |
|---|
| 107 | proj="merc" |
|---|
| 108 | else: |
|---|
| 109 | proj="merc" |
|---|
| 110 | return proj |
|---|
| 111 | |
|---|
| 112 | def ptitle (name): |
|---|
| 113 | from matplotlib.pyplot import title |
|---|
| 114 | title(name) |
|---|
| 115 | print name |
|---|
| 116 | |
|---|
| 117 | def simplinterv (lon2d,lat2d): |
|---|
| 118 | import numpy as np |
|---|
| 119 | return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]] |
|---|
| 120 | |
|---|
| 121 | def wrfinterv (lon2d,lat2d): |
|---|
| 122 | nx = len(lon2d[0,:])-1 |
|---|
| 123 | ny = len(lon2d[:,0])-1 |
|---|
| 124 | return [[lon2d[0,0],lon2d[nx,ny]],[lat2d[0,0],lat2d[nx,ny]]] |
|---|
| 125 | |
|---|
| 126 | def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True): |
|---|
| 127 | import matplotlib.pyplot as plt |
|---|
| 128 | res = int(res) |
|---|
| 129 | name = filename+"_"+str(res)+".png" |
|---|
| 130 | if folder != '': name = folder+'/'+name |
|---|
| 131 | plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value) |
|---|
| 132 | if disp: display(name) |
|---|
| 133 | return |
|---|
| 134 | |
|---|
| 135 | def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder=''): |
|---|
| 136 | makeplotpngres(filename,minres, pad_inches_value=pad_inches_value,folder=folder) |
|---|
| 137 | makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False) |
|---|
| 138 | return |
|---|
| 139 | |
|---|
| 140 | def dumpbdy (field): |
|---|
| 141 | nx = len(field[0,:])-1 |
|---|
| 142 | ny = len(field[:,0])-1 |
|---|
| 143 | return field[5:ny-5,5:nx-5] |
|---|
| 144 | |
|---|
| 145 | def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): |
|---|
| 146 | import numpy as np |
|---|
| 147 | if is1d: |
|---|
| 148 | lat = nc.variables[nlat][:] |
|---|
| 149 | lon = nc.variables[nlon][:] |
|---|
| 150 | [lon2d,lat2d] = np.meshgrid(lon,lat) |
|---|
| 151 | else: |
|---|
| 152 | lat = nc.variables[nlat][0,:,:] |
|---|
| 153 | lon = nc.variables[nlon][0,:,:] |
|---|
| 154 | [lon2d,lat2d] = [lon,lat] |
|---|
| 155 | return lon2d,lat2d |
|---|
| 156 | |
|---|
| 157 | def smooth (field, coeff): |
|---|
| 158 | ## actually blur_image could work with different coeff on x and y |
|---|
| 159 | if coeff > 1: result = blur_image(field,int(coeff)) |
|---|
| 160 | else: result = field |
|---|
| 161 | return result |
|---|
| 162 | |
|---|
| 163 | def gauss_kern(size, sizey=None): |
|---|
| 164 | import numpy as np |
|---|
| 165 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
|---|
| 166 | # Returns a normalized 2D gauss kernel array for convolutions |
|---|
| 167 | size = int(size) |
|---|
| 168 | if not sizey: |
|---|
| 169 | sizey = size |
|---|
| 170 | else: |
|---|
| 171 | sizey = int(sizey) |
|---|
| 172 | x, y = np.mgrid[-size:size+1, -sizey:sizey+1] |
|---|
| 173 | g = np.exp(-(x**2/float(size)+y**2/float(sizey))) |
|---|
| 174 | return g / g.sum() |
|---|
| 175 | |
|---|
| 176 | def blur_image(im, n, ny=None) : |
|---|
| 177 | from scipy.signal import convolve |
|---|
| 178 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
|---|
| 179 | # blurs the image by convolving with a gaussian kernel of typical size n. |
|---|
| 180 | # The optional keyword argument ny allows for a different size in the y direction. |
|---|
| 181 | g = gauss_kern(n, sizey=ny) |
|---|
| 182 | improc = convolve(im, g, mode='same') |
|---|
| 183 | return improc |
|---|
| 184 | |
|---|
| 185 | def getwinds (nc,charu='Um',charv='Vm'): |
|---|
| 186 | import numpy as np |
|---|
| 187 | u = nc.variables[charu] |
|---|
| 188 | v = nc.variables[charv] |
|---|
| 189 | if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1] |
|---|
| 190 | if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :] |
|---|
| 191 | ### ou alors prendre les coordonnees speciales |
|---|
| 192 | return u,v |
|---|
| 193 | |
|---|
| 194 | def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True): |
|---|
| 195 | ## scale regle la reference du vecteur |
|---|
| 196 | ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs. |
|---|
| 197 | import matplotlib.pyplot as plt |
|---|
| 198 | import numpy as np |
|---|
| 199 | posx = np.min(x) - np.std(x) / 10. |
|---|
| 200 | posy = np.min(y) - np.std(y) / 10. |
|---|
| 201 | u = smooth(u,csmooth) |
|---|
| 202 | v = smooth(v,csmooth) |
|---|
| 203 | widthvec = 0.003 #0.005 #0.003 |
|---|
| 204 | q = plt.quiver( x[::stride,::stride],\ |
|---|
| 205 | y[::stride,::stride],\ |
|---|
| 206 | u[::stride,::stride],\ |
|---|
| 207 | v[::stride,::stride],\ |
|---|
| 208 | angles='xy',color=color,\ |
|---|
| 209 | scale=factor,width=widthvec ) |
|---|
| 210 | if color in ['white','yellow']: kcolor='black' |
|---|
| 211 | else: kcolor=color |
|---|
| 212 | if key: p = plt.quiverkey(q,posx,posy,scale,\ |
|---|
| 213 | str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03) |
|---|
| 214 | return |
|---|
| 215 | |
|---|
| 216 | def display (name): |
|---|
| 217 | from os import system |
|---|
| 218 | system("display "+name+" > /dev/null 2> /dev/null &") |
|---|
| 219 | return name |
|---|
| 220 | |
|---|
| 221 | def findstep (wlon): |
|---|
| 222 | steplon = int((wlon[1]-wlon[0])/4.) #3 |
|---|
| 223 | step = 120. |
|---|
| 224 | while step > steplon and step > 15. : step = step / 2. |
|---|
| 225 | if step <= 15.: |
|---|
| 226 | while step > steplon and step > 5. : step = step - 5. |
|---|
| 227 | if step <= 5.: |
|---|
| 228 | while step > steplon and step > 1. : step = step - 1. |
|---|
| 229 | if step <= 1.: |
|---|
| 230 | step = 1. |
|---|
| 231 | return step |
|---|
| 232 | |
|---|
| 233 | def define_proj (char,wlon,wlat,back="."): |
|---|
| 234 | from mpl_toolkits.basemap import Basemap |
|---|
| 235 | import numpy as np |
|---|
| 236 | import matplotlib as mpl |
|---|
| 237 | meanlon = 0.5*(wlon[0]+wlon[1]) |
|---|
| 238 | meanlat = 0.5*(wlat[0]+wlat[1]) |
|---|
| 239 | if wlat[0] >= 80.: blat = 40. |
|---|
| 240 | elif wlat[0] <= -80.: blat = -40. |
|---|
| 241 | else: blat = wlat[0] |
|---|
| 242 | h = 2000. |
|---|
| 243 | radius = 3397200. |
|---|
| 244 | if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\ |
|---|
| 245 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 246 | elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon) |
|---|
| 247 | elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat) |
|---|
| 248 | elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\ |
|---|
| 249 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 250 | elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.) |
|---|
| 251 | elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.) |
|---|
| 252 | elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.) |
|---|
| 253 | elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\ |
|---|
| 254 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 255 | fontsizemer = int(mpl.rcParams['font.size']*3./4.) |
|---|
| 256 | if char in ["cyl","lcc","merc"]: step = findstep(wlon) |
|---|
| 257 | else: step = 10. |
|---|
| 258 | m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer) |
|---|
| 259 | m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer) |
|---|
| 260 | if back == ".": m.warpimage(marsmap(),scale=0.75) |
|---|
| 261 | elif back == None: pass |
|---|
| 262 | else: m.warpimage(marsmap(back),scale=0.75) |
|---|
| 263 | return m |
|---|
| 264 | |
|---|
| 265 | def fmtvar (whichvar="def"): |
|---|
| 266 | fmtvar = { \ |
|---|
| 267 | "tk": "%.0f",\ |
|---|
| 268 | "tpot": "%.0f",\ |
|---|
| 269 | "def": "%.1e",\ |
|---|
| 270 | "PTOT": "%.0f",\ |
|---|
| 271 | "HGT": "%.1e",\ |
|---|
| 272 | "USTM": "%.2f",\ |
|---|
| 273 | } |
|---|
| 274 | if whichvar not in fmtvar: |
|---|
| 275 | whichvar = "def" |
|---|
| 276 | return fmtvar[whichvar] |
|---|
| 277 | |
|---|
| 278 | def defcolorb (whichone="def"): |
|---|
| 279 | whichcolorb = { \ |
|---|
| 280 | "def": "spectral",\ |
|---|
| 281 | "HGT": "spectral",\ |
|---|
| 282 | "tk": "gist_heat",\ |
|---|
| 283 | "QH2O": "PuBu",\ |
|---|
| 284 | "USTM": "YlOrRd",\ |
|---|
| 285 | #"RdPu",\ |
|---|
| 286 | } |
|---|
| 287 | if whichone not in whichcolorb: |
|---|
| 288 | whichone = "def" |
|---|
| 289 | return whichcolorb[whichone] |
|---|
| 290 | |
|---|
| 291 | def definecolorvec (whichone="def"): |
|---|
| 292 | whichcolor = { \ |
|---|
| 293 | "def": "black",\ |
|---|
| 294 | "vis": "yellow",\ |
|---|
| 295 | "vishires": "yellow",\ |
|---|
| 296 | "molabw": "yellow",\ |
|---|
| 297 | "mola": "black",\ |
|---|
| 298 | "gist_heat": "white",\ |
|---|
| 299 | "hot": "tk",\ |
|---|
| 300 | "gist_rainbow": "black",\ |
|---|
| 301 | "spectral": "black",\ |
|---|
| 302 | "gray": "red",\ |
|---|
| 303 | "PuBu": "black",\ |
|---|
| 304 | } |
|---|
| 305 | if whichone not in whichcolor: |
|---|
| 306 | whichone = "def" |
|---|
| 307 | return whichcolor[whichone] |
|---|
| 308 | |
|---|
| 309 | def marsmap (whichone="vishires"): |
|---|
| 310 | whichlink = { \ |
|---|
| 311 | "vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ |
|---|
| 312 | "vishires": "http://dl.dropbox.com/u/11078310/MarsMap_2500x1250.jpg",\ |
|---|
| 313 | "mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ |
|---|
| 314 | "molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ |
|---|
| 315 | } |
|---|
| 316 | if whichone not in whichlink: |
|---|
| 317 | print "marsmap: choice not defined... you'll get the default one... " |
|---|
| 318 | whichone = "vishires" |
|---|
| 319 | return whichlink[whichone] |
|---|
| 320 | |
|---|
| 321 | def earthmap (whichone): |
|---|
| 322 | if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg" |
|---|
| 323 | elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg" |
|---|
| 324 | elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg" |
|---|
| 325 | return whichlink |
|---|
| 326 | |
|---|
