[252] | 1 | def errormess(text,printvar=None): |
---|
[233] | 2 | print text |
---|
[252] | 3 | if printvar: print printvar |
---|
[233] | 4 | exit() |
---|
| 5 | return |
---|
| 6 | |
---|
[252] | 7 | def getname(var=False,winds=False,anomaly=False): |
---|
| 8 | if var and winds: basename = var + '_UV' |
---|
| 9 | elif var: basename = var |
---|
| 10 | elif winds: basename = 'UV' |
---|
| 11 | else: errormess("please set at least winds or var",printvar=nc.variables) |
---|
| 12 | if anomaly: basename = 'd' + basename |
---|
| 13 | return basename |
---|
| 14 | |
---|
| 15 | def localtime(utc,lon): |
---|
| 16 | ltst = utc + lon / 15. |
---|
| 17 | ltst = int (ltst * 10) / 10. |
---|
| 18 | ltst = ltst % 24 |
---|
| 19 | return ltst |
---|
| 20 | |
---|
[233] | 21 | def whatkindfile (nc): |
---|
| 22 | if 'controle' in nc.variables: typefile = 'gcm' |
---|
[254] | 23 | elif 'phisinit' in nc.variables: typefile = 'gcmex' |
---|
[233] | 24 | elif 'vert' in nc.variables: typefile = 'mesoapi' |
---|
| 25 | elif 'U' in nc.variables: typefile = 'meso' |
---|
| 26 | elif 'HGT_M' in nc.variables: typefile = 'geo' |
---|
| 27 | else: errormess("whatkindfile: typefile not supported.") |
---|
| 28 | return typefile |
---|
| 29 | |
---|
| 30 | def getfield (nc,var): |
---|
| 31 | ## this allows to get much faster (than simply referring to nc.variables[var]) |
---|
| 32 | dimension = len(nc.variables[var].dimensions) |
---|
| 33 | if dimension == 2: field = nc.variables[var][:,:] |
---|
| 34 | elif dimension == 3: field = nc.variables[var][:,:,:] |
---|
| 35 | elif dimension == 4: field = nc.variables[var][:,:,:,:] |
---|
| 36 | return field |
---|
| 37 | |
---|
| 38 | def reducefield (input,d4=None,d3=None,d2=None,d1=None): |
---|
[252] | 39 | ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x" |
---|
[233] | 40 | ### it would be actually better to name d4 d3 d2 d1 as t z y x |
---|
| 41 | import numpy as np |
---|
| 42 | dimension = np.array(input).ndim |
---|
| 43 | shape = np.array(input).shape |
---|
| 44 | print 'dim,shape: ',dimension,shape |
---|
| 45 | output = input |
---|
| 46 | error = False |
---|
| 47 | if dimension == 2: |
---|
| 48 | if d2 >= shape[0]: error = True |
---|
| 49 | elif d1 >= shape[1]: error = True |
---|
| 50 | elif d1 is not None and d2 is not None: output = input[d2,d1] |
---|
| 51 | elif d1 is not None: output = input[:,d1] |
---|
| 52 | elif d2 is not None: output = input[d2,:] |
---|
| 53 | elif dimension == 3: |
---|
| 54 | if d4 >= shape[0]: error = True |
---|
| 55 | elif d2 >= shape[1]: error = True |
---|
| 56 | elif d1 >= shape[2]: error = True |
---|
| 57 | elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,d2,d1] |
---|
| 58 | elif d4 is not None and d2 is not None: output = input[d4,d2,:] |
---|
| 59 | elif d4 is not None and d1 is not None: output = input[d4,:,d1] |
---|
| 60 | elif d2 is not None and d1 is not None: output = input[:,d2,d1] |
---|
| 61 | elif d1 is not None: output = input[:,:,d1] |
---|
| 62 | elif d2 is not None: output = input[:,d2,:] |
---|
| 63 | elif d4 is not None: output = input[d4,:,:] |
---|
| 64 | elif dimension == 4: |
---|
| 65 | if d4 >= shape[0]: error = True |
---|
| 66 | elif d3 >= shape[1]: error = True |
---|
| 67 | elif d2 >= shape[2]: error = True |
---|
| 68 | elif d1 >= shape[3]: error = True |
---|
| 69 | elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: output = input[d4,d3,d2,d1] |
---|
| 70 | elif d4 is not None and d3 is not None and d2 is not None: output = input[d4,d3,d2,:] |
---|
| 71 | elif d4 is not None and d3 is not None and d1 is not None: output = input[d4,d3,:,d1] |
---|
| 72 | elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,:,d2,d1] |
---|
| 73 | elif d3 is not None and d2 is not None and d1 is not None: output = input[:,d3,d2,d1] |
---|
| 74 | elif d4 is not None and d3 is not None: output = input[d4,d3,:,:] |
---|
| 75 | elif d4 is not None and d2 is not None: output = input[d4,:,d2,:] |
---|
| 76 | elif d4 is not None and d1 is not None: output = input[d4,:,:,d1] |
---|
| 77 | elif d3 is not None and d2 is not None: output = input[:,d3,d2,:] |
---|
| 78 | elif d3 is not None and d1 is not None: output = input[:,d3,:,d1] |
---|
| 79 | elif d2 is not None and d1 is not None: output = input[:,:,d2,d1] |
---|
| 80 | elif d1 is not None: output = input[:,:,:,d1] |
---|
| 81 | elif d2 is not None: output = input[:,:,d2,:] |
---|
| 82 | elif d3 is not None: output = input[:,d3,:,:] |
---|
| 83 | elif d4 is not None: output = input[d4,:,:,:] |
---|
| 84 | dimension = np.array(output).ndim |
---|
| 85 | shape = np.array(output).shape |
---|
| 86 | print 'dim,shape: ',dimension,shape |
---|
| 87 | return output, error |
---|
| 88 | |
---|
| 89 | def definesubplot ( numplot, fig ): |
---|
| 90 | from matplotlib.pyplot import rcParams |
---|
| 91 | rcParams['font.size'] = 12. ## default (important for multiple calls) |
---|
| 92 | if numplot == 4: |
---|
| 93 | sub = 221 |
---|
| 94 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 95 | rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) |
---|
| 96 | elif numplot == 2: |
---|
| 97 | sub = 121 |
---|
| 98 | fig.subplots_adjust(wspace = 0.35) |
---|
| 99 | rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. ) |
---|
| 100 | elif numplot == 3: |
---|
| 101 | sub = 131 |
---|
| 102 | fig.subplots_adjust(wspace = 0.5) |
---|
| 103 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
| 104 | elif numplot == 6: |
---|
| 105 | sub = 231 |
---|
| 106 | fig.subplots_adjust(wspace = 0.4, hspace = 0.0) |
---|
| 107 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
| 108 | elif numplot == 8: |
---|
| 109 | sub = 331 #241 |
---|
| 110 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 111 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
| 112 | elif numplot == 9: |
---|
| 113 | sub = 331 |
---|
| 114 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 115 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
| 116 | elif numplot == 1: |
---|
| 117 | sub = 99999 |
---|
[241] | 118 | elif numplot <= 0: |
---|
[233] | 119 | sub = 99999 |
---|
| 120 | else: |
---|
| 121 | print "supported: 1,2,3,4,6,8,9" |
---|
| 122 | exit() |
---|
| 123 | return sub |
---|
| 124 | |
---|
| 125 | def getstralt(nc,nvert): |
---|
| 126 | typefile = whatkindfile(nc) |
---|
| 127 | if typefile is 'meso': |
---|
| 128 | stralt = "_lvl" + str(nvert) |
---|
| 129 | elif typefile is 'mesoapi': |
---|
| 130 | zelevel = int(nc.variables['vert'][nvert]) |
---|
| 131 | if abs(zelevel) < 10000.: strheight=str(zelevel)+"m" |
---|
| 132 | else: strheight=str(int(zelevel/1000.))+"km" |
---|
| 133 | if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR" |
---|
| 134 | elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS" |
---|
| 135 | elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight |
---|
| 136 | elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" |
---|
| 137 | else: stralt = "" |
---|
| 138 | else: |
---|
| 139 | stralt = "" |
---|
| 140 | return stralt |
---|
| 141 | |
---|
[195] | 142 | def getlschar ( namefile ): |
---|
| 143 | from netCDF4 import Dataset |
---|
| 144 | from timestuff import sol2ls |
---|
[233] | 145 | from numpy import array |
---|
[195] | 146 | nc = Dataset(namefile) |
---|
[237] | 147 | zetime = None |
---|
[233] | 148 | if 'Times' in nc.variables: |
---|
| 149 | zetime = nc.variables['Times'][0] |
---|
| 150 | shape = array(nc.variables['Times']).shape |
---|
| 151 | if shape[0] < 2: zetime = None |
---|
| 152 | if zetime is not None \ |
---|
[225] | 153 | and 'vert' not in nc.variables: |
---|
[233] | 154 | #### strangely enough this does not work for api or ncrcat results! |
---|
[195] | 155 | zetimestart = getattr(nc, 'START_DATE') |
---|
| 156 | zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) |
---|
| 157 | if zeday < 0: lschar="" ## might have crossed a month... fix soon |
---|
[241] | 158 | else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. ) |
---|
[197] | 159 | ### |
---|
| 160 | zetime2 = nc.variables['Times'][1] |
---|
| 161 | one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37. |
---|
| 162 | next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. |
---|
| 163 | zehour = one |
---|
| 164 | zehourin = abs ( next - one ) |
---|
[195] | 165 | else: |
---|
| 166 | lschar="" |
---|
[197] | 167 | zehour = 0 |
---|
| 168 | zehourin = 1 |
---|
| 169 | return lschar, zehour, zehourin |
---|
[195] | 170 | |
---|
[202] | 171 | def getprefix (nc): |
---|
| 172 | prefix = 'LMD_MMM_' |
---|
| 173 | prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_' |
---|
| 174 | prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_' |
---|
| 175 | return prefix |
---|
| 176 | |
---|
[184] | 177 | def getproj (nc): |
---|
[233] | 178 | typefile = whatkindfile(nc) |
---|
| 179 | if typefile in ['mesoapi','meso','geo']: |
---|
| 180 | ### (il faudrait passer CEN_LON dans la projection ?) |
---|
| 181 | map_proj = getattr(nc, 'MAP_PROJ') |
---|
| 182 | cen_lat = getattr(nc, 'CEN_LAT') |
---|
| 183 | if map_proj == 2: |
---|
| 184 | if cen_lat > 10.: |
---|
| 185 | proj="npstere" |
---|
| 186 | print "NP stereographic polar domain" |
---|
| 187 | else: |
---|
| 188 | proj="spstere" |
---|
| 189 | print "SP stereographic polar domain" |
---|
| 190 | elif map_proj == 1: |
---|
| 191 | print "lambert projection domain" |
---|
| 192 | proj="lcc" |
---|
| 193 | elif map_proj == 3: |
---|
| 194 | print "mercator projection" |
---|
| 195 | proj="merc" |
---|
| 196 | else: |
---|
| 197 | proj="merc" |
---|
[252] | 198 | elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?) |
---|
[233] | 199 | else: proj="ortho" |
---|
[184] | 200 | return proj |
---|
| 201 | |
---|
[180] | 202 | def ptitle (name): |
---|
| 203 | from matplotlib.pyplot import title |
---|
| 204 | title(name) |
---|
| 205 | print name |
---|
| 206 | |
---|
[252] | 207 | def polarinterv (lon2d,lat2d): |
---|
| 208 | import numpy as np |
---|
| 209 | wlon = [np.min(lon2d),np.max(lon2d)] |
---|
| 210 | ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole |
---|
| 211 | wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])] |
---|
| 212 | return [wlon,wlat] |
---|
| 213 | |
---|
[180] | 214 | def simplinterv (lon2d,lat2d): |
---|
| 215 | import numpy as np |
---|
| 216 | return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]] |
---|
| 217 | |
---|
[184] | 218 | def wrfinterv (lon2d,lat2d): |
---|
| 219 | nx = len(lon2d[0,:])-1 |
---|
| 220 | ny = len(lon2d[:,0])-1 |
---|
[225] | 221 | lon1 = lon2d[0,0] |
---|
| 222 | lon2 = lon2d[nx,ny] |
---|
| 223 | lat1 = lat2d[0,0] |
---|
| 224 | lat2 = lat2d[nx,ny] |
---|
[233] | 225 | if abs(0.5*(lat1+lat2)) > 60.: wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 |
---|
| 226 | else: wider = 0. |
---|
| 227 | if lon1 < lon2: wlon = [lon1, lon2 + wider] |
---|
[225] | 228 | else: wlon = [lon2, lon1 + wider] |
---|
| 229 | if lat1 < lat2: wlat = [lat1, lat2] |
---|
| 230 | else: wlat = [lat2, lat1] |
---|
| 231 | return [wlon,wlat] |
---|
[184] | 232 | |
---|
[240] | 233 | def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False): |
---|
[180] | 234 | import matplotlib.pyplot as plt |
---|
[240] | 235 | from os import system |
---|
| 236 | addstr = "" |
---|
| 237 | if res is not None: |
---|
| 238 | res = int(res) |
---|
| 239 | addstr = "_"+str(res) |
---|
| 240 | name = filename+addstr+"."+ext |
---|
[186] | 241 | if folder != '': name = folder+'/'+name |
---|
[180] | 242 | plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value) |
---|
[240] | 243 | if disp: display(name) |
---|
| 244 | if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name) |
---|
| 245 | if erase: system("mv "+name+" to_be_erased") |
---|
[180] | 246 | return |
---|
| 247 | |
---|
[240] | 248 | def dumpbdy (field,n,stag=None): |
---|
[184] | 249 | nx = len(field[0,:])-1 |
---|
| 250 | ny = len(field[:,0])-1 |
---|
[233] | 251 | if stag == 'U': nx = nx-1 |
---|
| 252 | if stag == 'V': ny = ny-1 |
---|
[240] | 253 | return field[n:ny-n,n:nx-n] |
---|
[180] | 254 | |
---|
[233] | 255 | def getcoorddef ( nc ): |
---|
| 256 | ## getcoord2d for predefined types |
---|
| 257 | typefile = whatkindfile(nc) |
---|
| 258 | if typefile in ['mesoapi','meso']: |
---|
| 259 | [lon2d,lat2d] = getcoord2d(nc) |
---|
[240] | 260 | lon2d = dumpbdy(lon2d,6) |
---|
| 261 | lat2d = dumpbdy(lat2d,6) |
---|
[254] | 262 | elif typefile in ['gcm','gcmex']: |
---|
[233] | 263 | [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) |
---|
| 264 | elif typefile in ['geo']: |
---|
| 265 | [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') |
---|
| 266 | return lon2d,lat2d |
---|
| 267 | |
---|
[184] | 268 | def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): |
---|
| 269 | import numpy as np |
---|
| 270 | if is1d: |
---|
| 271 | lat = nc.variables[nlat][:] |
---|
| 272 | lon = nc.variables[nlon][:] |
---|
| 273 | [lon2d,lat2d] = np.meshgrid(lon,lat) |
---|
| 274 | else: |
---|
| 275 | lat = nc.variables[nlat][0,:,:] |
---|
| 276 | lon = nc.variables[nlon][0,:,:] |
---|
| 277 | [lon2d,lat2d] = [lon,lat] |
---|
| 278 | return lon2d,lat2d |
---|
| 279 | |
---|
[180] | 280 | def smooth (field, coeff): |
---|
| 281 | ## actually blur_image could work with different coeff on x and y |
---|
| 282 | if coeff > 1: result = blur_image(field,int(coeff)) |
---|
| 283 | else: result = field |
---|
| 284 | return result |
---|
| 285 | |
---|
| 286 | def gauss_kern(size, sizey=None): |
---|
| 287 | import numpy as np |
---|
| 288 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
---|
| 289 | # Returns a normalized 2D gauss kernel array for convolutions |
---|
| 290 | size = int(size) |
---|
| 291 | if not sizey: |
---|
| 292 | sizey = size |
---|
| 293 | else: |
---|
| 294 | sizey = int(sizey) |
---|
| 295 | x, y = np.mgrid[-size:size+1, -sizey:sizey+1] |
---|
| 296 | g = np.exp(-(x**2/float(size)+y**2/float(sizey))) |
---|
| 297 | return g / g.sum() |
---|
| 298 | |
---|
| 299 | def blur_image(im, n, ny=None) : |
---|
| 300 | from scipy.signal import convolve |
---|
| 301 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
---|
| 302 | # blurs the image by convolving with a gaussian kernel of typical size n. |
---|
| 303 | # The optional keyword argument ny allows for a different size in the y direction. |
---|
| 304 | g = gauss_kern(n, sizey=ny) |
---|
| 305 | improc = convolve(im, g, mode='same') |
---|
| 306 | return improc |
---|
| 307 | |
---|
[233] | 308 | def getwinddef (nc): |
---|
| 309 | ## getwinds for predefined types |
---|
| 310 | typefile = whatkindfile(nc) |
---|
| 311 | ### |
---|
| 312 | if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm'] |
---|
| 313 | elif typefile is 'gcm': [uchar,vchar] = ['u','v'] |
---|
| 314 | elif typefile is 'meso': [uchar,vchar] = ['U','V'] |
---|
| 315 | else: [uchar,vchar] = ['not found','not found'] |
---|
| 316 | ### |
---|
| 317 | if typefile in ['meso']: metwind = False ## geometrical (wrt grid) |
---|
| 318 | else: metwind = True ## meteorological (zon/mer) |
---|
| 319 | if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)" |
---|
| 320 | ### |
---|
| 321 | return uchar,vchar,metwind |
---|
[202] | 322 | |
---|
[184] | 323 | def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True): |
---|
| 324 | ## scale regle la reference du vecteur |
---|
| 325 | ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs. |
---|
| 326 | import matplotlib.pyplot as plt |
---|
| 327 | import numpy as np |
---|
[187] | 328 | posx = np.min(x) - np.std(x) / 10. |
---|
| 329 | posy = np.min(y) - np.std(y) / 10. |
---|
[184] | 330 | u = smooth(u,csmooth) |
---|
| 331 | v = smooth(v,csmooth) |
---|
[188] | 332 | widthvec = 0.003 #0.005 #0.003 |
---|
[184] | 333 | q = plt.quiver( x[::stride,::stride],\ |
---|
| 334 | y[::stride,::stride],\ |
---|
| 335 | u[::stride,::stride],\ |
---|
| 336 | v[::stride,::stride],\ |
---|
[228] | 337 | angles='xy',color=color,pivot='middle',\ |
---|
[184] | 338 | scale=factor,width=widthvec ) |
---|
[202] | 339 | if color in ['white','yellow']: kcolor='black' |
---|
| 340 | else: kcolor=color |
---|
[184] | 341 | if key: p = plt.quiverkey(q,posx,posy,scale,\ |
---|
[194] | 342 | str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03) |
---|
[184] | 343 | return |
---|
[180] | 344 | |
---|
| 345 | def display (name): |
---|
[184] | 346 | from os import system |
---|
| 347 | system("display "+name+" > /dev/null 2> /dev/null &") |
---|
| 348 | return name |
---|
[180] | 349 | |
---|
| 350 | def findstep (wlon): |
---|
[184] | 351 | steplon = int((wlon[1]-wlon[0])/4.) #3 |
---|
| 352 | step = 120. |
---|
| 353 | while step > steplon and step > 15. : step = step / 2. |
---|
| 354 | if step <= 15.: |
---|
| 355 | while step > steplon and step > 5. : step = step - 5. |
---|
| 356 | if step <= 5.: |
---|
| 357 | while step > steplon and step > 1. : step = step - 1. |
---|
| 358 | if step <= 1.: |
---|
| 359 | step = 1. |
---|
[180] | 360 | return step |
---|
| 361 | |
---|
[233] | 362 | def define_proj (char,wlon,wlat,back=None): |
---|
[180] | 363 | from mpl_toolkits.basemap import Basemap |
---|
| 364 | import numpy as np |
---|
| 365 | import matplotlib as mpl |
---|
[240] | 366 | from mymath import max |
---|
[180] | 367 | meanlon = 0.5*(wlon[0]+wlon[1]) |
---|
| 368 | meanlat = 0.5*(wlat[0]+wlat[1]) |
---|
[184] | 369 | if wlat[0] >= 80.: blat = 40. |
---|
[247] | 370 | elif wlat[1] <= -80.: blat = -40. |
---|
| 371 | elif wlat[1] >= 0.: blat = wlat[0] |
---|
| 372 | elif wlat[0] <= 0.: blat = wlat[1] |
---|
[252] | 373 | print "blat ", blat |
---|
[207] | 374 | h = 50. ## en km |
---|
[202] | 375 | radius = 3397200. |
---|
[184] | 376 | if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\ |
---|
[180] | 377 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
[184] | 378 | elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon) |
---|
| 379 | elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat) |
---|
| 380 | elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\ |
---|
| 381 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
| 382 | elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.) |
---|
| 383 | elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.) |
---|
[207] | 384 | elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon) |
---|
| 385 | elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\ |
---|
| 386 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
[184] | 387 | elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.) |
---|
| 388 | elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\ |
---|
| 389 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
| 390 | fontsizemer = int(mpl.rcParams['font.size']*3./4.) |
---|
[207] | 391 | if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon) |
---|
| 392 | else: step = 10. |
---|
[238] | 393 | steplon = step*2. |
---|
| 394 | #if back in ["geolocal"]: |
---|
| 395 | # step = np.min([5.,step]) |
---|
| 396 | # steplon = step |
---|
[225] | 397 | print step |
---|
[238] | 398 | m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer) |
---|
[180] | 399 | m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer) |
---|
[233] | 400 | if back: m.warpimage(marsmap(back),scale=0.75) |
---|
| 401 | #if not back: |
---|
| 402 | # if not var: back = "mola" ## if no var: draw mola |
---|
| 403 | # elif typefile in ['mesoapi','meso','geo'] \ |
---|
| 404 | # and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw |
---|
| 405 | # else: pass ## else: draw None |
---|
[180] | 406 | return m |
---|
| 407 | |
---|
[232] | 408 | #### test temporaire |
---|
| 409 | def putpoints (map,plot): |
---|
| 410 | #### from http://www.scipy.org/Cookbook/Matplotlib/Maps |
---|
| 411 | # lat/lon coordinates of five cities. |
---|
| 412 | lats = [18.4] |
---|
| 413 | lons = [-134.0] |
---|
| 414 | points=['Olympus Mons'] |
---|
| 415 | # compute the native map projection coordinates for cities. |
---|
| 416 | x,y = map(lons,lats) |
---|
| 417 | # plot filled circles at the locations of the cities. |
---|
| 418 | map.plot(x,y,'bo') |
---|
| 419 | # plot the names of those five cities. |
---|
| 420 | wherept = 0 #1000 #50000 |
---|
| 421 | for name,xpt,ypt in zip(points,x,y): |
---|
| 422 | plot.text(xpt+wherept,ypt+wherept,name) |
---|
| 423 | ## le nom ne s'affiche pas... |
---|
| 424 | return |
---|
| 425 | |
---|
[233] | 426 | def calculate_bounds(field,vmin=None,vmax=None): |
---|
| 427 | import numpy as np |
---|
| 428 | from mymath import max,min,mean |
---|
| 429 | ind = np.where(field < 9e+35) |
---|
| 430 | fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple |
---|
| 431 | ### |
---|
| 432 | dev = np.std(fieldcalc)*3.0 |
---|
| 433 | ### |
---|
| 434 | if vmin is None: |
---|
| 435 | zevmin = mean(fieldcalc) - dev |
---|
| 436 | else: zevmin = vmin |
---|
| 437 | ### |
---|
| 438 | if vmax is None: zevmax = mean(fieldcalc) + dev |
---|
| 439 | else: zevmax = vmax |
---|
| 440 | if vmin == vmax: |
---|
| 441 | zevmin = mean(fieldcalc) - dev ### for continuity |
---|
| 442 | zevmax = mean(fieldcalc) + dev ### for continuity |
---|
| 443 | ### |
---|
| 444 | if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0. |
---|
| 445 | print "field ", min(fieldcalc), max(fieldcalc) |
---|
| 446 | print "bounds ", zevmin, zevmax |
---|
| 447 | return zevmin, zevmax |
---|
[232] | 448 | |
---|
[233] | 449 | def bounds(what_I_plot,zevmin,zevmax): |
---|
[247] | 450 | from mymath import max,min,mean |
---|
[233] | 451 | ### might be convenient to add the missing value in arguments |
---|
| 452 | what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) |
---|
[247] | 453 | print "new min ", min(what_I_plot) |
---|
[233] | 454 | what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 |
---|
| 455 | what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) |
---|
[247] | 456 | print "new max ", max(what_I_plot) |
---|
[233] | 457 | return what_I_plot |
---|
| 458 | |
---|
[241] | 459 | def nolow(what_I_plot): |
---|
| 460 | from mymath import max,min |
---|
| 461 | lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot))) |
---|
| 462 | print "on vire en dessous de ", lim |
---|
| 463 | what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 |
---|
| 464 | return what_I_plot |
---|
| 465 | |
---|
[233] | 466 | def zoomset (wlon,wlat,zoom): |
---|
| 467 | dlon = abs(wlon[1]-wlon[0])/2. |
---|
| 468 | dlat = abs(wlat[1]-wlat[0])/2. |
---|
| 469 | [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\ |
---|
| 470 | [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ] |
---|
| 471 | print "zoom %",zoom,wlon,wlat |
---|
| 472 | return wlon,wlat |
---|
| 473 | |
---|
[201] | 474 | def fmtvar (whichvar="def"): |
---|
[204] | 475 | fmtvar = { \ |
---|
| 476 | "tk": "%.0f",\ |
---|
| 477 | "tpot": "%.0f",\ |
---|
| 478 | "def": "%.1e",\ |
---|
| 479 | "PTOT": "%.0f",\ |
---|
| 480 | "HGT": "%.1e",\ |
---|
| 481 | "USTM": "%.2f",\ |
---|
[225] | 482 | "HFX": "%.0f",\ |
---|
[232] | 483 | "ICETOT": "%.1e",\ |
---|
[237] | 484 | "TAU_ICE": "%.2f",\ |
---|
[252] | 485 | "VMR_ICE": "%.1e",\ |
---|
| 486 | "MTOT": "%.0f",\ |
---|
[240] | 487 | "anomaly": "%.1f",\ |
---|
[241] | 488 | "W": "%.1f",\ |
---|
[287] | 489 | "WMAX_TH": "%.1f",\ |
---|
| 490 | "QSURFICE": "%.0f",\ |
---|
[204] | 491 | } |
---|
| 492 | if whichvar not in fmtvar: |
---|
| 493 | whichvar = "def" |
---|
| 494 | return fmtvar[whichvar] |
---|
[201] | 495 | |
---|
[233] | 496 | #################################################################################################################### |
---|
| 497 | ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png |
---|
[202] | 498 | def defcolorb (whichone="def"): |
---|
[204] | 499 | whichcolorb = { \ |
---|
| 500 | "def": "spectral",\ |
---|
| 501 | "HGT": "spectral",\ |
---|
| 502 | "tk": "gist_heat",\ |
---|
| 503 | "QH2O": "PuBu",\ |
---|
| 504 | "USTM": "YlOrRd",\ |
---|
[225] | 505 | "HFX": "RdYlBu",\ |
---|
[232] | 506 | "ICETOT": "YlGnBu",\ |
---|
[252] | 507 | "MTOT": "PuBu",\ |
---|
[238] | 508 | "TAU_ICE": "Blues",\ |
---|
[252] | 509 | "VMR_ICE": "Blues",\ |
---|
[241] | 510 | "W": "jet",\ |
---|
[287] | 511 | "WMAX_TH": "spectral",\ |
---|
[240] | 512 | "anomaly": "RdBu_r",\ |
---|
[287] | 513 | "QSURFICE": "hot_r",\ |
---|
[204] | 514 | } |
---|
[241] | 515 | #W --> spectral ou jet |
---|
[240] | 516 | #spectral BrBG RdBu_r |
---|
[241] | 517 | print "predefined colorbars" |
---|
[204] | 518 | if whichone not in whichcolorb: |
---|
| 519 | whichone = "def" |
---|
| 520 | return whichcolorb[whichone] |
---|
[202] | 521 | |
---|
| 522 | def definecolorvec (whichone="def"): |
---|
| 523 | whichcolor = { \ |
---|
| 524 | "def": "black",\ |
---|
| 525 | "vis": "yellow",\ |
---|
| 526 | "vishires": "yellow",\ |
---|
| 527 | "molabw": "yellow",\ |
---|
| 528 | "mola": "black",\ |
---|
| 529 | "gist_heat": "white",\ |
---|
| 530 | "hot": "tk",\ |
---|
| 531 | "gist_rainbow": "black",\ |
---|
| 532 | "spectral": "black",\ |
---|
| 533 | "gray": "red",\ |
---|
| 534 | "PuBu": "black",\ |
---|
| 535 | } |
---|
| 536 | if whichone not in whichcolor: |
---|
| 537 | whichone = "def" |
---|
| 538 | return whichcolor[whichone] |
---|
| 539 | |
---|
[180] | 540 | def marsmap (whichone="vishires"): |
---|
[233] | 541 | from os import uname |
---|
| 542 | mymachine = uname()[1] |
---|
| 543 | ### not sure about speed-up with this method... looks the same |
---|
| 544 | if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/" |
---|
| 545 | else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/" |
---|
[180] | 546 | whichlink = { \ |
---|
[233] | 547 | #"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ |
---|
| 548 | #"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\ |
---|
| 549 | #"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\ |
---|
| 550 | #"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ |
---|
| 551 | #"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ |
---|
| 552 | "vis": domain+"mar0kuu2.jpg",\ |
---|
| 553 | "vishires": domain+"MarsMap_2500x1250.jpg",\ |
---|
| 554 | "geolocal": domain+"geolocal.jpg",\ |
---|
| 555 | "mola": domain+"mars-mola-2k.jpg",\ |
---|
| 556 | "molabw": domain+"MarsElevation_2500x1250.jpg",\ |
---|
[238] | 557 | "clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\ |
---|
| 558 | "jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\ |
---|
| 559 | "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\ |
---|
[273] | 560 | "bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\ |
---|
| 561 | "contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\ |
---|
| 562 | "nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\ |
---|
| 563 | "blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\ |
---|
[180] | 564 | } |
---|
[238] | 565 | ### see http://www.mmedia.is/~bjj/planetary_maps.html |
---|
[180] | 566 | if whichone not in whichlink: |
---|
| 567 | print "marsmap: choice not defined... you'll get the default one... " |
---|
| 568 | whichone = "vishires" |
---|
| 569 | return whichlink[whichone] |
---|
| 570 | |
---|
[273] | 571 | #def earthmap (whichone): |
---|
| 572 | # if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg" |
---|
| 573 | # elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg" |
---|
| 574 | # elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg" |
---|
| 575 | # return whichlink |
---|
[180] | 576 | |
---|
[241] | 577 | def latinterv (area="Whole"): |
---|
| 578 | list = { \ |
---|
| 579 | "Europe": [[ 20., 80.],[- 50., 50.]],\ |
---|
| 580 | "Central_America": [[-10., 40.],[ 230., 300.]],\ |
---|
| 581 | "Africa": [[-20., 50.],[- 50., 50.]],\ |
---|
[273] | 582 | "Whole": [[-90., 90.],[-180., 180.]],\ |
---|
| 583 | "Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\ |
---|
| 584 | "Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\ |
---|
[241] | 585 | "Tharsis": [[-30., 60.],[-170.,- 10.]],\ |
---|
| 586 | "Whole_No_High": [[-60., 60.],[-180., 180.]],\ |
---|
| 587 | "Chryse": [[-60., 60.],[- 60., 60.]],\ |
---|
| 588 | "North_Pole": [[ 50., 90.],[-180., 180.]],\ |
---|
| 589 | "Close_North_Pole": [[ 75., 90.],[-180., 180.]],\ |
---|
| 590 | "Far_South_Pole": [[-90.,-40.],[-180., 180.]],\ |
---|
| 591 | "South_Pole": [[-90.,-50.],[-180., 180.]],\ |
---|
| 592 | "Close_South_Pole": [[-90.,-75.],[-180., 180.]],\ |
---|
| 593 | } |
---|
| 594 | if area not in list: area = "Whole" |
---|
| 595 | [olat,olon] = list[area] |
---|
| 596 | return olon,olat |
---|
| 597 | |
---|