[345] | 1 | ## Author: AS |
---|
[252] | 2 | def errormess(text,printvar=None): |
---|
[233] | 3 | print text |
---|
[399] | 4 | if printvar is not None: print printvar |
---|
[233] | 5 | exit() |
---|
| 6 | return |
---|
| 7 | |
---|
[345] | 8 | ## Author: AS |
---|
[349] | 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 |
---|
[252] | 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 | |
---|
[345] | 30 | ## Author: AS |
---|
[252] | 31 | def localtime(utc,lon): |
---|
| 32 | ltst = utc + lon / 15. |
---|
| 33 | ltst = int (ltst * 10) / 10. |
---|
| 34 | ltst = ltst % 24 |
---|
| 35 | return ltst |
---|
| 36 | |
---|
[428] | 37 | ## Author: AS, AC |
---|
[233] | 38 | def whatkindfile (nc): |
---|
[429] | 39 | if 'controle' in nc.variables: typefile = 'gcm' |
---|
| 40 | elif 'phisinit' in nc.variables: typefile = 'gcm' |
---|
[432] | 41 | elif hasattr(nc,'START_DATE'): |
---|
| 42 | if '9999' in getattr(nc,'START_DATE') : typefile = 'mesoideal' |
---|
[434] | 43 | elif 'vert' in nc.variables: typefile = 'mesoapi' |
---|
| 44 | elif 'U' in nc.variables: typefile = 'meso' |
---|
[429] | 45 | elif 'HGT_M' in nc.variables: typefile = 'geo' |
---|
[345] | 46 | #else: errormess("whatkindfile: typefile not supported.") |
---|
[429] | 47 | else: typefile = 'gcm' # for lslin-ed files from gcm |
---|
[233] | 48 | return typefile |
---|
| 49 | |
---|
[345] | 50 | ## Author: AS |
---|
[233] | 51 | def getfield (nc,var): |
---|
| 52 | ## this allows to get much faster (than simply referring to nc.variables[var]) |
---|
[395] | 53 | import numpy as np |
---|
[233] | 54 | dimension = len(nc.variables[var].dimensions) |
---|
[392] | 55 | #print " Opening variable with", dimension, "dimensions ..." |
---|
[233] | 56 | if dimension == 2: field = nc.variables[var][:,:] |
---|
| 57 | elif dimension == 3: field = nc.variables[var][:,:,:] |
---|
| 58 | elif dimension == 4: field = nc.variables[var][:,:,:,:] |
---|
[395] | 59 | # if there are NaNs in the ncdf, they should be loaded as a masked array which will be |
---|
| 60 | # recasted as a regular array later in reducefield |
---|
| 61 | if (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')): |
---|
| 62 | print "Warning: netcdf as nan values but is not loaded as a Masked Array." |
---|
| 63 | print "recasting array type" |
---|
| 64 | out=np.ma.masked_invalid(field) |
---|
| 65 | out.set_fill_value([np.NaN]) |
---|
| 66 | else: |
---|
| 67 | out=field |
---|
| 68 | return out |
---|
[233] | 69 | |
---|
[382] | 70 | ## Author: AS + TN + AC |
---|
[405] | 71 | def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False): |
---|
[252] | 72 | ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x" |
---|
[233] | 73 | ### it would be actually better to name d4 d3 d2 d1 as t z y x |
---|
[405] | 74 | ### ... note, anomaly is only computed over d1 and d2 for the moment |
---|
[233] | 75 | import numpy as np |
---|
[349] | 76 | from mymath import max,mean |
---|
[422] | 77 | csmooth = 12 ## a fair amount of grid points (too high results in high computation time) |
---|
[233] | 78 | dimension = np.array(input).ndim |
---|
| 79 | shape = np.array(input).shape |
---|
[349] | 80 | #print 'd1,d2,d3,d4: ',d1,d2,d3,d4 |
---|
[392] | 81 | print 'IN REDUCEFIELD dim,shape: ',dimension,shape |
---|
[405] | 82 | if anomaly: print 'ANOMALY ANOMALY' |
---|
[233] | 83 | output = input |
---|
| 84 | error = False |
---|
[350] | 85 | #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays |
---|
[345] | 86 | if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4] |
---|
| 87 | if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3] |
---|
| 88 | if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2] |
---|
| 89 | if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1] |
---|
| 90 | ### now the main part |
---|
[233] | 91 | if dimension == 2: |
---|
| 92 | if d2 >= shape[0]: error = True |
---|
| 93 | elif d1 >= shape[1]: error = True |
---|
[350] | 94 | elif d1 is not None and d2 is not None: output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0) |
---|
| 95 | elif d1 is not None: output = mean(input[:,d1],axis=1) |
---|
| 96 | elif d2 is not None: output = mean(input[d2,:],axis=0) |
---|
[233] | 97 | elif dimension == 3: |
---|
[345] | 98 | if max(d4) >= shape[0]: error = True |
---|
| 99 | elif max(d2) >= shape[1]: error = True |
---|
| 100 | elif max(d1) >= shape[2]: error = True |
---|
[350] | 101 | elif d4 is not None and d2 is not None and d1 is not None: |
---|
| 102 | output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0) |
---|
[349] | 103 | elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0) |
---|
| 104 | elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1) |
---|
| 105 | elif d2 is not None and d1 is not None: output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1) |
---|
| 106 | elif d1 is not None: output = mean(input[:,:,d1],axis=2) |
---|
| 107 | elif d2 is not None: output = mean(input[:,d2,:],axis=1) |
---|
| 108 | elif d4 is not None: output = mean(input[d4,:,:],axis=0) |
---|
[233] | 109 | elif dimension == 4: |
---|
[345] | 110 | if max(d4) >= shape[0]: error = True |
---|
| 111 | elif max(d3) >= shape[1]: error = True |
---|
| 112 | elif max(d2) >= shape[2]: error = True |
---|
| 113 | elif max(d1) >= shape[3]: error = True |
---|
[382] | 114 | elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: |
---|
[392] | 115 | output = mean(input[d4,:,:,:],axis=0) |
---|
| 116 | output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) |
---|
[427] | 117 | if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) |
---|
[392] | 118 | output = mean(output[d2,:],axis=0) |
---|
| 119 | output = mean(output[d1],axis=0) |
---|
[350] | 120 | elif d4 is not None and d3 is not None and d2 is not None: |
---|
[392] | 121 | output = mean(input[d4,:,:,:],axis=0) |
---|
| 122 | output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) |
---|
[405] | 123 | if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) |
---|
[392] | 124 | output = mean(output[d2,:],axis=0) |
---|
[350] | 125 | elif d4 is not None and d3 is not None and d1 is not None: |
---|
[392] | 126 | output = mean(input[d4,:,:,:],axis=0) |
---|
| 127 | output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) |
---|
[405] | 128 | if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) |
---|
[392] | 129 | output = mean(output[:,d1],axis=1) |
---|
[350] | 130 | elif d4 is not None and d2 is not None and d1 is not None: |
---|
[392] | 131 | output = mean(input[d4,:,:,:],axis=0) |
---|
[405] | 132 | if anomaly: |
---|
| 133 | for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.) |
---|
[392] | 134 | output = mean(output[:,d2,:],axis=1) |
---|
| 135 | output = mean(output[:,d1],axis=1) |
---|
[405] | 136 | #noperturb = smooth1d(output,window_len=7) |
---|
| 137 | #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4] |
---|
| 138 | #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show() |
---|
[350] | 139 | elif d3 is not None and d2 is not None and d1 is not None: |
---|
[392] | 140 | output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) |
---|
[405] | 141 | if anomaly: |
---|
| 142 | for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.) |
---|
[392] | 143 | output = mean(output[:,d2,:],axis=1) |
---|
| 144 | output = mean(output[:,d1],axis=1) |
---|
| 145 | elif d4 is not None and d3 is not None: |
---|
| 146 | output = mean(input[d4,:,:,:],axis=0) |
---|
| 147 | output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) |
---|
[405] | 148 | if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) |
---|
[392] | 149 | elif d4 is not None and d2 is not None: |
---|
| 150 | output = mean(input[d4,:,:,:],axis=0) |
---|
| 151 | output = mean(output[:,d2,:],axis=1) |
---|
| 152 | elif d4 is not None and d1 is not None: |
---|
| 153 | output = mean(input[d4,:,:,:],axis=0) |
---|
| 154 | output = mean(output[:,:,d1],axis=2) |
---|
| 155 | elif d3 is not None and d2 is not None: |
---|
| 156 | output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) |
---|
| 157 | output = mean(output[:,d2,:],axis=1) |
---|
| 158 | elif d3 is not None and d1 is not None: |
---|
| 159 | output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) |
---|
[448] | 160 | output = mean(output[:,:,d1],axis=2) |
---|
[392] | 161 | elif d2 is not None and d1 is not None: |
---|
| 162 | output = mean(input[:,:,d2,:],axis=2) |
---|
| 163 | output = mean(output[:,:,d1],axis=2) |
---|
| 164 | elif d1 is not None: output = mean(input[:,:,:,d1],axis=3) |
---|
| 165 | elif d2 is not None: output = mean(input[:,:,d2,:],axis=2) |
---|
[437] | 166 | elif d3 is not None: output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) |
---|
[392] | 167 | elif d4 is not None: output = mean(input[d4,:,:,:],axis=0) |
---|
[233] | 168 | dimension = np.array(output).ndim |
---|
| 169 | shape = np.array(output).shape |
---|
[392] | 170 | print 'OUT REDUCEFIELD dim,shape: ',dimension,shape |
---|
[233] | 171 | return output, error |
---|
| 172 | |
---|
[392] | 173 | ## Author: AC + AS |
---|
| 174 | def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None): |
---|
[382] | 175 | from mymath import max,mean |
---|
| 176 | from scipy import integrate |
---|
[392] | 177 | if yint and vert is not None and indice is not None: |
---|
[391] | 178 | if type(input).__name__=='MaskedArray': |
---|
| 179 | input.set_fill_value([np.NaN]) |
---|
[392] | 180 | output = integrate.trapz(input.filled(),x=vert[indice],axis=ax) |
---|
[391] | 181 | else: |
---|
[396] | 182 | output = integrate.trapz(input,x=vert[indice],axis=ax) |
---|
[382] | 183 | else: |
---|
| 184 | output = mean(input,axis=ax) |
---|
| 185 | return output |
---|
| 186 | |
---|
[345] | 187 | ## Author: AS + TN |
---|
[233] | 188 | def definesubplot ( numplot, fig ): |
---|
| 189 | from matplotlib.pyplot import rcParams |
---|
| 190 | rcParams['font.size'] = 12. ## default (important for multiple calls) |
---|
[345] | 191 | if numplot <= 0: |
---|
| 192 | subv = 99999 |
---|
| 193 | subh = 99999 |
---|
| 194 | elif numplot == 1: |
---|
| 195 | subv = 99999 |
---|
| 196 | subh = 99999 |
---|
[233] | 197 | elif numplot == 2: |
---|
[453] | 198 | subv = 2 |
---|
| 199 | subh = 1 |
---|
[233] | 200 | fig.subplots_adjust(wspace = 0.35) |
---|
| 201 | rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. ) |
---|
| 202 | elif numplot == 3: |
---|
[453] | 203 | subv = 3 |
---|
| 204 | subh = 1 |
---|
[233] | 205 | fig.subplots_adjust(wspace = 0.5) |
---|
| 206 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
[345] | 207 | elif numplot == 4: |
---|
| 208 | subv = 2 |
---|
| 209 | subh = 2 |
---|
| 210 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 211 | rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) |
---|
| 212 | elif numplot <= 6: |
---|
| 213 | subv = 2 |
---|
| 214 | subh = 3 |
---|
[233] | 215 | fig.subplots_adjust(wspace = 0.4, hspace = 0.0) |
---|
| 216 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
[345] | 217 | elif numplot <= 8: |
---|
| 218 | subv = 2 |
---|
| 219 | subh = 4 |
---|
[233] | 220 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 221 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
[345] | 222 | elif numplot <= 9: |
---|
| 223 | subv = 3 |
---|
| 224 | subh = 3 |
---|
[233] | 225 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 226 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
[345] | 227 | elif numplot <= 12: |
---|
| 228 | subv = 3 |
---|
| 229 | subh = 4 |
---|
[453] | 230 | fig.subplots_adjust(wspace = 0, hspace = 0.1) |
---|
[345] | 231 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
| 232 | elif numplot <= 16: |
---|
| 233 | subv = 4 |
---|
| 234 | subh = 4 |
---|
| 235 | fig.subplots_adjust(wspace = 0.3, hspace = 0.3) |
---|
| 236 | rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) |
---|
[233] | 237 | else: |
---|
[345] | 238 | print "number of plot supported: 1 to 16" |
---|
[233] | 239 | exit() |
---|
[345] | 240 | return subv,subh |
---|
[233] | 241 | |
---|
[345] | 242 | ## Author: AS |
---|
[233] | 243 | def getstralt(nc,nvert): |
---|
| 244 | typefile = whatkindfile(nc) |
---|
| 245 | if typefile is 'meso': |
---|
| 246 | stralt = "_lvl" + str(nvert) |
---|
| 247 | elif typefile is 'mesoapi': |
---|
| 248 | zelevel = int(nc.variables['vert'][nvert]) |
---|
| 249 | if abs(zelevel) < 10000.: strheight=str(zelevel)+"m" |
---|
| 250 | else: strheight=str(int(zelevel/1000.))+"km" |
---|
| 251 | if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR" |
---|
| 252 | elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS" |
---|
| 253 | elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight |
---|
| 254 | elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" |
---|
| 255 | else: stralt = "" |
---|
| 256 | else: |
---|
| 257 | stralt = "" |
---|
| 258 | return stralt |
---|
| 259 | |
---|
[345] | 260 | ## Author: AS |
---|
[195] | 261 | def getlschar ( namefile ): |
---|
| 262 | from netCDF4 import Dataset |
---|
| 263 | from timestuff import sol2ls |
---|
[233] | 264 | from numpy import array |
---|
[400] | 265 | from string import rstrip |
---|
[405] | 266 | namefile = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p") |
---|
[400] | 267 | #### we assume that wrfout is next to wrfout_z and wrfout_zabg |
---|
[195] | 268 | nc = Dataset(namefile) |
---|
[237] | 269 | zetime = None |
---|
[400] | 270 | if 'Times' in nc.variables: |
---|
[233] | 271 | zetime = nc.variables['Times'][0] |
---|
| 272 | shape = array(nc.variables['Times']).shape |
---|
| 273 | if shape[0] < 2: zetime = None |
---|
| 274 | if zetime is not None \ |
---|
[225] | 275 | and 'vert' not in nc.variables: |
---|
[317] | 276 | #### strangely enough this does not work for api or ncrcat results! |
---|
[195] | 277 | zetimestart = getattr(nc, 'START_DATE') |
---|
| 278 | zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) |
---|
| 279 | if zeday < 0: lschar="" ## might have crossed a month... fix soon |
---|
[241] | 280 | else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. ) |
---|
[197] | 281 | ### |
---|
| 282 | zetime2 = nc.variables['Times'][1] |
---|
| 283 | one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37. |
---|
| 284 | next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. |
---|
| 285 | zehour = one |
---|
| 286 | zehourin = abs ( next - one ) |
---|
[195] | 287 | else: |
---|
| 288 | lschar="" |
---|
[197] | 289 | zehour = 0 |
---|
| 290 | zehourin = 1 |
---|
| 291 | return lschar, zehour, zehourin |
---|
[195] | 292 | |
---|
[345] | 293 | ## Author: AS |
---|
[202] | 294 | def getprefix (nc): |
---|
| 295 | prefix = 'LMD_MMM_' |
---|
| 296 | prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_' |
---|
| 297 | prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_' |
---|
| 298 | return prefix |
---|
| 299 | |
---|
[345] | 300 | ## Author: AS |
---|
[184] | 301 | def getproj (nc): |
---|
[233] | 302 | typefile = whatkindfile(nc) |
---|
| 303 | if typefile in ['mesoapi','meso','geo']: |
---|
| 304 | ### (il faudrait passer CEN_LON dans la projection ?) |
---|
| 305 | map_proj = getattr(nc, 'MAP_PROJ') |
---|
| 306 | cen_lat = getattr(nc, 'CEN_LAT') |
---|
| 307 | if map_proj == 2: |
---|
| 308 | if cen_lat > 10.: |
---|
| 309 | proj="npstere" |
---|
[392] | 310 | #print "NP stereographic polar domain" |
---|
[233] | 311 | else: |
---|
| 312 | proj="spstere" |
---|
[392] | 313 | #print "SP stereographic polar domain" |
---|
[233] | 314 | elif map_proj == 1: |
---|
[392] | 315 | #print "lambert projection domain" |
---|
[233] | 316 | proj="lcc" |
---|
| 317 | elif map_proj == 3: |
---|
[392] | 318 | #print "mercator projection" |
---|
[233] | 319 | proj="merc" |
---|
| 320 | else: |
---|
| 321 | proj="merc" |
---|
[252] | 322 | elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?) |
---|
[233] | 323 | else: proj="ortho" |
---|
[184] | 324 | return proj |
---|
| 325 | |
---|
[345] | 326 | ## Author: AS |
---|
[180] | 327 | def ptitle (name): |
---|
| 328 | from matplotlib.pyplot import title |
---|
| 329 | title(name) |
---|
| 330 | print name |
---|
| 331 | |
---|
[345] | 332 | ## Author: AS |
---|
[252] | 333 | def polarinterv (lon2d,lat2d): |
---|
| 334 | import numpy as np |
---|
| 335 | wlon = [np.min(lon2d),np.max(lon2d)] |
---|
| 336 | ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole |
---|
| 337 | wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])] |
---|
| 338 | return [wlon,wlat] |
---|
| 339 | |
---|
[345] | 340 | ## Author: AS |
---|
[180] | 341 | def simplinterv (lon2d,lat2d): |
---|
| 342 | import numpy as np |
---|
| 343 | return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]] |
---|
| 344 | |
---|
[345] | 345 | ## Author: AS |
---|
[184] | 346 | def wrfinterv (lon2d,lat2d): |
---|
| 347 | nx = len(lon2d[0,:])-1 |
---|
| 348 | ny = len(lon2d[:,0])-1 |
---|
[225] | 349 | lon1 = lon2d[0,0] |
---|
| 350 | lon2 = lon2d[nx,ny] |
---|
| 351 | lat1 = lat2d[0,0] |
---|
| 352 | lat2 = lat2d[nx,ny] |
---|
[233] | 353 | if abs(0.5*(lat1+lat2)) > 60.: wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 |
---|
| 354 | else: wider = 0. |
---|
| 355 | if lon1 < lon2: wlon = [lon1, lon2 + wider] |
---|
[225] | 356 | else: wlon = [lon2, lon1 + wider] |
---|
| 357 | if lat1 < lat2: wlat = [lat1, lat2] |
---|
| 358 | else: wlat = [lat2, lat1] |
---|
| 359 | return [wlon,wlat] |
---|
[184] | 360 | |
---|
[345] | 361 | ## Author: AS |
---|
[240] | 362 | def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False): |
---|
[180] | 363 | import matplotlib.pyplot as plt |
---|
[240] | 364 | from os import system |
---|
| 365 | addstr = "" |
---|
| 366 | if res is not None: |
---|
| 367 | res = int(res) |
---|
| 368 | addstr = "_"+str(res) |
---|
| 369 | name = filename+addstr+"."+ext |
---|
[186] | 370 | if folder != '': name = folder+'/'+name |
---|
[180] | 371 | plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value) |
---|
[240] | 372 | if disp: display(name) |
---|
| 373 | if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name) |
---|
| 374 | if erase: system("mv "+name+" to_be_erased") |
---|
[180] | 375 | return |
---|
| 376 | |
---|
[430] | 377 | ## Author: AS + AC |
---|
[451] | 378 | def dumpbdy (field,n,stag=None,condition=False,onlyx=False,onlyy=False): |
---|
[447] | 379 | nx = len(field[0,:])-1 |
---|
| 380 | ny = len(field[:,0])-1 |
---|
[444] | 381 | if condition: |
---|
| 382 | if stag == 'U': nx = nx-1 |
---|
| 383 | if stag == 'V': ny = ny-1 |
---|
| 384 | if stag == 'W': nx = nx+1 #special les case when we dump stag on W |
---|
[451] | 385 | if onlyx: result = field[:,n:nx-n] |
---|
| 386 | elif onlyy: result = field[n:ny-n,:] |
---|
| 387 | else: result = field[n:ny-n,n:nx-n] |
---|
| 388 | return result |
---|
[180] | 389 | |
---|
[444] | 390 | ## Author: AS + AC |
---|
[233] | 391 | def getcoorddef ( nc ): |
---|
[317] | 392 | import numpy as np |
---|
[233] | 393 | ## getcoord2d for predefined types |
---|
| 394 | typefile = whatkindfile(nc) |
---|
| 395 | if typefile in ['mesoapi','meso']: |
---|
| 396 | [lon2d,lat2d] = getcoord2d(nc) |
---|
[240] | 397 | lon2d = dumpbdy(lon2d,6) |
---|
| 398 | lat2d = dumpbdy(lat2d,6) |
---|
[317] | 399 | elif typefile in ['gcm']: |
---|
[233] | 400 | [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) |
---|
| 401 | elif typefile in ['geo']: |
---|
| 402 | [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') |
---|
[429] | 403 | elif typefile in ['mesoideal']: |
---|
[428] | 404 | nx=getattr(nc,'WEST-EAST_GRID_DIMENSION') |
---|
| 405 | ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION') |
---|
| 406 | [lon2d,lat2d] = np.meshgrid(np.arange(nx),np.arange(ny)) |
---|
[233] | 407 | return lon2d,lat2d |
---|
| 408 | |
---|
[345] | 409 | ## Author: AS |
---|
[184] | 410 | def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): |
---|
| 411 | import numpy as np |
---|
| 412 | if is1d: |
---|
| 413 | lat = nc.variables[nlat][:] |
---|
| 414 | lon = nc.variables[nlon][:] |
---|
| 415 | [lon2d,lat2d] = np.meshgrid(lon,lat) |
---|
| 416 | else: |
---|
| 417 | lat = nc.variables[nlat][0,:,:] |
---|
| 418 | lon = nc.variables[nlon][0,:,:] |
---|
| 419 | [lon2d,lat2d] = [lon,lat] |
---|
| 420 | return lon2d,lat2d |
---|
| 421 | |
---|
[405] | 422 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
---|
| 423 | def smooth1d(x,window_len=11,window='hanning'): |
---|
| 424 | import numpy |
---|
| 425 | """smooth the data using a window with requested size. |
---|
| 426 | This method is based on the convolution of a scaled window with the signal. |
---|
| 427 | The signal is prepared by introducing reflected copies of the signal |
---|
| 428 | (with the window size) in both ends so that transient parts are minimized |
---|
| 429 | in the begining and end part of the output signal. |
---|
| 430 | input: |
---|
| 431 | x: the input signal |
---|
| 432 | window_len: the dimension of the smoothing window; should be an odd integer |
---|
| 433 | window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman' |
---|
| 434 | flat window will produce a moving average smoothing. |
---|
| 435 | output: |
---|
| 436 | the smoothed signal |
---|
| 437 | example: |
---|
| 438 | t=linspace(-2,2,0.1) |
---|
| 439 | x=sin(t)+randn(len(t))*0.1 |
---|
| 440 | y=smooth(x) |
---|
| 441 | see also: |
---|
| 442 | numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve |
---|
| 443 | scipy.signal.lfilter |
---|
| 444 | TODO: the window parameter could be the window itself if an array instead of a string |
---|
| 445 | """ |
---|
| 446 | x = numpy.array(x) |
---|
| 447 | if x.ndim != 1: |
---|
| 448 | raise ValueError, "smooth only accepts 1 dimension arrays." |
---|
| 449 | if x.size < window_len: |
---|
| 450 | raise ValueError, "Input vector needs to be bigger than window size." |
---|
| 451 | if window_len<3: |
---|
| 452 | return x |
---|
| 453 | if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: |
---|
| 454 | raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'" |
---|
| 455 | s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]] |
---|
| 456 | #print(len(s)) |
---|
| 457 | if window == 'flat': #moving average |
---|
| 458 | w=numpy.ones(window_len,'d') |
---|
| 459 | else: |
---|
| 460 | w=eval('numpy.'+window+'(window_len)') |
---|
| 461 | y=numpy.convolve(w/w.sum(),s,mode='valid') |
---|
| 462 | return y |
---|
| 463 | |
---|
[345] | 464 | ## Author: AS |
---|
[180] | 465 | def smooth (field, coeff): |
---|
| 466 | ## actually blur_image could work with different coeff on x and y |
---|
| 467 | if coeff > 1: result = blur_image(field,int(coeff)) |
---|
| 468 | else: result = field |
---|
| 469 | return result |
---|
| 470 | |
---|
[345] | 471 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
---|
[180] | 472 | def gauss_kern(size, sizey=None): |
---|
| 473 | import numpy as np |
---|
| 474 | # Returns a normalized 2D gauss kernel array for convolutions |
---|
| 475 | size = int(size) |
---|
| 476 | if not sizey: |
---|
| 477 | sizey = size |
---|
| 478 | else: |
---|
| 479 | sizey = int(sizey) |
---|
| 480 | x, y = np.mgrid[-size:size+1, -sizey:sizey+1] |
---|
| 481 | g = np.exp(-(x**2/float(size)+y**2/float(sizey))) |
---|
| 482 | return g / g.sum() |
---|
| 483 | |
---|
[345] | 484 | ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth |
---|
[180] | 485 | def blur_image(im, n, ny=None) : |
---|
| 486 | from scipy.signal import convolve |
---|
| 487 | # blurs the image by convolving with a gaussian kernel of typical size n. |
---|
| 488 | # The optional keyword argument ny allows for a different size in the y direction. |
---|
| 489 | g = gauss_kern(n, sizey=ny) |
---|
| 490 | improc = convolve(im, g, mode='same') |
---|
| 491 | return improc |
---|
| 492 | |
---|
[345] | 493 | ## Author: AS |
---|
[233] | 494 | def getwinddef (nc): |
---|
| 495 | ## getwinds for predefined types |
---|
| 496 | typefile = whatkindfile(nc) |
---|
| 497 | ### |
---|
[429] | 498 | if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm'] |
---|
| 499 | elif typefile is 'gcm': [uchar,vchar] = ['u','v'] |
---|
| 500 | elif typefile in ['meso','mesoideal']: [uchar,vchar] = ['U','V'] |
---|
| 501 | else: [uchar,vchar] = ['not found','not found'] |
---|
[233] | 502 | ### |
---|
| 503 | if typefile in ['meso']: metwind = False ## geometrical (wrt grid) |
---|
| 504 | else: metwind = True ## meteorological (zon/mer) |
---|
| 505 | if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)" |
---|
| 506 | ### |
---|
| 507 | return uchar,vchar,metwind |
---|
[202] | 508 | |
---|
[345] | 509 | ## Author: AS |
---|
[184] | 510 | def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True): |
---|
| 511 | ## scale regle la reference du vecteur |
---|
| 512 | ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs. |
---|
| 513 | import matplotlib.pyplot as plt |
---|
| 514 | import numpy as np |
---|
[187] | 515 | posx = np.min(x) - np.std(x) / 10. |
---|
| 516 | posy = np.min(y) - np.std(y) / 10. |
---|
[184] | 517 | u = smooth(u,csmooth) |
---|
| 518 | v = smooth(v,csmooth) |
---|
[188] | 519 | widthvec = 0.003 #0.005 #0.003 |
---|
[184] | 520 | q = plt.quiver( x[::stride,::stride],\ |
---|
| 521 | y[::stride,::stride],\ |
---|
| 522 | u[::stride,::stride],\ |
---|
| 523 | v[::stride,::stride],\ |
---|
[228] | 524 | angles='xy',color=color,pivot='middle',\ |
---|
[184] | 525 | scale=factor,width=widthvec ) |
---|
[202] | 526 | if color in ['white','yellow']: kcolor='black' |
---|
| 527 | else: kcolor=color |
---|
[184] | 528 | if key: p = plt.quiverkey(q,posx,posy,scale,\ |
---|
[194] | 529 | str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03) |
---|
[184] | 530 | return |
---|
[180] | 531 | |
---|
[345] | 532 | ## Author: AS |
---|
[180] | 533 | def display (name): |
---|
[184] | 534 | from os import system |
---|
| 535 | system("display "+name+" > /dev/null 2> /dev/null &") |
---|
| 536 | return name |
---|
[180] | 537 | |
---|
[345] | 538 | ## Author: AS |
---|
[180] | 539 | def findstep (wlon): |
---|
[184] | 540 | steplon = int((wlon[1]-wlon[0])/4.) #3 |
---|
| 541 | step = 120. |
---|
| 542 | while step > steplon and step > 15. : step = step / 2. |
---|
| 543 | if step <= 15.: |
---|
| 544 | while step > steplon and step > 5. : step = step - 5. |
---|
| 545 | if step <= 5.: |
---|
| 546 | while step > steplon and step > 1. : step = step - 1. |
---|
| 547 | if step <= 1.: |
---|
| 548 | step = 1. |
---|
[180] | 549 | return step |
---|
| 550 | |
---|
[345] | 551 | ## Author: AS |
---|
[451] | 552 | def define_proj (char,wlon,wlat,back=None,blat=None,blon=None): |
---|
[180] | 553 | from mpl_toolkits.basemap import Basemap |
---|
| 554 | import numpy as np |
---|
| 555 | import matplotlib as mpl |
---|
[240] | 556 | from mymath import max |
---|
[180] | 557 | meanlon = 0.5*(wlon[0]+wlon[1]) |
---|
| 558 | meanlat = 0.5*(wlat[0]+wlat[1]) |
---|
[385] | 559 | if blat is None: |
---|
[398] | 560 | ortholat=meanlat |
---|
[453] | 561 | if wlat[0] >= 80.: blat = -40. |
---|
[345] | 562 | elif wlat[1] <= -80.: blat = -40. |
---|
| 563 | elif wlat[1] >= 0.: blat = wlat[0] |
---|
| 564 | elif wlat[0] <= 0.: blat = wlat[1] |
---|
[398] | 565 | else: ortholat=blat |
---|
[451] | 566 | if blon is None: ortholon=meanlon |
---|
| 567 | else: ortholon=blon |
---|
[207] | 568 | h = 50. ## en km |
---|
[202] | 569 | radius = 3397200. |
---|
[184] | 570 | if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\ |
---|
[180] | 571 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
[184] | 572 | elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon) |
---|
[451] | 573 | elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=ortholon,lat_0=ortholat) |
---|
[184] | 574 | elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\ |
---|
| 575 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
| 576 | elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.) |
---|
[395] | 577 | elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=180.) |
---|
[207] | 578 | elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon) |
---|
| 579 | elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\ |
---|
| 580 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
[184] | 581 | elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.) |
---|
| 582 | elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\ |
---|
| 583 | llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) |
---|
| 584 | fontsizemer = int(mpl.rcParams['font.size']*3./4.) |
---|
[207] | 585 | if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon) |
---|
| 586 | else: step = 10. |
---|
[238] | 587 | steplon = step*2. |
---|
[453] | 588 | zecolor ='grey' |
---|
| 589 | zelinewidth = 1 |
---|
| 590 | zelatmax = 80 |
---|
| 591 | # to show gcm grid: |
---|
| 592 | #zecolor = 'r' |
---|
| 593 | #zelinewidth = 1 |
---|
| 594 | #step = 5.625 |
---|
| 595 | #steplon = 5.625 |
---|
| 596 | #zelatmax = 89.9 |
---|
| 597 | m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax) |
---|
| 598 | m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax) |
---|
[233] | 599 | if back: m.warpimage(marsmap(back),scale=0.75) |
---|
| 600 | #if not back: |
---|
| 601 | # if not var: back = "mola" ## if no var: draw mola |
---|
| 602 | # elif typefile in ['mesoapi','meso','geo'] \ |
---|
| 603 | # and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw |
---|
| 604 | # else: pass ## else: draw None |
---|
[180] | 605 | return m |
---|
| 606 | |
---|
[345] | 607 | ## Author: AS |
---|
[232] | 608 | #### test temporaire |
---|
| 609 | def putpoints (map,plot): |
---|
| 610 | #### from http://www.scipy.org/Cookbook/Matplotlib/Maps |
---|
| 611 | # lat/lon coordinates of five cities. |
---|
| 612 | lats = [18.4] |
---|
| 613 | lons = [-134.0] |
---|
| 614 | points=['Olympus Mons'] |
---|
| 615 | # compute the native map projection coordinates for cities. |
---|
| 616 | x,y = map(lons,lats) |
---|
| 617 | # plot filled circles at the locations of the cities. |
---|
| 618 | map.plot(x,y,'bo') |
---|
| 619 | # plot the names of those five cities. |
---|
| 620 | wherept = 0 #1000 #50000 |
---|
| 621 | for name,xpt,ypt in zip(points,x,y): |
---|
| 622 | plot.text(xpt+wherept,ypt+wherept,name) |
---|
| 623 | ## le nom ne s'affiche pas... |
---|
| 624 | return |
---|
| 625 | |
---|
[345] | 626 | ## Author: AS |
---|
[233] | 627 | def calculate_bounds(field,vmin=None,vmax=None): |
---|
| 628 | import numpy as np |
---|
| 629 | from mymath import max,min,mean |
---|
| 630 | ind = np.where(field < 9e+35) |
---|
| 631 | fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple |
---|
| 632 | ### |
---|
| 633 | dev = np.std(fieldcalc)*3.0 |
---|
| 634 | ### |
---|
| 635 | if vmin is None: |
---|
| 636 | zevmin = mean(fieldcalc) - dev |
---|
| 637 | else: zevmin = vmin |
---|
| 638 | ### |
---|
| 639 | if vmax is None: zevmax = mean(fieldcalc) + dev |
---|
| 640 | else: zevmax = vmax |
---|
| 641 | if vmin == vmax: |
---|
| 642 | zevmin = mean(fieldcalc) - dev ### for continuity |
---|
| 643 | zevmax = mean(fieldcalc) + dev ### for continuity |
---|
| 644 | ### |
---|
| 645 | if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0. |
---|
[392] | 646 | print "BOUNDS field ", min(fieldcalc), max(fieldcalc) |
---|
| 647 | print "BOUNDS adopted ", zevmin, zevmax |
---|
[233] | 648 | return zevmin, zevmax |
---|
[232] | 649 | |
---|
[345] | 650 | ## Author: AS |
---|
[233] | 651 | def bounds(what_I_plot,zevmin,zevmax): |
---|
[247] | 652 | from mymath import max,min,mean |
---|
[233] | 653 | ### might be convenient to add the missing value in arguments |
---|
[310] | 654 | #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7) |
---|
| 655 | if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7) |
---|
| 656 | else: what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7) |
---|
[451] | 657 | #print "NEW MIN ", min(what_I_plot) |
---|
[233] | 658 | what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 |
---|
[310] | 659 | what_I_plot[ what_I_plot > zevmax ] = zevmax |
---|
[451] | 660 | #print "NEW MAX ", max(what_I_plot) |
---|
[233] | 661 | return what_I_plot |
---|
| 662 | |
---|
[345] | 663 | ## Author: AS |
---|
[241] | 664 | def nolow(what_I_plot): |
---|
| 665 | from mymath import max,min |
---|
| 666 | lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot))) |
---|
[392] | 667 | print "NO PLOT BELOW VALUE ", lim |
---|
[241] | 668 | what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 |
---|
| 669 | return what_I_plot |
---|
| 670 | |
---|
[418] | 671 | |
---|
| 672 | ## Author : AC |
---|
| 673 | def hole_bounds(what_I_plot,zevmin,zevmax): |
---|
| 674 | import numpy as np |
---|
| 675 | zi=0 |
---|
| 676 | for i in what_I_plot: |
---|
| 677 | zj=0 |
---|
| 678 | for j in i: |
---|
| 679 | if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN |
---|
| 680 | zj=zj+1 |
---|
| 681 | zi=zi+1 |
---|
| 682 | |
---|
| 683 | return what_I_plot |
---|
| 684 | |
---|
[345] | 685 | ## Author: AS |
---|
[233] | 686 | def zoomset (wlon,wlat,zoom): |
---|
| 687 | dlon = abs(wlon[1]-wlon[0])/2. |
---|
| 688 | dlat = abs(wlat[1]-wlat[0])/2. |
---|
| 689 | [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\ |
---|
| 690 | [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ] |
---|
[392] | 691 | print "ZOOM %",zoom,wlon,wlat |
---|
[233] | 692 | return wlon,wlat |
---|
| 693 | |
---|
[345] | 694 | ## Author: AS |
---|
[201] | 695 | def fmtvar (whichvar="def"): |
---|
[204] | 696 | fmtvar = { \ |
---|
[405] | 697 | "TK": "%.0f",\ |
---|
[425] | 698 | # Variables from TES ncdf format |
---|
[363] | 699 | "T_NADIR_DAY": "%.0f",\ |
---|
[376] | 700 | "T_NADIR_NIT": "%.0f",\ |
---|
[425] | 701 | # Variables from tes.py ncdf format |
---|
[398] | 702 | "TEMP_DAY": "%.0f",\ |
---|
| 703 | "TEMP_NIGHT": "%.0f",\ |
---|
[425] | 704 | # Variables from MCS and mcs.py ncdf format |
---|
[427] | 705 | "DTEMP": "%.0f",\ |
---|
| 706 | "NTEMP": "%.0f",\ |
---|
| 707 | "DNUMBINTEMP": "%.0f",\ |
---|
| 708 | "NNUMBINTEMP": "%.0f",\ |
---|
[425] | 709 | # other stuff |
---|
[405] | 710 | "TPOT": "%.0f",\ |
---|
[295] | 711 | "TSURF": "%.0f",\ |
---|
[204] | 712 | "def": "%.1e",\ |
---|
| 713 | "PTOT": "%.0f",\ |
---|
| 714 | "HGT": "%.1e",\ |
---|
| 715 | "USTM": "%.2f",\ |
---|
[225] | 716 | "HFX": "%.0f",\ |
---|
[232] | 717 | "ICETOT": "%.1e",\ |
---|
[237] | 718 | "TAU_ICE": "%.2f",\ |
---|
[451] | 719 | "TAUICE": "%.2f",\ |
---|
[252] | 720 | "VMR_ICE": "%.1e",\ |
---|
[345] | 721 | "MTOT": "%.1f",\ |
---|
[405] | 722 | "ANOMALY": "%.1f",\ |
---|
[241] | 723 | "W": "%.1f",\ |
---|
[287] | 724 | "WMAX_TH": "%.1f",\ |
---|
| 725 | "QSURFICE": "%.0f",\ |
---|
[405] | 726 | "UM": "%.0f",\ |
---|
[295] | 727 | "ALBBARE": "%.2f",\ |
---|
[317] | 728 | "TAU": "%.1f",\ |
---|
[382] | 729 | "CO2": "%.2f",\ |
---|
[345] | 730 | #### T.N. |
---|
| 731 | "TEMP": "%.0f",\ |
---|
| 732 | "VMR_H2OICE": "%.0f",\ |
---|
| 733 | "VMR_H2OVAP": "%.0f",\ |
---|
| 734 | "TAUTES": "%.2f",\ |
---|
| 735 | "TAUTESAP": "%.2f",\ |
---|
| 736 | |
---|
[204] | 737 | } |
---|
| 738 | if whichvar not in fmtvar: |
---|
| 739 | whichvar = "def" |
---|
| 740 | return fmtvar[whichvar] |
---|
[201] | 741 | |
---|
[345] | 742 | ## Author: AS |
---|
[233] | 743 | #################################################################################################################### |
---|
| 744 | ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png |
---|
[202] | 745 | def defcolorb (whichone="def"): |
---|
[204] | 746 | whichcolorb = { \ |
---|
| 747 | "def": "spectral",\ |
---|
| 748 | "HGT": "spectral",\ |
---|
[426] | 749 | "HGT_M": "spectral",\ |
---|
[405] | 750 | "TK": "gist_heat",\ |
---|
[425] | 751 | "TPOT": "Paired",\ |
---|
[295] | 752 | "TSURF": "RdBu_r",\ |
---|
[204] | 753 | "QH2O": "PuBu",\ |
---|
| 754 | "USTM": "YlOrRd",\ |
---|
[363] | 755 | #"T_nadir_nit": "RdBu_r",\ |
---|
| 756 | #"T_nadir_day": "RdBu_r",\ |
---|
[225] | 757 | "HFX": "RdYlBu",\ |
---|
[310] | 758 | "ICETOT": "YlGnBu_r",\ |
---|
[345] | 759 | #"MTOT": "PuBu",\ |
---|
| 760 | "CCNQ": "YlOrBr",\ |
---|
| 761 | "CCNN": "YlOrBr",\ |
---|
| 762 | "TEMP": "Jet",\ |
---|
[238] | 763 | "TAU_ICE": "Blues",\ |
---|
[451] | 764 | "TAUICE": "Blues",\ |
---|
[252] | 765 | "VMR_ICE": "Blues",\ |
---|
[241] | 766 | "W": "jet",\ |
---|
[287] | 767 | "WMAX_TH": "spectral",\ |
---|
[405] | 768 | "ANOMALY": "RdBu_r",\ |
---|
[287] | 769 | "QSURFICE": "hot_r",\ |
---|
[295] | 770 | "ALBBARE": "spectral",\ |
---|
[317] | 771 | "TAU": "YlOrBr_r",\ |
---|
[382] | 772 | "CO2": "YlOrBr_r",\ |
---|
[345] | 773 | #### T.N. |
---|
| 774 | "MTOT": "Jet",\ |
---|
| 775 | "H2O_ICE_S": "RdBu",\ |
---|
| 776 | "VMR_H2OICE": "PuBu",\ |
---|
| 777 | "VMR_H2OVAP": "PuBu",\ |
---|
[453] | 778 | "WATERCAPTAG": "Blues",\ |
---|
[204] | 779 | } |
---|
[241] | 780 | #W --> spectral ou jet |
---|
[240] | 781 | #spectral BrBG RdBu_r |
---|
[392] | 782 | #print "predefined colorbars" |
---|
[204] | 783 | if whichone not in whichcolorb: |
---|
| 784 | whichone = "def" |
---|
| 785 | return whichcolorb[whichone] |
---|
[202] | 786 | |
---|
[345] | 787 | ## Author: AS |
---|
[202] | 788 | def definecolorvec (whichone="def"): |
---|
| 789 | whichcolor = { \ |
---|
| 790 | "def": "black",\ |
---|
| 791 | "vis": "yellow",\ |
---|
| 792 | "vishires": "yellow",\ |
---|
| 793 | "molabw": "yellow",\ |
---|
| 794 | "mola": "black",\ |
---|
| 795 | "gist_heat": "white",\ |
---|
| 796 | "hot": "tk",\ |
---|
| 797 | "gist_rainbow": "black",\ |
---|
| 798 | "spectral": "black",\ |
---|
| 799 | "gray": "red",\ |
---|
| 800 | "PuBu": "black",\ |
---|
| 801 | } |
---|
| 802 | if whichone not in whichcolor: |
---|
| 803 | whichone = "def" |
---|
| 804 | return whichcolor[whichone] |
---|
| 805 | |
---|
[345] | 806 | ## Author: AS |
---|
[180] | 807 | def marsmap (whichone="vishires"): |
---|
[233] | 808 | from os import uname |
---|
| 809 | mymachine = uname()[1] |
---|
| 810 | ### not sure about speed-up with this method... looks the same |
---|
| 811 | if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/" |
---|
| 812 | else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/" |
---|
[180] | 813 | whichlink = { \ |
---|
[233] | 814 | #"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ |
---|
| 815 | #"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\ |
---|
| 816 | #"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\ |
---|
| 817 | #"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ |
---|
| 818 | #"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ |
---|
[453] | 819 | "thermalday": domain+"thermalday.jpg",\ |
---|
| 820 | "thermalnight": domain+"thermalnight.jpg",\ |
---|
| 821 | "tesalbedo": domain+"tesalbedo.jpg",\ |
---|
[233] | 822 | "vis": domain+"mar0kuu2.jpg",\ |
---|
| 823 | "vishires": domain+"MarsMap_2500x1250.jpg",\ |
---|
| 824 | "geolocal": domain+"geolocal.jpg",\ |
---|
| 825 | "mola": domain+"mars-mola-2k.jpg",\ |
---|
| 826 | "molabw": domain+"MarsElevation_2500x1250.jpg",\ |
---|
[238] | 827 | "clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\ |
---|
| 828 | "jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\ |
---|
| 829 | "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\ |
---|
[273] | 830 | "bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\ |
---|
| 831 | "contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\ |
---|
| 832 | "nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\ |
---|
| 833 | "blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\ |
---|
[296] | 834 | "blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\ |
---|
| 835 | "justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\ |
---|
[180] | 836 | } |
---|
[238] | 837 | ### see http://www.mmedia.is/~bjj/planetary_maps.html |
---|
[180] | 838 | if whichone not in whichlink: |
---|
| 839 | print "marsmap: choice not defined... you'll get the default one... " |
---|
| 840 | whichone = "vishires" |
---|
| 841 | return whichlink[whichone] |
---|
| 842 | |
---|
[273] | 843 | #def earthmap (whichone): |
---|
| 844 | # if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg" |
---|
| 845 | # elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg" |
---|
| 846 | # elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg" |
---|
| 847 | # return whichlink |
---|
[180] | 848 | |
---|
[345] | 849 | ## Author: AS |
---|
[241] | 850 | def latinterv (area="Whole"): |
---|
| 851 | list = { \ |
---|
| 852 | "Europe": [[ 20., 80.],[- 50., 50.]],\ |
---|
| 853 | "Central_America": [[-10., 40.],[ 230., 300.]],\ |
---|
| 854 | "Africa": [[-20., 50.],[- 50., 50.]],\ |
---|
[273] | 855 | "Whole": [[-90., 90.],[-180., 180.]],\ |
---|
| 856 | "Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\ |
---|
| 857 | "Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\ |
---|
[241] | 858 | "Tharsis": [[-30., 60.],[-170.,- 10.]],\ |
---|
| 859 | "Whole_No_High": [[-60., 60.],[-180., 180.]],\ |
---|
| 860 | "Chryse": [[-60., 60.],[- 60., 60.]],\ |
---|
| 861 | "North_Pole": [[ 50., 90.],[-180., 180.]],\ |
---|
| 862 | "Close_North_Pole": [[ 75., 90.],[-180., 180.]],\ |
---|
| 863 | "Far_South_Pole": [[-90.,-40.],[-180., 180.]],\ |
---|
| 864 | "South_Pole": [[-90.,-50.],[-180., 180.]],\ |
---|
| 865 | "Close_South_Pole": [[-90.,-75.],[-180., 180.]],\ |
---|
| 866 | } |
---|
| 867 | if area not in list: area = "Whole" |
---|
| 868 | [olat,olon] = list[area] |
---|
| 869 | return olon,olat |
---|
| 870 | |
---|
[345] | 871 | ## Author: TN |
---|
| 872 | def separatenames (name): |
---|
| 873 | from numpy import concatenate |
---|
| 874 | # look for comas in the input name to separate different names (files, variables,etc ..) |
---|
| 875 | if name is None: |
---|
| 876 | names = None |
---|
| 877 | else: |
---|
| 878 | names = [] |
---|
| 879 | stop = 0 |
---|
| 880 | currentname = name |
---|
| 881 | while stop == 0: |
---|
| 882 | indexvir = currentname.find(',') |
---|
| 883 | if indexvir == -1: |
---|
| 884 | stop = 1 |
---|
| 885 | name1 = currentname |
---|
| 886 | else: |
---|
| 887 | name1 = currentname[0:indexvir] |
---|
| 888 | names = concatenate((names,[name1])) |
---|
| 889 | currentname = currentname[indexvir+1:len(currentname)] |
---|
| 890 | return names |
---|
| 891 | |
---|
| 892 | ## Author: TN [old] |
---|
| 893 | def getopmatrix (kind,n): |
---|
| 894 | import numpy as np |
---|
| 895 | # compute matrix of operations between files |
---|
| 896 | # the matrix is 'number of files'-square |
---|
| 897 | # 1: difference (row minus column), 2: add |
---|
| 898 | # not 0 in diag : just plot |
---|
| 899 | if n == 1: |
---|
| 900 | opm = np.eye(1) |
---|
| 901 | elif kind == 'basic': |
---|
| 902 | opm = np.eye(n) |
---|
| 903 | elif kind == 'difference1': # show differences with 1st file |
---|
| 904 | opm = np.zeros((n,n)) |
---|
| 905 | opm[0,:] = 1 |
---|
| 906 | opm[0,0] = 0 |
---|
| 907 | elif kind == 'difference2': # show differences with 1st file AND show 1st file |
---|
| 908 | opm = np.zeros((n,n)) |
---|
| 909 | opm[0,:] = 1 |
---|
| 910 | else: |
---|
| 911 | opm = np.eye(n) |
---|
| 912 | return opm |
---|
| 913 | |
---|
| 914 | ## Author: TN [old] |
---|
| 915 | def checkcoherence (nfiles,nlat,nlon,ntime): |
---|
| 916 | if (nfiles > 1): |
---|
| 917 | if (nlat > 1): |
---|
| 918 | errormess("what you asked is not possible !") |
---|
| 919 | return 1 |
---|
| 920 | |
---|
| 921 | ## Author: TN |
---|
| 922 | def readslices(saxis): |
---|
| 923 | from numpy import empty |
---|
| 924 | if saxis == None: |
---|
| 925 | zesaxis = None |
---|
| 926 | else: |
---|
| 927 | zesaxis = empty((len(saxis),2)) |
---|
| 928 | for i in range(len(saxis)): |
---|
| 929 | a = separatenames(saxis[i]) |
---|
| 930 | if len(a) == 1: |
---|
| 931 | zesaxis[i,:] = float(a[0]) |
---|
| 932 | else: |
---|
| 933 | zesaxis[i,0] = float(a[0]) |
---|
| 934 | zesaxis[i,1] = float(a[1]) |
---|
| 935 | |
---|
| 936 | return zesaxis |
---|
| 937 | |
---|
[399] | 938 | ## Author: AS |
---|
| 939 | def bidimfind(lon2d,lat2d,vlon,vlat): |
---|
| 940 | import numpy as np |
---|
| 941 | if vlat is None: array = (lon2d - vlon)**2 |
---|
| 942 | elif vlon is None: array = (lat2d - vlat)**2 |
---|
| 943 | else: array = (lon2d - vlon)**2 + (lat2d - vlat)**2 |
---|
| 944 | idy,idx = np.unravel_index( np.argmin(array), lon2d.shape ) |
---|
| 945 | if vlon is not None: |
---|
| 946 | #print lon2d[idy,idx],vlon |
---|
| 947 | if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d) |
---|
| 948 | if vlat is not None: |
---|
| 949 | #print lat2d[idy,idx],vlat |
---|
| 950 | if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d) |
---|
| 951 | return idx,idy |
---|
| 952 | |
---|
[345] | 953 | ## Author: TN |
---|
[399] | 954 | def getsindex(saxis,index,axis): |
---|
[345] | 955 | # input : all the desired slices and the good index |
---|
| 956 | # output : all indexes to be taken into account for reducing field |
---|
| 957 | import numpy as np |
---|
[425] | 958 | ### added by AS to treat the case of stime = - LT |
---|
| 959 | if saxis is not None: |
---|
| 960 | if saxis[0][0] < 0: saxis = - saxis |
---|
| 961 | ### |
---|
[349] | 962 | if ( np.array(axis).ndim == 2): |
---|
| 963 | axis = axis[:,0] |
---|
[345] | 964 | if saxis is None: |
---|
| 965 | zeindex = None |
---|
| 966 | else: |
---|
| 967 | aaa = int(np.argmin(abs(saxis[index,0] - axis))) |
---|
| 968 | bbb = int(np.argmin(abs(saxis[index,1] - axis))) |
---|
| 969 | [imin,imax] = np.sort(np.array([aaa,bbb])) |
---|
| 970 | zeindex = np.array(range(imax-imin+1))+imin |
---|
| 971 | # because -180 and 180 are the same point in longitude, |
---|
| 972 | # we get rid of one for averaging purposes. |
---|
| 973 | if axis[imin] == -180 and axis[imax] == 180: |
---|
| 974 | zeindex = zeindex[0:len(zeindex)-1] |
---|
[392] | 975 | print "INFO: whole longitude averaging asked, so last point is not taken into account." |
---|
[345] | 976 | return zeindex |
---|
| 977 | |
---|
| 978 | ## Author: TN |
---|
| 979 | def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode): |
---|
| 980 | # Purpose of define_axis is to find x and y axis scales in a smart way |
---|
| 981 | # x axis priority: 1/time 2/lon 3/lat 4/vertical |
---|
| 982 | # To be improved !!!... |
---|
| 983 | from numpy import array,swapaxes |
---|
| 984 | x = None |
---|
| 985 | y = None |
---|
| 986 | count = 0 |
---|
| 987 | what_I_plot = array(what_I_plot) |
---|
| 988 | shape = what_I_plot.shape |
---|
| 989 | if indextime is None: |
---|
[392] | 990 | print "AXIS is time" |
---|
[345] | 991 | x = time |
---|
| 992 | count = count+1 |
---|
| 993 | if indexlon is None: |
---|
[392] | 994 | print "AXIS is lon" |
---|
[345] | 995 | if count == 0: x = lon |
---|
| 996 | else: y = lon |
---|
| 997 | count = count+1 |
---|
| 998 | if indexlat is None: |
---|
[392] | 999 | print "AXIS is lat" |
---|
[345] | 1000 | if count == 0: x = lat |
---|
| 1001 | else: y = lat |
---|
| 1002 | count = count+1 |
---|
| 1003 | if indexvert is None and dim0 is 4: |
---|
[392] | 1004 | print "AXIS is vert" |
---|
[345] | 1005 | if vertmode == 0: # vertical axis is as is (GCM grid) |
---|
| 1006 | if count == 0: x=range(len(vert)) |
---|
| 1007 | else: y=range(len(vert)) |
---|
| 1008 | count = count+1 |
---|
| 1009 | else: # vertical axis is in kms |
---|
| 1010 | if count == 0: x = vert |
---|
| 1011 | else: y = vert |
---|
| 1012 | count = count+1 |
---|
| 1013 | x = array(x) |
---|
| 1014 | y = array(y) |
---|
[392] | 1015 | print "CHECK: what_I_plot.shape", what_I_plot.shape |
---|
| 1016 | print "CHECK: x.shape, y.shape", x.shape, y.shape |
---|
[345] | 1017 | if len(shape) == 1: |
---|
[350] | 1018 | if shape[0] != len(x): |
---|
[345] | 1019 | print "WARNING HERE !!!" |
---|
| 1020 | x = y |
---|
| 1021 | elif len(shape) == 2: |
---|
| 1022 | if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]: |
---|
| 1023 | what_I_plot = swapaxes(what_I_plot,0,1) |
---|
[392] | 1024 | print "INFO: swapaxes", what_I_plot.shape, shape |
---|
[345] | 1025 | return what_I_plot,x,y |
---|
[349] | 1026 | |
---|
[428] | 1027 | # Author: TN + AS |
---|
[349] | 1028 | def determineplot(slon, slat, svert, stime): |
---|
| 1029 | nlon = 1 # number of longitudinal slices -- 1 is None |
---|
| 1030 | nlat = 1 |
---|
| 1031 | nvert = 1 |
---|
| 1032 | ntime = 1 |
---|
| 1033 | nslices = 1 |
---|
| 1034 | if slon is not None: |
---|
| 1035 | nslices = nslices*len(slon) |
---|
| 1036 | nlon = len(slon) |
---|
| 1037 | if slat is not None: |
---|
| 1038 | nslices = nslices*len(slat) |
---|
| 1039 | nlat = len(slat) |
---|
| 1040 | if svert is not None: |
---|
| 1041 | nslices = nslices*len(svert) |
---|
| 1042 | nvert = len(svert) |
---|
| 1043 | if stime is not None: |
---|
| 1044 | nslices = nslices*len(stime) |
---|
| 1045 | ntime = len(stime) |
---|
| 1046 | #else: |
---|
| 1047 | # nslices = 2 |
---|
| 1048 | mapmode = 0 |
---|
| 1049 | if slon is None and slat is None: |
---|
| 1050 | mapmode = 1 # in this case we plot a map, with the given projection |
---|
| 1051 | |
---|
| 1052 | return nlon, nlat, nvert, ntime, mapmode, nslices |
---|
[440] | 1053 | |
---|
[448] | 1054 | ## Author: AC |
---|
| 1055 | ## Reduce complexity of main script by moving the contour part here. Also allow to call it from elsewhere |
---|
| 1056 | ## which can be usefull |
---|
| 1057 | # |
---|
| 1058 | #def call_contour(what_I_plot,error,x,y,m,lon,lat,vert,time,vertmode,ze_var2,indextime,indexlon,indexlat,indexvert,yintegral,mapmode,typefile,var2,ticks): |
---|
| 1059 | # from matplotlib.pyplot import contour, plot, clabel |
---|
| 1060 | # import numpy as np |
---|
| 1061 | # #what_I_plot = what_I_plot*mult |
---|
| 1062 | # if not error: |
---|
| 1063 | # if mapmode == 1: |
---|
| 1064 | # if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) |
---|
| 1065 | # zevmin, zevmax = calculate_bounds(what_I_plot) |
---|
| 1066 | # zelevels = np.linspace(zevmin,zevmax,ticks) #20) |
---|
| 1067 | # if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.) |
---|
| 1068 | # if mapmode == 0: |
---|
| 1069 | # #if typefile in ['mesoideal']: what_I_plot = dumpbdy(what_I_plot,0,stag='W') |
---|
| 1070 | # itime=indextime |
---|
| 1071 | # if len(what_I_plot.shape) is 3: itime=[0] |
---|
| 1072 | # what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ |
---|
| 1073 | # itime,what_I_plot, len(ze_var2.shape),vertmode) |
---|
| 1074 | # ### If we plot a 2-D field |
---|
| 1075 | # if len(what_I_plot.shape) is 2: |
---|
| 1076 | # #zelevels=[1.] |
---|
| 1077 | # if mapmode == 0:cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5) |
---|
| 1078 | # elif mapmode == 1:cs = m.contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5) |
---|
| 1079 | # #clabel(cs,zelevels,inline=3,fmt='%1.1f',fontsize=7) |
---|
| 1080 | # ### If we plot a 1-D field |
---|
| 1081 | # elif len(what_I_plot.shape) is 1: |
---|
| 1082 | # plot(what_I_plot,x) |
---|
| 1083 | # else: |
---|
| 1084 | # errormess("There is an error in reducing field !") |
---|
| 1085 | # return error |
---|
[440] | 1086 | |
---|