Changeset 725


Ignore:
Timestamp:
Jul 17, 2012, 12:17:17 AM (12 years ago)
Author:
aslmd
Message:

UTIL PYTHON : improved comments. otherwise purely cosmetic changes to improve lisibility.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/planetoplot.py

    r724 r725  
    1111### J. Leconte   -- LMD --    02/2012 -- Added area weighted averaging. Compatibility with terrestrial gcm.
    1212### A. Spiga     -- LMD --    03/2012 -- Cleaning and improved comments
     13### T. Navarro   -- LMD --    04/2012 -- Added capabilities (e.g. histograms for difference maps)
     14### A. Colaitis  -- LMD -- 05~06/2012 -- Added capabilities for analysis of mesoscale files (wind speed, etc...)
     15### A. Spiga     -- LMD -- 04~07/2012 -- Added larger support of files + planets. Corrected a few bugs. Cleaning and improved comments
     16
    1317def planetoplot (namefiles,\
    1418           level=0,\
    … …  
    8387    from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
    8488    import matplotlib as mpl
    85     from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, subplots_adjust, axes, clabel
     89    from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, \
     90                                  pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, \
     91                                  subplots_adjust, axes, clabel
    8692    from matplotlib.cm import get_cmap
    8793    #from mpl_toolkits.basemap import cm
    … …  
    156162             print "WELL... nothing about time axis. I took default: first time reference stored in file."
    157163    ### we get the names of variables to be read. in case only one available, we choose this one.
     164    ### (we take out of the test specific names e.g. UV is not in the file but used to ask a wind speed computation)
    158165      varname=var[vvv]
    159       if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT'])) and not ((typefile in ['gcm']) and (varname in ['enfact']))):
     166      logical_novarname = varname not in nc.variables
     167      specificname_meso = ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT']
     168      logical_nospecificname_meso = not ((typefile in ['meso']) and (varname in specificname_meso))
     169      specificname_gcm = ['enfact']
     170      logical_nospecificname_gcm = not ((typefile in ['gcm']) and (varname in specificname_gcm))
     171      if ( logical_novarname and logical_nospecificname_meso and logical_nospecificname_gcm ):
    160172          if len(varinfile) == 1:   varname = varinfile[0]
    161173          else:                     varname = False
    … …  
    178190         elif zarea is not None:           [wlon,wlat] = latinterv(area=zarea)
    179191
    180 ##########################################################
    181 ############ LOAD 4D DIMENSIONS : x, y, z, t #############
    182 ##########################################################
     192#############################################################
     193############ WE LOAD 4D DIMENSIONS : x, y, z, t #############
     194#############################################################
     195
     196    ### TYPE 1 : GCM files or simple files
    183197      if typefile in ["gcm","earthgcm","ecmwf"]:
     198      ### this is needed for continuity
    184199          if slon is not None: sslon = slon 
    185200          if slat is not None: sslat = slat 
    186           ### LAT : done in getcoorddef
    187           lat = lat2d[0,:]
    188           ### LON : done in getcoorddef
    189           lon = lon2d[:,0]
    190           ### ALT
     201      ### we define lat/lon vectors. we get what was done in getcoorddef.
     202          lat = lat2d[0,:] ; lon = lon2d[:,0]
     203      ### we define areas. this is needed for calculate means and weight with area. this is not compulsory (see reduce_field).
     204          if "aire" in nc.variables:      area = nc.variables["aire"][:,:]
     205          ### --> add a line here if your reference is not present
     206          else:                           area = None
     207      ### we define altitude vector. either it is referenced or it is guessed based on last variable's dimensions.
    191208          if "altitude" in nc.variables:   vert = nc.variables["altitude"][:]
    192209          elif "Alt" in nc.variables:      vert = nc.variables["Alt"][:]
    193210          elif "lev" in nc.variables:      vert = nc.variables["lev"][:]
    194211          elif "presnivs" in nc.variables: vert = nc.variables["presnivs"][:]
     212          ### --> add a line here if your reference is not present
    195213          else:
    196               print "No altitude found. Try to build a simple axis."
    197               dadim = getdimfromvar(nc)
     214              print "No altitude found. Try to build a simple axis." ; dadim = getdimfromvar(nc)
    198215              if   len(dadim) == 4:  print "-- 4D field. Assume z is dim 2." ; vert = np.arange(dadim[-3])
    199216              elif len(dadim) == 3:  print "-- 3D field. Assume z is dim 1." ; vert = [0.]
    200               else:                  errormess("-- 2D or 1D field. Unclear. Stop.")
    201           #else:                            vert = [0.]
    202           ### AIRE (to weight means with the area)
    203           if "aire" in nc.variables:      area = nc.variables["aire"][:,:] 
    204           else:                           area = None
    205           ### TIME
     217              else:                  vert = [0.]
     218      ### we define time vector. either it is referenced or it is guessed based on last variable's dimensions.
    206219          if "Time" in nc.variables:            time = nc.variables["Time"][:]
    207220          elif "time" in nc.variables:          time = nc.variables["time"][:]
    208           elif "time_counter" in nc.variables:  time = nc.variables["time_counter"][:]/86400. #### time counter cinverstion from s-> days
     221          elif "time_counter" in nc.variables:  time = nc.variables["time_counter"][:]/86400. #### convert from s to days
     222          ### --> add a line here if your reference is not present
    209223          else:
    210224              print "No time found. Try to build a simple axis. Assume t is dim 1."
    … …  
    212226              if   len(dadim) == 4:  time = np.arange(dadim[-4])
    213227              elif len(dadim) == 3:  time = np.arange(dadim[-3])
    214               else:                  errormess("-- I am sorry this is 2D or 1D field. I failed.")
    215           #else:                                time = [0.] #errormess("no time axis found.")
     228              else:                  time = [0.] #errormess("no time axis found.")
     229          ### (SPECIFIC. convert to Ls for Martian GCM files.)
    216230          if axtime in ["ls"]:
    217231              print "converting to Ls ..."
    … …  
    219233                time[iii] = sol2ls(time[iii])
    220234                if iii > 0:
    221                   while abs(time[iii]-time[iii-1]) > 300:
    222                     time[iii] = time[iii]+360
    223           # for 1D plots (no need for longitude computation):
     235                  while abs(time[iii]-time[iii-1]) > 300: time[iii] = time[iii]+360
     236          ### (a case where the user would like to set local times. e.g. : 1D plot with no longitude reference)
    224237          if axtime in ["lt"]:
    225238              if initime == -1: initime=input("Please type initial local time:")
    226               time = (initime+time*24)%24
    227               print "LOCAL TIMES.... ", time
     239              time = (initime+time*24)%24 ; print "LOCAL TIMES.... ", time
     240
     241    ### TYPE 2 : MESOSCALE FILES
    228242      elif typefile in ['meso','geo']:
    229           area = None ## not active for the moment
    230           ###### STUFF TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
    231           ###### principle: calculate correct indices then repopulate slon and slat
     243          ### area not active with mesoscale files
     244          area = None
     245          ### HACK TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
     246          ### principle: calculate correct indices then repopulate slon and slat
    232247          if slon is not None or slat is not None:
    233               if firstfile and save == 'png' and typefile == 'meso' and "HGT" in varinfile and display:   iwantawhereplot = nc     #show a topo map with a cross on the chosen point
    234               else:                                                                           iwantawhereplot = None   #do not show anything, just select indices
     248              show_topo_map_user = ((save == 'png') and (typefile == 'meso') and ("HGT" in varinfile) and display)
     249              if firstfile and show_topo_map_user:  iwantawhereplot = nc     #show a topo map with a cross on the chosen point
     250              else:                                 iwantawhereplot = None   #do not show anything, just select indices
    235251              numlon = 1 ; numlat = 1
    236252              if slon is not None:   numlon = slon.shape[1]   
    … …  
    249265                 if slat is not None: sslat[0][jjj] = indices[0,jjj,0] #...this is idy
    250266              lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] )
    251           ######
    252           if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  ### important to do that now and not before
    253           ######
     267          ### we get rid of boundary relaxation zone for plots. important to do that now and not before.
     268          if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6) 
     269          ### we read the keyword for vertical dimension. we take care for vertical staggering.
    254270          if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
    255271          else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
    256272          if (var2 is not None and var2 not in ['PHTOT','W']): dumped_vert_stag=True ; vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
    257273          else:                                                dumped_vert_stag=False
     274          ### we read the keyword for horizontal dimensions. we take care for horizontal staggering.
    258275          if varname in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
    259276          else:                 latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
    … …  
    261278          else:                 londim='WEST-EAST_PATCH_END_UNSTAG'
    262279          lon = np.arange(0,getattr(nc,londim),1) ; lat = np.arange(0,getattr(nc,latdim),1)
    263           ###
     280          ### we define the time axis and take care of various specificities (lt, ls, sol) or issues (concatenation)
    264281          if axtime in ["ls","sol"]:
    265282              lstab, soltab, lttab = getlschar ( namefile, getaxis = True )
    … …  
    279296              time=check_localtime(time)
    280297              print "LOCAL TIMES.... ", time
    281 
    282           ###
     298          ### we define the vertical axis (or lack thereof) and cover possibilities for it to be altitude, pressure, geopotential. quite SPECIFIC.
    283299          if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
    284300          else:
    … …  
    290306              elif vertmode == 1 or vertmode == 2:  vert = nc.variables["vert"][:]        ## pressure in Pa
    291307              else:                                 vert = nc.variables["vert"][:]/1000.  ## altitude in km
    292        #if firstfile:
    293        #   lat0 = lat
    294        #elif len(lat0) != len(lat):
    295        #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
    296        #elif sum((lat == lat0) == False) != 0:
    297        #   errormess("Not the same latitudes !", [lat,lat0])
    298        ## Faire d'autre checks sur les compatibilites entre fichiers!!
    299 ##########################################################
    300 ##########################################################
    301 ##########################################################
    302 
     308####################################################################
     309############ END of WE LOAD 4D DIMENSIONS : x, y, z, t #############
     310####################################################################
     311
     312    ### we fill the arrays of varname, namefile, time at the current step considered (NB: why use both k and nnn ?)
    303313      all_varname[k] = varname
    304314      all_namefile[k] = namefile
    … …  
    306316      if var2:
    307317         all_var2[k] = getfield(nc,var2)
     318         ### v--- too SPECIFIC (see below)
    308319         if ((var2 in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (var2 not in nc.variables)): all_var2[k] = slopeamplitude(nc)
    309320      if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
    310       #####################
    311       ##### GETFIELDS #####
    312       #####################
    313       ##### SPECIFIC CASES
    314       ##### note : we could probably call those via a "toolbox" in myplot.
    315       ##### 1. saturation temperature
     321    ### we fill the arrays of fields to be plotted at the current step considered
     322      ### ----------- SPECIFIC CASES. NOT HAPPY WITH THIS. note : we could probably call those via a "toolbox" in myplot.
     323      ### ----------- 1. saturation temperature
    316324      if varname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and tsat:
    317325          tt=getfield(nc,varname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
    … …  
    319327          else:                                 tinput=tt
    320328          all_var[k]=get_tsat(vert,tinput,zlon=lon,zlat=lat,zalt=vert,ztime=time)
    321       ##### 2. wind amplitude
     329      ### ----------- 2. wind amplitude
    322330      elif ((varname in ['UV','uv','uvmet']) and (typefile in ['meso']) and (varname not in nc.variables)):
    323331          all_var[k]=windamplitude(nc)
    324332      elif ((varname in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (varname not in nc.variables)):
    325333          all_var[k]=slopeamplitude(nc)
    326       #### 3. Near surface instability
     334      ### ------------ 3. Near surface instability
    327335      elif ((varname in ['DELTAT','deltat']) and (typefile in ['meso']) and (varname not in nc.variables)):
    328336          all_var[k]=deltat0t1(nc)
    329       #### 4. Enrichment factor
     337      ### ------------ 4. Enrichment factor
    330338      elif ((typefile in ['gcm']) and (varname in ['enfact'])):
    331339          all_var[k]=enrichment_factor(nc)
    332340      else:
    333       ##### GENERAL STUFF HERE
     341      ### ideally only this line should be here
    334342          all_var[k] = getfield(nc,varname)
     343    ### we inform the user about the loop then increment the loop. this is the last line of "for namefile in namefiles"
    335344      print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
    336       #### End of for namefile in namefiles
     345
     346    ### note to self : stopped comment rewriting here.
    337347
    338348    ##################################
Note: See TracChangeset for help on using the changeset viewer.