Changeset 525
- Timestamp:
- Feb 13, 2012, 12:27:22 PM (13 years ago)
- Location:
- trunk/UTIL/PYTHON
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/UTIL/PYTHON/mymath.py
r430 r525 38 38 else: 39 39 return np.array(field).mean(axis=axis) 40 41 def sum (field,axis=None): 42 import numpy as np 43 if field is None: return None 44 else: 45 if type(field).__name__=='MaskedArray': 46 field.set_fill_value(np.NaN) 47 zout=np.ma.array(field).sum(axis=axis) 48 if axis is not None: 49 zout.set_fill_value(np.NaN) 50 return zout.filled() 51 else:return zout 52 elif (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')): 53 zout=np.ma.masked_invalid(field).sum(axis=axis) 54 if axis is not None: 55 zout.set_fill_value([np.NaN]) 56 return zout.filled() 57 else:return zout 58 else: 59 return np.array(field).sum(axis=axis) 40 60 41 61 def deg (): -
trunk/UTIL/PYTHON/myplot.py
r518 r525 36 36 return ltst 37 37 38 ## Author: AS, AC 38 ## Author: AS, AC, JL 39 39 def whatkindfile (nc): 40 40 if 'controle' in nc.variables: typefile = 'gcm' 41 41 elif 'phisinit' in nc.variables: typefile = 'gcm' 42 elif 'time_counter' in nc.variables: typefile = 'earthgcm' 42 43 elif hasattr(nc,'START_DATE'): 43 44 if '9999' in getattr(nc,'START_DATE') : typefile = 'mesoideal' … … 77 78 78 79 ## Author: AS + TN + AC 79 def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None ):80 def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None,mesharea=None): 80 81 ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x" 81 82 ### it would be actually better to name d4 d3 d2 d1 as t z y x 82 83 ### ... note, anomaly is only computed over d1 and d2 for the moment 83 84 import numpy as np 84 from mymath import max,mean,min 85 from mymath import max,mean,min,sum 85 86 csmooth = 12 ## a fair amount of grid points (too high results in high computation time) 86 87 if redope is not None: … … 95 96 #elif redope == "maxx": input = max(input,axis=3) ; d4 = None 96 97 dimension = np.array(input).ndim 97 shape = np.array( input).shape98 shape = np.array(np.array(input).shape) 98 99 #print 'd1,d2,d3,d4: ',d1,d2,d3,d4 99 100 if anomaly: print 'ANOMALY ANOMALY' … … 107 108 ### now the main part 108 109 if dimension == 2: 109 if d2 >= shape[0]: error = True 110 elif d1 >= shape[1]: error = True 111 elif d1 is not None and d2 is not None: output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0) 110 if mesharea is None: mesharea=np.ones(shape) 111 if max(d2) >= shape[0]: error = True 112 elif max(d1) >= shape[1]: error = True 113 elif d1 is not None and d2 is not None: 114 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 115 output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea 112 116 elif d1 is not None: output = mean(input[:,d1],axis=1) 113 elif d2 is not None: output = mean(input[d2,:],axis=0)117 elif d2 is not None: totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea 114 118 elif dimension == 3: 119 if mesharea is None: mesharea=np.ones(shape[[1,2]]) 115 120 if max(d4) >= shape[0]: error = True 116 121 elif max(d2) >= shape[1]: error = True 117 122 elif max(d1) >= shape[2]: error = True 118 123 elif d4 is not None and d2 is not None and d1 is not None: 119 output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0) 120 elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0) 124 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 125 output = mean(input[d4,:,:],axis=0) 126 output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea 127 elif d4 is not None and d2 is not None: 128 output = mean(input[d4,:,:],axis=0) 129 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea 121 130 elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1) 122 elif d2 is not None and d1 is not None: output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1) 123 elif d1 is not None: output = mean(input[:,:,d1],axis=2) 124 elif d2 is not None: output = mean(input[:,d2,:],axis=1) 125 elif d4 is not None: output = mean(input[d4,:,:],axis=0) 131 elif d2 is not None and d1 is not None: 132 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 133 output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea 134 elif d1 is not None: output = mean(input[:,:,d1],axis=2) 135 elif d2 is not None: totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea 136 elif d4 is not None: output = mean(input[d4,:,:],axis=0) 126 137 elif dimension == 4: 138 if mesharea is None: mesharea=np.ones(shape[[2,3]]) 127 139 if max(d4) >= shape[0]: error = True 128 140 elif max(d3) >= shape[1]: error = True … … 133 145 output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) 134 146 if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 135 output = mean(output[d2,:],axis=0) 136 output = mean(output[d1],axis=0) 147 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 148 output = output*mesharea 149 output = sum(output[d2,:],axis=0) 150 output = sum(output[d1],axis=0)/totalarea 137 151 elif d4 is not None and d3 is not None and d2 is not None: 138 152 output = mean(input[d4,:,:,:],axis=0) 139 153 output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) 140 if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 141 output = mean(output[d2,:],axis=0) 154 if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 155 totalarea=sum(mesharea[d2,:],axis=0) 156 output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea 142 157 elif d4 is not None and d3 is not None and d1 is not None: 143 158 output = mean(input[d4,:,:,:],axis=0) … … 149 164 if anomaly: 150 165 for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.) 151 output = mean(output[:,d2,:],axis=1)152 output = mean(output[:,d1],axis=1) 166 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 167 output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea 153 168 #noperturb = smooth1d(output,window_len=7) 154 169 #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4] … … 158 173 if anomaly: 159 174 for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.) 160 output = mean(output[:,d2,:],axis=1)161 output = mean(output[:,d1],axis=1) 175 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 176 output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea 162 177 elif d4 is not None and d3 is not None: 163 178 output = mean(input[d4,:,:,:],axis=0) … … 166 181 elif d4 is not None and d2 is not None: 167 182 output = mean(input[d4,:,:,:],axis=0) 168 output = mean(output[:,d2,:],axis=1) 183 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea 169 184 elif d4 is not None and d1 is not None: 170 185 output = mean(input[d4,:,:,:],axis=0) … … 172 187 elif d3 is not None and d2 is not None: 173 188 output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 174 output = mean(output[:,d2,:],axis=1) 189 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea 175 190 elif d3 is not None and d1 is not None: 176 191 output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 177 192 output = mean(output[:,:,d1],axis=2) 178 193 elif d2 is not None and d1 is not None: 179 output = mean(input[:,:,d2,:],axis=2)180 output = mean(output[:,:,d1],axis=2) 194 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 195 output = output*mesharea; output = sum(output[:,:,d2,:],axis=2); output = sum(output[:,:,d1],axis=2)/totalarea 181 196 elif d1 is not None: output = mean(input[:,:,:,d1],axis=3) 182 elif d2 is not None: output = mean(input[:,:,d2,:],axis=2) 197 elif d2 is not None: 198 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,:,d2,:],axis=2)/totalarea 183 199 elif d3 is not None: output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 184 200 elif d4 is not None: output = mean(input[d4,:,:,:],axis=0) … … 427 443 elif typefile in ['gcm']: 428 444 [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) 445 elif typefile in ['earthgcm']: 446 [lon2d,lat2d] = getcoord2d(nc,nlat="lat",nlon="lon",is1d=True) 429 447 elif typefile in ['geo']: 430 448 [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') -
trunk/UTIL/PYTHON/planetoplot.py
r518 r525 9 9 ### A. Colaitis -- LMD -- 12/2011 -- Added movie capability [mencoder must be installed] 10 10 ### A. Spiga -- LMD -- 12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...] + consistent generic movie loop 11 11 ### J. Leconte -- LMD -- 02/2012 -- Added area weighted averaging. Compatibility with terrestrial gcm. 12 12 def planetoplot (namefiles,\ 13 13 level=0,\ … … 126 126 ### ... TYPEFILE 127 127 typefile = whatkindfile(nc) 128 if typefile in ['mesoideal']: mapmode=0;winds=False 129 elif typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1: mapmode=0;winds=False 128 if typefile in ['mesoideal']: mapmode=0 129 elif typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1: mapmode=0 130 elif typefile in ['earthgcm'] and len(nc.variables["lon"][:]) is 1 and len(nc.variables["lat"][:]) is 1: mapmode=0 130 131 if mapmode == 0: winds=False 131 elif mapmode == 1 :132 elif mapmode == 1 and redope is None: 132 133 if svert is None: svert = readslices(str(level)) ; nvert=1 133 134 if stime is None and mrate is None: … … 166 167 ############ LOAD 4D DIMENSIONS : x, y, z, t ############# 167 168 ########################################################## 168 if typefile == "gcm": 169 lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:] ; vert = nc.variables["altitude"][:] 170 if "Time" in nc.variables: time = nc.variables["Time"][:] 171 elif "time" in nc.variables: time = nc.variables["time"][:] 172 else: errormess("no time axis found.") 169 if typefile in ["gcm","earthgcm"]: 170 ### SPACE 171 if typefile == "gcm": lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:] ; vert = nc.variables["altitude"][:] 172 elif typefile == "earthgcm": lat = nc.variables["lat"][:] ; lon = nc.variables["lon"][:] ; vert = nc.variables["Alt"][:] 173 if "aire" in nc.variables: area = nc.variables["aire"][:,:] #JL to weight means with the area 174 else: area = None 175 ### TIME 176 if "Time" in nc.variables: time = nc.variables["Time"][:] 177 elif "time" in nc.variables: time = nc.variables["time"][:] 178 elif "time_counter" in nc.variables: time = nc.variables["time_counter"][:]/86400. #### time counter cinverstion from s-> days 179 else: errormess("no time axis found.") 173 180 if axtime in ["ls","sol"]: errormess("not supported. should not be too difficult though.") 174 181 # for 1D plots (no need for longitude computation): … … 178 185 print "LOCAL TIMES.... ", time 179 186 elif typefile in ['meso','mesoapi','geo','mesoideal']: 187 area = None ## not active for the moment 180 188 ###### STUFF TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS 181 189 ###### principle: calculate correct indices then repopulate slon and slat … … 341 349 if varname: ### what is shaded. 342 350 what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \ 343 yint=yintegral, alt=vert, anomaly=anomaly, redope=redope )351 yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area ) 344 352 if mult != 2718.: what_I_plot = what_I_plot*mult 345 353 else: what_I_plot = np.log10(what_I_plot) ; print "log plot"
Note: See TracChangeset
for help on using the changeset viewer.
