| 1 | ## Author: AS |
|---|
| 2 | def errormess(text,printvar=None): |
|---|
| 3 | print text |
|---|
| 4 | if printvar: print printvar |
|---|
| 5 | exit() |
|---|
| 6 | return |
|---|
| 7 | |
|---|
| 8 | ## Author: AS |
|---|
| 9 | def adjust_length (tab, zelen): |
|---|
| 10 | from numpy import ones |
|---|
| 11 | if tab is None: |
|---|
| 12 | outtab = ones(zelen) * -999999 |
|---|
| 13 | else: |
|---|
| 14 | if zelen != len(tab): |
|---|
| 15 | print "not enough or too much values... setting same values all variables" |
|---|
| 16 | outtab = ones(zelen) * tab[0] |
|---|
| 17 | else: |
|---|
| 18 | outtab = tab |
|---|
| 19 | return outtab |
|---|
| 20 | |
|---|
| 21 | ## Author: AS |
|---|
| 22 | def getname(var=False,winds=False,anomaly=False): |
|---|
| 23 | if var and winds: basename = var + '_UV' |
|---|
| 24 | elif var: basename = var |
|---|
| 25 | elif winds: basename = 'UV' |
|---|
| 26 | else: errormess("please set at least winds or var",printvar=nc.variables) |
|---|
| 27 | if anomaly: basename = 'd' + basename |
|---|
| 28 | return basename |
|---|
| 29 | |
|---|
| 30 | ## Author: AS |
|---|
| 31 | def localtime(utc,lon): |
|---|
| 32 | ltst = utc + lon / 15. |
|---|
| 33 | ltst = int (ltst * 10) / 10. |
|---|
| 34 | ltst = ltst % 24 |
|---|
| 35 | return ltst |
|---|
| 36 | |
|---|
| 37 | ## Author: AS |
|---|
| 38 | def whatkindfile (nc): |
|---|
| 39 | if 'controle' in nc.variables: typefile = 'gcm' |
|---|
| 40 | elif 'phisinit' in nc.variables: typefile = 'gcm' |
|---|
| 41 | elif 'vert' in nc.variables: typefile = 'mesoapi' |
|---|
| 42 | elif 'U' in nc.variables: typefile = 'meso' |
|---|
| 43 | elif 'HGT_M' in nc.variables: typefile = 'geo' |
|---|
| 44 | #else: errormess("whatkindfile: typefile not supported.") |
|---|
| 45 | else: typefile = 'gcm' # for lslin-ed files from gcm |
|---|
| 46 | return typefile |
|---|
| 47 | |
|---|
| 48 | ## Author: AS |
|---|
| 49 | def getfield (nc,var): |
|---|
| 50 | ## this allows to get much faster (than simply referring to nc.variables[var]) |
|---|
| 51 | dimension = len(nc.variables[var].dimensions) |
|---|
| 52 | print " Opening variable with", dimension, "dimensions ..." |
|---|
| 53 | if dimension == 2: field = nc.variables[var][:,:] |
|---|
| 54 | elif dimension == 3: field = nc.variables[var][:,:,:] |
|---|
| 55 | elif dimension == 4: field = nc.variables[var][:,:,:,:] |
|---|
| 56 | return field |
|---|
| 57 | |
|---|
| 58 | ## Author: AS + TN |
|---|
| 59 | def reducefield (input,d4=None,d3=None,d2=None,d1=None): |
|---|
| 60 | ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x" |
|---|
| 61 | ### it would be actually better to name d4 d3 d2 d1 as t z y x |
|---|
| 62 | import numpy as np |
|---|
| 63 | from mymath import max,mean |
|---|
| 64 | dimension = np.array(input).ndim |
|---|
| 65 | shape = np.array(input).shape |
|---|
| 66 | #print 'd1,d2,d3,d4: ',d1,d2,d3,d4 |
|---|
| 67 | print 'dim,shape: ',dimension,shape |
|---|
| 68 | output = input |
|---|
| 69 | error = False |
|---|
| 70 | #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays |
|---|
| 71 | if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4] |
|---|
| 72 | if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3] |
|---|
| 73 | if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2] |
|---|
| 74 | if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1] |
|---|
| 75 | ### now the main part |
|---|
| 76 | if dimension == 2: |
|---|
| 77 | if d2 >= shape[0]: error = True |
|---|
| 78 | elif d1 >= shape[1]: error = True |
|---|
| 79 | elif d1 is not None and d2 is not None: output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0) |
|---|
| 80 | elif d1 is not None: output = mean(input[:,d1],axis=1) |
|---|
| 81 | elif d2 is not None: output = mean(input[d2,:],axis=0) |
|---|
| 82 | elif dimension == 3: |
|---|
| 83 | if max(d4) >= shape[0]: error = True |
|---|
| 84 | elif max(d2) >= shape[1]: error = True |
|---|
| 85 | elif max(d1) >= shape[2]: error = True |
|---|
| 86 | elif d4 is not None and d2 is not None and d1 is not None: |
|---|
| 87 | output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0) |
|---|
| 88 | elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0) |
|---|
| 89 | elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1) |
|---|
| 90 | elif d2 is not None and d1 is not None: output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1) |
|---|
| 91 | elif d1 is not None: output = mean(input[:,:,d1],axis=2) |
|---|
| 92 | elif d2 is not None: output = mean(input[:,d2,:],axis=1) |
|---|
| 93 | elif d4 is not None: output = mean(input[d4,:,:],axis=0) |
|---|
| 94 | elif dimension == 4: |
|---|
| 95 | if max(d4) >= shape[0]: error = True |
|---|
| 96 | elif max(d3) >= shape[1]: error = True |
|---|
| 97 | elif max(d2) >= shape[2]: error = True |
|---|
| 98 | elif max(d1) >= shape[3]: error = True |
|---|
| 99 | elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: |
|---|
| 100 | output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0) |
|---|
| 101 | elif d4 is not None and d3 is not None and d2 is not None: |
|---|
| 102 | output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0) |
|---|
| 103 | elif d4 is not None and d3 is not None and d1 is not None: |
|---|
| 104 | output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[:,d1],axis=1) |
|---|
| 105 | elif d4 is not None and d2 is not None and d1 is not None: |
|---|
| 106 | output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1) |
|---|
| 107 | elif d3 is not None and d2 is not None and d1 is not None: |
|---|
| 108 | output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1) |
|---|
| 109 | elif d4 is not None and d3 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0) |
|---|
| 110 | elif d4 is not None and d2 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1) |
|---|
| 111 | elif d4 is not None and d1 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,:,d1],axis=2) |
|---|
| 112 | elif d3 is not None and d2 is not None: output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1) |
|---|
| 113 | elif d3 is not None and d1 is not None: output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0) |
|---|
| 114 | elif d2 is not None and d1 is not None: output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2) |
|---|
| 115 | elif d1 is not None: output = mean(input[:,:,:,d1],axis=3) |
|---|
| 116 | elif d2 is not None: output = mean(input[:,:,d2,:],axis=2) |
|---|
| 117 | elif d3 is not None: output = mean(input[:,d3,:,:],axis=1) |
|---|
| 118 | elif d4 is not None: output = mean(input[d4,:,:,:],axis=0) |
|---|
| 119 | dimension = np.array(output).ndim |
|---|
| 120 | shape = np.array(output).shape |
|---|
| 121 | print 'dim,shape: ',dimension,shape |
|---|
| 122 | return output, error |
|---|
| 123 | |
|---|
| 124 | ## Author: AS + TN |
|---|
| 125 | def definesubplot ( numplot, fig ): |
|---|
| 126 | from matplotlib.pyplot import rcParams |
|---|
| 127 | rcParams['font.size'] = 12. ## default (important for multiple calls) |
|---|
| 128 | if numplot <= 0: |
|---|
| 129 | subv = 99999 |
|---|
| 130 | subh = 99999 |
|---|
| 131 | elif numplot == 1: |
|---|
| 132 | subv = 99999 |
|---|
| 133 | subh = 99999 |
|---|
| 134 | elif numplot == 2: |
|---|
| 135 | subv = 1 |
|---|
| 136 | subh = 2 |
|---|
| 137 | fig.subplots_adjust(wspace = 0.35) |
|---|
| 138 | rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. ) |
|---|
| 139 | elif numplot == 3: |
|---|
| 140 | subv = 3 |
|---|
| 141 | subh = 1 |
|---|
| 142 | fig.subplots_adjust(wspace = 0.5) |
|---|
| 143 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
|---|
| 144 | elif numplot == 4: |
|---|
| 145 | subv = 2 |
|---|
| 146 | subh = 2 |
|---|
| 147 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
|---|
| 148 | rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) |
|---|
| 149 | elif numplot <= 6: |
|---|
| 150 | subv = 2 |
|---|
| 151 | subh = 3 |
|---|
| 152 | fig.subplots_adjust(wspace = 0.4, hspace = 0.0) |
|---|
| 153 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
|---|
| 154 | elif numplot <= 8: |
|---|
| 155 | subv = 2 |
|---|
| 156 | subh = 4 |
|---|
| 157 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
|---|
| 158 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
|---|
| 159 | elif numplot <= 9: |
|---|
| 160 | subv = 3 |
|---|
| 161 | subh = 3 |
|---|
| 162 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
|---|
| 163 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
|---|
| 164 | elif numplot <= 12: |
|---|
| 165 | subv = 3 |
|---|
| 166 | subh = 4 |
|---|
| 167 | fig.subplots_adjust(wspace = 0.1, hspace = 0.1) |
|---|
| 168 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
|---|
| 169 | elif numplot <= 16: |
|---|
| 170 | subv = 4 |
|---|
| 171 | subh = 4 |
|---|
| 172 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
|---|
| 173 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
|---|
| 174 | else: |
|---|
| 175 | print "number of plot supported: 1 to 16" |
|---|
| 176 | exit() |
|---|
| 177 | return subv,subh |
|---|
| 178 | |
|---|
| 179 | ## Author: AS |
|---|
| 180 | def getstralt(nc,nvert): |
|---|
| 181 | typefile = whatkindfile(nc) |
|---|
| 182 | if typefile is 'meso': |
|---|
| 183 | stralt = "_lvl" + str(nvert) |
|---|
| 184 | elif typefile is 'mesoapi': |
|---|
| 185 | zelevel = int(nc.variables['vert'][nvert]) |
|---|
| 186 | if abs(zelevel) < 10000.: strheight=str(zelevel)+"m" |
|---|
| 187 | else: strheight=str(int(zelevel/1000.))+"km" |
|---|
| 188 | if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR" |
|---|
| 189 | elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS" |
|---|
| 190 | elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight |
|---|
| 191 | elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" |
|---|
| 192 | else: stralt = "" |
|---|
| 193 | else: |
|---|
| 194 | stralt = "" |
|---|
| 195 | return stralt |
|---|
| 196 | |
|---|
| 197 | ## Author: AS |
|---|
| 198 | def getlschar ( namefile ): |
|---|
| 199 | from netCDF4 import Dataset |
|---|
| 200 | from timestuff import sol2ls |
|---|
| 201 | from numpy import array |
|---|
| 202 | nc = Dataset(namefile) |
|---|
| 203 | zetime = None |
|---|
| 204 | if 'Times' in nc.variables: |
|---|
| 205 | zetime = nc.variables['Times'][0] |
|---|
| 206 | shape = array(nc.variables['Times']).shape |
|---|
| 207 | if shape[0] < 2: zetime = None |
|---|
| 208 | if zetime is not None \ |
|---|
| 209 | and 'vert' not in nc.variables: |
|---|
| 210 | #### strangely enough this does not work for api or ncrcat results! |
|---|
| 211 | zetimestart = getattr(nc, 'START_DATE') |
|---|
| 212 | zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) |
|---|
| 213 | if zeday < 0: lschar="" ## might have crossed a month... fix soon |
|---|
| 214 | else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. ) |
|---|
| 215 | ### |
|---|
| 216 | zetime2 = nc.variables['Times'][1] |
|---|
| 217 | one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37. |
|---|
| 218 | next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. |
|---|
| 219 | zehour = one |
|---|
| 220 | zehourin = abs ( next - one ) |
|---|
| 221 | else: |
|---|
| 222 | lschar="" |
|---|
| 223 | zehour = 0 |
|---|
| 224 | zehourin = 1 |
|---|
| 225 | return lschar, zehour, zehourin |
|---|
| 226 | |
|---|
| 227 | ## Author: AS |
|---|
| 228 | def getprefix (nc): |
|---|
| 229 | prefix = 'LMD_MMM_' |
|---|
| 230 | prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_' |
|---|
| 231 | prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_' |
|---|
| 232 | return prefix |
|---|
| 233 | |
|---|
| 234 | ## Author: AS |
|---|
| 235 | def getproj (nc): |
|---|
| 236 | typefile = whatkindfile(nc) |
|---|
| 237 | if typefile in ['mesoapi','meso','geo']: |
|---|
| 238 | ### (il faudrait passer CEN_LON dans la projection ?) |
|---|
| 239 | map_proj = getattr(nc, 'MAP_PROJ') |
|---|
| 240 | cen_lat = getattr(nc, 'CEN_LAT') |
|---|
| 241 | if map_proj == 2: |
|---|
| 242 | if cen_lat > 10.: |
|---|
| 243 | proj="npstere" |
|---|
| 244 | print "NP stereographic polar domain" |
|---|
| 245 | else: |
|---|
| 246 | proj="spstere" |
|---|
| 247 | print "SP stereographic polar domain" |
|---|
| 248 | elif map_proj == 1: |
|---|
| 249 | print "lambert projection domain" |
|---|
| 250 | proj="lcc" |
|---|
| 251 | elif map_proj == 3: |
|---|
| 252 | print "mercator projection" |
|---|
| 253 | proj="merc" |
|---|
| 254 | else: |
|---|
| 255 | proj="merc" |
|---|
| 256 | elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?) |
|---|
| 257 | else: proj="ortho" |
|---|
| 258 | return proj |
|---|
| 259 | |
|---|
| 260 | ## Author: AS |
|---|
| 261 | def ptitle (name): |
|---|
| 262 | from matplotlib.pyplot import title |
|---|
| 263 | title(name) |
|---|
| 264 | print name |
|---|
| 265 | |
|---|
| 266 | ## Author: AS |
|---|
| 267 | def polarinterv (lon2d,lat2d): |
|---|
| 268 | import numpy as np |
|---|
| 269 | wlon = [np.min(lon2d),np.max(lon2d)] |
|---|
| 270 | ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole |
|---|
| 271 | wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])] |
|---|
| 272 | return [wlon,wlat] |
|---|
| 273 | |
|---|
| 274 | ## Author: AS |
|---|
| 275 | def simplinterv (lon2d,lat2d): |
|---|
| 276 | import numpy as np |
|---|
| 277 | return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]] |
|---|
| 278 | |
|---|
| 279 | ## Author: AS |
|---|
| 280 | def wrfinterv (lon2d,lat2d): |
|---|
| 281 | nx = len(lon2d[0,:])-1 |
|---|
| 282 | ny = len(lon2d[:,0])-1 |
|---|
| 283 | lon1 = lon2d[0,0] |
|---|
| 284 | lon2 = lon2d[nx,ny] |
|---|
| 285 | lat1 = lat2d[0,0] |
|---|
| 286 | lat2 = lat2d[nx,ny] |
|---|
| 287 | if abs(0.5*(lat1+lat2)) > 60.: wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 |
|---|
| 288 | else: wider = 0. |
|---|
| 289 | if lon1 < lon2: wlon = [lon1, lon2 + wider] |
|---|
| 290 | else: wlon = [lon2, lon1 + wider] |
|---|
| 291 | if lat1 < lat2: wlat = [lat1, lat2] |
|---|
| 292 | else: wlat = [lat2, lat1] |
|---|
| 293 | return [wlon,wlat] |
|---|
| 294 | |
|---|
| 295 | ## Author: AS |
|---|
| 296 | def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False): |
|---|
| 297 | import matplotlib.pyplot as plt |
|---|
| 298 | from os import system |
|---|
| 299 | addstr = "" |
|---|
| 300 | if res is not None: |
|---|
| 301 | res = int(res) |
|---|
| 302 | addstr = "_"+str(res) |
|---|
| 303 | name = filename+addstr+"."+ext |
|---|
| 304 | if folder != '': name = folder+'/'+name |
|---|
| 305 | plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value) |
|---|
| 306 | if disp: display(name) |
|---|
| 307 | if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name) |
|---|
| 308 | if erase: system("mv "+name+" to_be_erased") |
|---|
| 309 | return |
|---|
| 310 | |
|---|
| 311 | ## Author: AS |
|---|
| 312 | def dumpbdy (field,n,stag=None): |
|---|
| 313 | nx = len(field[0,:])-1 |
|---|
| 314 | ny = len(field[:,0])-1 |
|---|
| 315 | if stag == 'U': nx = nx-1 |
|---|
| 316 | if stag == 'V': ny = ny-1 |
|---|
| 317 | return field[n:ny-n,n:nx-n] |
|---|
| 318 | |
|---|
| 319 | ## Author: AS |
|---|
| 320 | def getcoorddef ( nc ): |
|---|
| 321 | import numpy as np |
|---|
| 322 | ## getcoord2d for predefined types |
|---|
| 323 | typefile = whatkindfile(nc) |
|---|
| 324 | if typefile in ['mesoapi','meso']: |
|---|
| 325 | [lon2d,lat2d] = getcoord2d(nc) |
|---|
| 326 | lon2d = dumpbdy(lon2d,6) |
|---|
| 327 | lat2d = dumpbdy(lat2d,6) |
|---|
| 328 | elif typefile in ['gcm']: |
|---|
| 329 | [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) |
|---|
| 330 | elif typefile in ['geo']: |
|---|
| 331 | [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') |
|---|
| 332 | return lon2d,lat2d |
|---|
| 333 | |
|---|
| 334 | ## Author: AS |
|---|
| 335 | def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): |
|---|
| 336 | import numpy as np |
|---|
| 337 | if is1d: |
|---|
| 338 | lat = nc.variables[nlat][:] |
|---|
| 339 | lon = nc.variables[nlon][:] |
|---|
| 340 | [lon2d,lat2d] = np.meshgrid(lon,lat) |
|---|
| 341 | else: |
|---|
| 342 | lat = nc.variables[nlat][0,:,:] |
|---|
| 343 | lon = nc.variables[nlon][0,:,:] |
|---|
| 344 | [lon2d,lat2d] = [lon,lat] |
|---|
| 345 | return lon2d,lat2d |
|---|
| 346 | |
|---|
| 347 | ## Author: AS |
|---|
| 348 | def smooth (field, coeff): |
|---|
| 349 | ## actually blur_image could work with different coeff on x and y |
|---|
| 350 | if coeff > 1: result = blur_image(field,int(coeff)) |
|---|
| 351 | else: result = field |
|---|
| 352 | return result |
|---|
| 353 | |
|---|
| 354 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
|---|
| 355 | def gauss_kern(size, sizey=None): |
|---|
| 356 | import numpy as np |
|---|
| 357 | # Returns a normalized 2D gauss kernel array for convolutions |
|---|
| 358 | size = int(size) |
|---|
| 359 | if not sizey: |
|---|
| 360 | sizey = size |
|---|
| 361 | else: |
|---|
| 362 | sizey = int(sizey) |
|---|
| 363 | x, y = np.mgrid[-size:size+1, -sizey:sizey+1] |
|---|
| 364 | g = np.exp(-(x**2/float(size)+y**2/float(sizey))) |
|---|
| 365 | return g / g.sum() |
|---|
| 366 | |
|---|
| 367 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
|---|
| 368 | def blur_image(im, n, ny=None) : |
|---|
| 369 | from scipy.signal import convolve |
|---|
| 370 | # blurs the image by convolving with a gaussian kernel of typical size n. |
|---|
| 371 | # The optional keyword argument ny allows for a different size in the y direction. |
|---|
| 372 | g = gauss_kern(n, sizey=ny) |
|---|
| 373 | improc = convolve(im, g, mode='same') |
|---|
| 374 | return improc |
|---|
| 375 | |
|---|
| 376 | ## Author: AS |
|---|
| 377 | def getwinddef (nc): |
|---|
| 378 | ## getwinds for predefined types |
|---|
| 379 | typefile = whatkindfile(nc) |
|---|
| 380 | ### |
|---|
| 381 | if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm'] |
|---|
| 382 | elif typefile is 'gcm': [uchar,vchar] = ['u','v'] |
|---|
| 383 | elif typefile is 'meso': [uchar,vchar] = ['U','V'] |
|---|
| 384 | else: [uchar,vchar] = ['not found','not found'] |
|---|
| 385 | ### |
|---|
| 386 | if typefile in ['meso']: metwind = False ## geometrical (wrt grid) |
|---|
| 387 | else: metwind = True ## meteorological (zon/mer) |
|---|
| 388 | if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)" |
|---|
| 389 | ### |
|---|
| 390 | return uchar,vchar,metwind |
|---|
| 391 | |
|---|
| 392 | ## Author: AS |
|---|
| 393 | def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True): |
|---|
| 394 | ## scale regle la reference du vecteur |
|---|
| 395 | ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs. |
|---|
| 396 | import matplotlib.pyplot as plt |
|---|
| 397 | import numpy as np |
|---|
| 398 | posx = np.min(x) - np.std(x) / 10. |
|---|
| 399 | posy = np.min(y) - np.std(y) / 10. |
|---|
| 400 | u = smooth(u,csmooth) |
|---|
| 401 | v = smooth(v,csmooth) |
|---|
| 402 | widthvec = 0.003 #0.005 #0.003 |
|---|
| 403 | q = plt.quiver( x[::stride,::stride],\ |
|---|
| 404 | y[::stride,::stride],\ |
|---|
| 405 | u[::stride,::stride],\ |
|---|
| 406 | v[::stride,::stride],\ |
|---|
| 407 | angles='xy',color=color,pivot='middle',\ |
|---|
| 408 | scale=factor,width=widthvec ) |
|---|
| 409 | if color in ['white','yellow']: kcolor='black' |
|---|
| 410 | else: kcolor=color |
|---|
| 411 | if key: p = plt.quiverkey(q,posx,posy,scale,\ |
|---|
| 412 | str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03) |
|---|
| 413 | return |
|---|
| 414 | |
|---|
| 415 | ## Author: AS |
|---|
| 416 | def display (name): |
|---|
| 417 | from os import system |
|---|
| 418 | system("display "+name+" > /dev/null 2> /dev/null &") |
|---|
| 419 | return name |
|---|
| 420 | |
|---|
| 421 | ## Author: AS |
|---|
| 422 | def findstep (wlon): |
|---|
| 423 | steplon = int((wlon[1]-wlon[0])/4.) #3 |
|---|
| 424 | step = 120. |
|---|
| 425 | while step > steplon and step > 15. : step = step / 2. |
|---|
| 426 | if step <= 15.: |
|---|
| 427 | while step > steplon and step > 5. : step = step - 5. |
|---|
| 428 | if step <= 5.: |
|---|
| 429 | while step > steplon and step > 1. : step = step - 1. |
|---|
| 430 | if step <= 1.: |
|---|
| 431 | step = 1. |
|---|
| 432 | return step |
|---|
| 433 | |
|---|
| 434 | ## Author: AS |
|---|
| 435 | def define_proj (char,wlon,wlat,back=None,blat=False): |
|---|
| 436 | from mpl_toolkits.basemap import Basemap |
|---|
| 437 | import numpy as np |
|---|
| 438 | import matplotlib as mpl |
|---|
| 439 | from mymath import max |
|---|
| 440 | meanlon = 0.5*(wlon[0]+wlon[1]) |
|---|
| 441 | meanlat = 0.5*(wlat[0]+wlat[1]) |
|---|
| 442 | if not blat: |
|---|
| 443 | if wlat[0] >= 80.: blat = 40. |
|---|
| 444 | elif wlat[1] <= -80.: blat = -40. |
|---|
| 445 | elif wlat[1] >= 0.: blat = wlat[0] |
|---|
| 446 | elif wlat[0] <= 0.: blat = wlat[1] |
|---|
| 447 | print "blat ", blat |
|---|
| 448 | h = 50. ## en km |
|---|
| 449 | radius = 3397200. |
|---|
| 450 | if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\ |
|---|
| 451 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 452 | elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon) |
|---|
| 453 | elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat) |
|---|
| 454 | elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\ |
|---|
| 455 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 456 | elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.) |
|---|
| 457 | elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.) |
|---|
| 458 | elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon) |
|---|
| 459 | elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\ |
|---|
| 460 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 461 | elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.) |
|---|
| 462 | elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\ |
|---|
| 463 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
|---|
| 464 | fontsizemer = int(mpl.rcParams['font.size']*3./4.) |
|---|
| 465 | if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon) |
|---|
| 466 | else: step = 10. |
|---|
| 467 | steplon = step*2. |
|---|
| 468 | #if back in ["geolocal"]: |
|---|
| 469 | # step = np.min([5.,step]) |
|---|
| 470 | # steplon = step |
|---|
| 471 | print step |
|---|
| 472 | m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer) |
|---|
| 473 | m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer) |
|---|
| 474 | if back: m.warpimage(marsmap(back),scale=0.75) |
|---|
| 475 | #if not back: |
|---|
| 476 | # if not var: back = "mola" ## if no var: draw mola |
|---|
| 477 | # elif typefile in ['mesoapi','meso','geo'] \ |
|---|
| 478 | # and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw |
|---|
| 479 | # else: pass ## else: draw None |
|---|
| 480 | return m |
|---|
| 481 | |
|---|
| 482 | ## Author: AS |
|---|
| 483 | #### test temporaire |
|---|
| 484 | def putpoints (map,plot): |
|---|
| 485 | #### from http://www.scipy.org/Cookbook/Matplotlib/Maps |
|---|
| 486 | # lat/lon coordinates of five cities. |
|---|
| 487 | lats = [18.4] |
|---|
| 488 | lons = [-134.0] |
|---|
| 489 | points=['Olympus Mons'] |
|---|
| 490 | # compute the native map projection coordinates for cities. |
|---|
| 491 | x,y = map(lons,lats) |
|---|
| 492 | # plot filled circles at the locations of the cities. |
|---|
| 493 | map.plot(x,y,'bo') |
|---|
| 494 | # plot the names of those five cities. |
|---|
| 495 | wherept = 0 #1000 #50000 |
|---|
| 496 | for name,xpt,ypt in zip(points,x,y): |
|---|
| 497 | plot.text(xpt+wherept,ypt+wherept,name) |
|---|
| 498 | ## le nom ne s'affiche pas... |
|---|
| 499 | return |
|---|
| 500 | |
|---|
| 501 | ## Author: AS |
|---|
| 502 | def calculate_bounds(field,vmin=None,vmax=None): |
|---|
| 503 | import numpy as np |
|---|
| 504 | from mymath import max,min,mean |
|---|
| 505 | ind = np.where(field < 9e+35) |
|---|
| 506 | fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple |
|---|
| 507 | ### |
|---|
| 508 | dev = np.std(fieldcalc)*3.0 |
|---|
| 509 | ### |
|---|
| 510 | if vmin is None: |
|---|
| 511 | zevmin = mean(fieldcalc) - dev |
|---|
| 512 | else: zevmin = vmin |
|---|
| 513 | ### |
|---|
| 514 | if vmax is None: zevmax = mean(fieldcalc) + dev |
|---|
| 515 | else: zevmax = vmax |
|---|
| 516 | if vmin == vmax: |
|---|
| 517 | zevmin = mean(fieldcalc) - dev ### for continuity |
|---|
| 518 | zevmax = mean(fieldcalc) + dev ### for continuity |
|---|
| 519 | ### |
|---|
| 520 | if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0. |
|---|
| 521 | print "field ", min(fieldcalc), max(fieldcalc) |
|---|
| 522 | print "bounds ", zevmin, zevmax |
|---|
| 523 | return zevmin, zevmax |
|---|
| 524 | |
|---|
| 525 | ## Author: AS |
|---|
| 526 | def bounds(what_I_plot,zevmin,zevmax): |
|---|
| 527 | from mymath import max,min,mean |
|---|
| 528 | ### might be convenient to add the missing value in arguments |
|---|
| 529 | #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7) |
|---|
| 530 | if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7) |
|---|
| 531 | else: what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7) |
|---|
| 532 | print "new min ", min(what_I_plot) |
|---|
| 533 | what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 |
|---|
| 534 | what_I_plot[ what_I_plot > zevmax ] = zevmax |
|---|
| 535 | print "new max ", max(what_I_plot) |
|---|
| 536 | |
|---|
| 537 | return what_I_plot |
|---|
| 538 | |
|---|
| 539 | ## Author: AS |
|---|
| 540 | def nolow(what_I_plot): |
|---|
| 541 | from mymath import max,min |
|---|
| 542 | lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot))) |
|---|
| 543 | print "on vire en dessous de ", lim |
|---|
| 544 | what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 |
|---|
| 545 | return what_I_plot |
|---|
| 546 | |
|---|
| 547 | ## Author: AS |
|---|
| 548 | def zoomset (wlon,wlat,zoom): |
|---|
| 549 | dlon = abs(wlon[1]-wlon[0])/2. |
|---|
| 550 | dlat = abs(wlat[1]-wlat[0])/2. |
|---|
| 551 | [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\ |
|---|
| 552 | [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ] |
|---|
| 553 | print "zoom %",zoom,wlon,wlat |
|---|
| 554 | return wlon,wlat |
|---|
| 555 | |
|---|
| 556 | ## Author: AS |
|---|
| 557 | def fmtvar (whichvar="def"): |
|---|
| 558 | fmtvar = { \ |
|---|
| 559 | "tk": "%.0f",\ |
|---|
| 560 | "tpot": "%.0f",\ |
|---|
| 561 | "TSURF": "%.0f",\ |
|---|
| 562 | "def": "%.1e",\ |
|---|
| 563 | "PTOT": "%.0f",\ |
|---|
| 564 | "HGT": "%.1e",\ |
|---|
| 565 | "USTM": "%.2f",\ |
|---|
| 566 | "HFX": "%.0f",\ |
|---|
| 567 | "ICETOT": "%.1e",\ |
|---|
| 568 | "TAU_ICE": "%.2f",\ |
|---|
| 569 | "VMR_ICE": "%.1e",\ |
|---|
| 570 | "MTOT": "%.1f",\ |
|---|
| 571 | "anomaly": "%.1f",\ |
|---|
| 572 | "W": "%.1f",\ |
|---|
| 573 | "WMAX_TH": "%.1f",\ |
|---|
| 574 | "QSURFICE": "%.0f",\ |
|---|
| 575 | "Um": "%.0f",\ |
|---|
| 576 | "ALBBARE": "%.2f",\ |
|---|
| 577 | "TAU": "%.1f",\ |
|---|
| 578 | #### T.N. |
|---|
| 579 | "TEMP": "%.0f",\ |
|---|
| 580 | "VMR_H2OICE": "%.0f",\ |
|---|
| 581 | "VMR_H2OVAP": "%.0f",\ |
|---|
| 582 | "TAUTES": "%.2f",\ |
|---|
| 583 | "TAUTESAP": "%.2f",\ |
|---|
| 584 | |
|---|
| 585 | } |
|---|
| 586 | if whichvar not in fmtvar: |
|---|
| 587 | whichvar = "def" |
|---|
| 588 | return fmtvar[whichvar] |
|---|
| 589 | |
|---|
| 590 | ## Author: AS |
|---|
| 591 | #################################################################################################################### |
|---|
| 592 | ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png |
|---|
| 593 | def defcolorb (whichone="def"): |
|---|
| 594 | whichcolorb = { \ |
|---|
| 595 | "def": "spectral",\ |
|---|
| 596 | "HGT": "spectral",\ |
|---|
| 597 | "tk": "gist_heat",\ |
|---|
| 598 | "TSURF": "RdBu_r",\ |
|---|
| 599 | "QH2O": "PuBu",\ |
|---|
| 600 | "USTM": "YlOrRd",\ |
|---|
| 601 | "HFX": "RdYlBu",\ |
|---|
| 602 | "ICETOT": "YlGnBu_r",\ |
|---|
| 603 | #"MTOT": "PuBu",\ |
|---|
| 604 | "CCNQ": "YlOrBr",\ |
|---|
| 605 | "CCNN": "YlOrBr",\ |
|---|
| 606 | "TEMP": "Jet",\ |
|---|
| 607 | "TAU_ICE": "Blues",\ |
|---|
| 608 | "VMR_ICE": "Blues",\ |
|---|
| 609 | "W": "jet",\ |
|---|
| 610 | "WMAX_TH": "spectral",\ |
|---|
| 611 | "anomaly": "RdBu_r",\ |
|---|
| 612 | "QSURFICE": "hot_r",\ |
|---|
| 613 | "ALBBARE": "spectral",\ |
|---|
| 614 | "TAU": "YlOrBr_r",\ |
|---|
| 615 | #### T.N. |
|---|
| 616 | "MTOT": "Jet",\ |
|---|
| 617 | "H2O_ICE_S": "RdBu",\ |
|---|
| 618 | "VMR_H2OICE": "PuBu",\ |
|---|
| 619 | "VMR_H2OVAP": "PuBu",\ |
|---|
| 620 | } |
|---|
| 621 | #W --> spectral ou jet |
|---|
| 622 | #spectral BrBG RdBu_r |
|---|
| 623 | print "predefined colorbars" |
|---|
| 624 | if whichone not in whichcolorb: |
|---|
| 625 | whichone = "def" |
|---|
| 626 | return whichcolorb[whichone] |
|---|
| 627 | |
|---|
| 628 | ## Author: AS |
|---|
| 629 | def definecolorvec (whichone="def"): |
|---|
| 630 | whichcolor = { \ |
|---|
| 631 | "def": "black",\ |
|---|
| 632 | "vis": "yellow",\ |
|---|
| 633 | "vishires": "yellow",\ |
|---|
| 634 | "molabw": "yellow",\ |
|---|
| 635 | "mola": "black",\ |
|---|
| 636 | "gist_heat": "white",\ |
|---|
| 637 | "hot": "tk",\ |
|---|
| 638 | "gist_rainbow": "black",\ |
|---|
| 639 | "spectral": "black",\ |
|---|
| 640 | "gray": "red",\ |
|---|
| 641 | "PuBu": "black",\ |
|---|
| 642 | } |
|---|
| 643 | if whichone not in whichcolor: |
|---|
| 644 | whichone = "def" |
|---|
| 645 | return whichcolor[whichone] |
|---|
| 646 | |
|---|
| 647 | ## Author: AS |
|---|
| 648 | def marsmap (whichone="vishires"): |
|---|
| 649 | from os import uname |
|---|
| 650 | mymachine = uname()[1] |
|---|
| 651 | ### not sure about speed-up with this method... looks the same |
|---|
| 652 | if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/" |
|---|
| 653 | else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/" |
|---|
| 654 | whichlink = { \ |
|---|
| 655 | #"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ |
|---|
| 656 | #"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\ |
|---|
| 657 | #"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\ |
|---|
| 658 | #"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ |
|---|
| 659 | #"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ |
|---|
| 660 | "vis": domain+"mar0kuu2.jpg",\ |
|---|
| 661 | "vishires": domain+"MarsMap_2500x1250.jpg",\ |
|---|
| 662 | "geolocal": domain+"geolocal.jpg",\ |
|---|
| 663 | "mola": domain+"mars-mola-2k.jpg",\ |
|---|
| 664 | "molabw": domain+"MarsElevation_2500x1250.jpg",\ |
|---|
| 665 | "clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\ |
|---|
| 666 | "jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\ |
|---|
| 667 | "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\ |
|---|
| 668 | "bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\ |
|---|
| 669 | "contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\ |
|---|
| 670 | "nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\ |
|---|
| 671 | "blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\ |
|---|
| 672 | "blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\ |
|---|
| 673 | "justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\ |
|---|
| 674 | } |
|---|
| 675 | ### see http://www.mmedia.is/~bjj/planetary_maps.html |
|---|
| 676 | if whichone not in whichlink: |
|---|
| 677 | print "marsmap: choice not defined... you'll get the default one... " |
|---|
| 678 | whichone = "vishires" |
|---|
| 679 | return whichlink[whichone] |
|---|
| 680 | |
|---|
| 681 | #def earthmap (whichone): |
|---|
| 682 | # if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg" |
|---|
| 683 | # elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg" |
|---|
| 684 | # elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg" |
|---|
| 685 | # return whichlink |
|---|
| 686 | |
|---|
| 687 | ## Author: AS |
|---|
| 688 | def latinterv (area="Whole"): |
|---|
| 689 | list = { \ |
|---|
| 690 | "Europe": [[ 20., 80.],[- 50., 50.]],\ |
|---|
| 691 | "Central_America": [[-10., 40.],[ 230., 300.]],\ |
|---|
| 692 | "Africa": [[-20., 50.],[- 50., 50.]],\ |
|---|
| 693 | "Whole": [[-90., 90.],[-180., 180.]],\ |
|---|
| 694 | "Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\ |
|---|
| 695 | "Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\ |
|---|
| 696 | "Tharsis": [[-30., 60.],[-170.,- 10.]],\ |
|---|
| 697 | "Whole_No_High": [[-60., 60.],[-180., 180.]],\ |
|---|
| 698 | "Chryse": [[-60., 60.],[- 60., 60.]],\ |
|---|
| 699 | "North_Pole": [[ 50., 90.],[-180., 180.]],\ |
|---|
| 700 | "Close_North_Pole": [[ 75., 90.],[-180., 180.]],\ |
|---|
| 701 | "Far_South_Pole": [[-90.,-40.],[-180., 180.]],\ |
|---|
| 702 | "South_Pole": [[-90.,-50.],[-180., 180.]],\ |
|---|
| 703 | "Close_South_Pole": [[-90.,-75.],[-180., 180.]],\ |
|---|
| 704 | } |
|---|
| 705 | if area not in list: area = "Whole" |
|---|
| 706 | [olat,olon] = list[area] |
|---|
| 707 | return olon,olat |
|---|
| 708 | |
|---|
| 709 | ## Author: TN |
|---|
| 710 | def separatenames (name): |
|---|
| 711 | from numpy import concatenate |
|---|
| 712 | # look for comas in the input name to separate different names (files, variables,etc ..) |
|---|
| 713 | if name is None: |
|---|
| 714 | names = None |
|---|
| 715 | else: |
|---|
| 716 | names = [] |
|---|
| 717 | stop = 0 |
|---|
| 718 | currentname = name |
|---|
| 719 | while stop == 0: |
|---|
| 720 | indexvir = currentname.find(',') |
|---|
| 721 | if indexvir == -1: |
|---|
| 722 | stop = 1 |
|---|
| 723 | name1 = currentname |
|---|
| 724 | else: |
|---|
| 725 | name1 = currentname[0:indexvir] |
|---|
| 726 | names = concatenate((names,[name1])) |
|---|
| 727 | currentname = currentname[indexvir+1:len(currentname)] |
|---|
| 728 | return names |
|---|
| 729 | |
|---|
| 730 | ## Author: TN [old] |
|---|
| 731 | def getopmatrix (kind,n): |
|---|
| 732 | import numpy as np |
|---|
| 733 | # compute matrix of operations between files |
|---|
| 734 | # the matrix is 'number of files'-square |
|---|
| 735 | # 1: difference (row minus column), 2: add |
|---|
| 736 | # not 0 in diag : just plot |
|---|
| 737 | if n == 1: |
|---|
| 738 | opm = np.eye(1) |
|---|
| 739 | elif kind == 'basic': |
|---|
| 740 | opm = np.eye(n) |
|---|
| 741 | elif kind == 'difference1': # show differences with 1st file |
|---|
| 742 | opm = np.zeros((n,n)) |
|---|
| 743 | opm[0,:] = 1 |
|---|
| 744 | opm[0,0] = 0 |
|---|
| 745 | elif kind == 'difference2': # show differences with 1st file AND show 1st file |
|---|
| 746 | opm = np.zeros((n,n)) |
|---|
| 747 | opm[0,:] = 1 |
|---|
| 748 | else: |
|---|
| 749 | opm = np.eye(n) |
|---|
| 750 | return opm |
|---|
| 751 | |
|---|
| 752 | ## Author: TN [old] |
|---|
| 753 | def checkcoherence (nfiles,nlat,nlon,ntime): |
|---|
| 754 | if (nfiles > 1): |
|---|
| 755 | if (nlat > 1): |
|---|
| 756 | errormess("what you asked is not possible !") |
|---|
| 757 | return 1 |
|---|
| 758 | |
|---|
| 759 | ## Author: TN |
|---|
| 760 | def readslices(saxis): |
|---|
| 761 | from numpy import empty |
|---|
| 762 | if saxis == None: |
|---|
| 763 | zesaxis = None |
|---|
| 764 | else: |
|---|
| 765 | zesaxis = empty((len(saxis),2)) |
|---|
| 766 | for i in range(len(saxis)): |
|---|
| 767 | a = separatenames(saxis[i]) |
|---|
| 768 | if len(a) == 1: |
|---|
| 769 | zesaxis[i,:] = float(a[0]) |
|---|
| 770 | else: |
|---|
| 771 | zesaxis[i,0] = float(a[0]) |
|---|
| 772 | zesaxis[i,1] = float(a[1]) |
|---|
| 773 | |
|---|
| 774 | return zesaxis |
|---|
| 775 | |
|---|
| 776 | ## Author: TN |
|---|
| 777 | def getsindex(saxis,index,axis): |
|---|
| 778 | # input : all the desired slices and the good index |
|---|
| 779 | # output : all indexes to be taken into account for reducing field |
|---|
| 780 | import numpy as np |
|---|
| 781 | if ( np.array(axis).ndim == 2): |
|---|
| 782 | axis = axis[:,0] |
|---|
| 783 | if saxis is None: |
|---|
| 784 | zeindex = None |
|---|
| 785 | else: |
|---|
| 786 | aaa = int(np.argmin(abs(saxis[index,0] - axis))) |
|---|
| 787 | bbb = int(np.argmin(abs(saxis[index,1] - axis))) |
|---|
| 788 | [imin,imax] = np.sort(np.array([aaa,bbb])) |
|---|
| 789 | zeindex = np.array(range(imax-imin+1))+imin |
|---|
| 790 | # because -180 and 180 are the same point in longitude, |
|---|
| 791 | # we get rid of one for averaging purposes. |
|---|
| 792 | if axis[imin] == -180 and axis[imax] == 180: |
|---|
| 793 | zeindex = zeindex[0:len(zeindex)-1] |
|---|
| 794 | print "whole longitude averaging asked, so last point is not taken into account." |
|---|
| 795 | return zeindex |
|---|
| 796 | |
|---|
| 797 | ## Author: TN |
|---|
| 798 | def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode): |
|---|
| 799 | # Purpose of define_axis is to find x and y axis scales in a smart way |
|---|
| 800 | # x axis priority: 1/time 2/lon 3/lat 4/vertical |
|---|
| 801 | # To be improved !!!... |
|---|
| 802 | from numpy import array,swapaxes |
|---|
| 803 | x = None |
|---|
| 804 | y = None |
|---|
| 805 | count = 0 |
|---|
| 806 | what_I_plot = array(what_I_plot) |
|---|
| 807 | shape = what_I_plot.shape |
|---|
| 808 | if indextime is None: |
|---|
| 809 | print "axis is time" |
|---|
| 810 | x = time |
|---|
| 811 | count = count+1 |
|---|
| 812 | if indexlon is None: |
|---|
| 813 | print "axis is lon" |
|---|
| 814 | if count == 0: x = lon |
|---|
| 815 | else: y = lon |
|---|
| 816 | count = count+1 |
|---|
| 817 | if indexlat is None: |
|---|
| 818 | print "axis is lat" |
|---|
| 819 | if count == 0: x = lat |
|---|
| 820 | else: y = lat |
|---|
| 821 | count = count+1 |
|---|
| 822 | if indexvert is None and dim0 is 4: |
|---|
| 823 | print "axis is vert" |
|---|
| 824 | if vertmode == 0: # vertical axis is as is (GCM grid) |
|---|
| 825 | if count == 0: x=range(len(vert)) |
|---|
| 826 | else: y=range(len(vert)) |
|---|
| 827 | count = count+1 |
|---|
| 828 | else: # vertical axis is in kms |
|---|
| 829 | if count == 0: x = vert |
|---|
| 830 | else: y = vert |
|---|
| 831 | count = count+1 |
|---|
| 832 | x = array(x) |
|---|
| 833 | y = array(y) |
|---|
| 834 | print "what_I_plot.shape", what_I_plot.shape |
|---|
| 835 | print "x.shape, y.shape", x.shape, y.shape |
|---|
| 836 | if len(shape) == 1: |
|---|
| 837 | print shape[0] |
|---|
| 838 | if shape[0] != len(x): |
|---|
| 839 | print "WARNING HERE !!!" |
|---|
| 840 | x = y |
|---|
| 841 | elif len(shape) == 2: |
|---|
| 842 | print shape[1], len(y), shape[0], len(x) |
|---|
| 843 | if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]: |
|---|
| 844 | what_I_plot = swapaxes(what_I_plot,0,1) |
|---|
| 845 | print "swapaxes", what_I_plot.shape, shape |
|---|
| 846 | #x0 = x |
|---|
| 847 | #x = y |
|---|
| 848 | #y = x0 |
|---|
| 849 | #print "define_axis", x, y |
|---|
| 850 | return what_I_plot,x,y |
|---|
| 851 | |
|---|
| 852 | # Author: TN + AS |
|---|
| 853 | def determineplot(slon, slat, svert, stime): |
|---|
| 854 | nlon = 1 # number of longitudinal slices -- 1 is None |
|---|
| 855 | nlat = 1 |
|---|
| 856 | nvert = 1 |
|---|
| 857 | ntime = 1 |
|---|
| 858 | nslices = 1 |
|---|
| 859 | if slon is not None: |
|---|
| 860 | nslices = nslices*len(slon) |
|---|
| 861 | nlon = len(slon) |
|---|
| 862 | if slat is not None: |
|---|
| 863 | nslices = nslices*len(slat) |
|---|
| 864 | nlat = len(slat) |
|---|
| 865 | if svert is not None: |
|---|
| 866 | nslices = nslices*len(svert) |
|---|
| 867 | nvert = len(svert) |
|---|
| 868 | if stime is not None: |
|---|
| 869 | nslices = nslices*len(stime) |
|---|
| 870 | ntime = len(stime) |
|---|
| 871 | #else: |
|---|
| 872 | # nslices = 2 |
|---|
| 873 | |
|---|
| 874 | mapmode = 0 |
|---|
| 875 | if slon is None and slat is None: |
|---|
| 876 | mapmode = 1 # in this case we plot a map, with the given projection |
|---|
| 877 | #elif proj is not None: |
|---|
| 878 | # print "WARNING: you specified a", proj,\ |
|---|
| 879 | # "projection but asked for slices", slon,"in longitude and",slat,"in latitude" |
|---|
| 880 | print "mapmode: ", mapmode |
|---|
| 881 | |
|---|
| 882 | return nlon, nlat, nvert, ntime, mapmode, nslices |
|---|
