Changeset 399


Ignore:
Timestamp:
Nov 18, 2011, 8:43:59 PM (13 years ago)
Author:
aslmd
Message:

GRAPHICS: improved the selection of lon,lat with pp.py for mesoscale files

Location:
trunk/UTIL/PYTHON
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/myplot.py

    r398 r399  
    22def errormess(text,printvar=None):
    33    print text
    4     if printvar: print printvar
     4    if printvar is not None: print printvar
    55    exit()
    66    return
    … …  
    829829  return zesaxis
    830830
     831## Author: AS
     832def bidimfind(lon2d,lat2d,vlon,vlat):
     833   import numpy as np
     834   if vlat is None:    array = (lon2d - vlon)**2
     835   elif vlon is None:  array = (lat2d - vlat)**2
     836   else:               array = (lon2d - vlon)**2 + (lat2d - vlat)**2
     837   idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
     838   if vlon is not None:
     839       #print lon2d[idy,idx],vlon
     840       if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
     841   if vlat is not None:
     842       #print lat2d[idy,idx],vlat
     843       if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
     844   return idx,idy
     845
    831846## Author: TN
    832 def  getsindex(saxis,index,axis):
     847def getsindex(saxis,index,axis):
    833848# input  : all the desired slices and the good index
    834849# output : all indexes to be taken into account for reducing field
    … …  
    890905   print "CHECK: x.shape, y.shape", x.shape, y.shape
    891906   if len(shape) == 1:
    892       print shape[0]
    893907      if shape[0] != len(x):
    894908         print "WARNING HERE !!!"
    895909         x = y
    896910   elif len(shape) == 2:
    897       print shape[1], len(y), shape[0], len(x)
    898911      if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
    899912         what_I_plot = swapaxes(what_I_plot,0,1)
    900913         print "INFO: swapaxes", what_I_plot.shape, shape
    901          #x0 = x
    902          #x = y
    903          #y = x0
    904    #print "define_axis", x, y
    905914   return what_I_plot,x,y
    906915

  • trunk/UTIL/PYTHON/myscript.py

    r394 r399  
    4444    parser.add_option('--lon',          action='append',dest='slon',   type="string",  default=None, help='slices along lon. 2 comma-separated values: averaging')
    4545    parser.add_option('--vert',         action='append',dest='svert',  type="string",  default=None, help='slices along vert. 2 comma-separated values: averaging')
    46     parser.add_option('--column',         action='store_true',dest='column',           default=False,help='changes the function of --vert z1,z2 from MEAN to INTEGRATE along z')
     46    parser.add_option('--column',       action='store_true',dest='column',             default=False,help='changes --vert z1,z2 from MEAN to INTEGRATE along z')
    4747    parser.add_option('--time',         action='append',dest='stime',  type="string",  default=None, help='slices along time. 2 comma-separated values: averaging')
    4848    parser.add_option('--xmax',         action='store',dest='xmax',    type="float",   default=None, help='max value for x-axis in contour-plots [max(xaxis)]')

  • trunk/UTIL/PYTHON/planetoplot.py

    r398 r399  
    1818           var=None,\
    1919           colorb="def",\
    20            winds=True,\
     20           winds=False,\
    2121           addchar=None,\
    2222           interv=[0,1],\
    … …  
    6464                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
    6565                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
    66                        getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices
     66                       getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar
    6767    from mymath import deg,max,min,mean,get_tsat
    6868    import matplotlib as mpl
    … …  
    8686    ################################
    8787    nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
     88    vlon = None ; vlat = None
     89    if slon is not None: vlon = slon[0][0]
     90    if slat is not None: vlat = slat[0][0]
    8891    if mapmode == 0:       winds=False
    8992    elif mapmode == 1:     
    … …  
    9396    print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
    9497    print "********** MAPMODE: ", mapmode
    95     if fileref is not None:       all_var  = [[]]*(zelen+2) ;  all_varname = [[]]*(zelen+2)
    96     else:                         all_var  = [[]]*zelen     ;  all_var2  = [[]]*zelen ;  all_title = [[]]*zelen ;  all_varname = [[]]*zelen
     98    if fileref is not None:       all_var  = [[]]*(zelen+2) ;  all_varname = [[]]*(zelen+2)  ;  all_namefile = [[]]*(zelen+2)
     99    else:                         all_var  = [[]]*zelen     ;  all_var2  = [[]]*zelen ;  all_title = [[]]*zelen ;  all_varname = [[]]*zelen ; all_namefile = [[]]*zelen
    97100 
    98101    #################################################################################################
    99102    ### Loop over the files + vars initially separated by commas to be plotted on the same figure ###
    100103    #################################################################################################
    101     k = 0
    102     firstfile = True
     104    k = 0 ; firstfile = True
    103105    for nnn in range(len(namefiles)):
    104106     for vvv in range(len(var)):
    … …  
    127129      if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)
    128130      ### ... COORDINATES, could be moved below
    129       [lon2d,lat2d] = getcoorddef(nc)                       
     131      [lon2d,lat2d] = getcoorddef(nc)
    130132      ### ... PROJECTION
    131133      if proj == None:   proj = getproj(nc)                 
    … …  
    144146          vert = nc.variables["altitude"][:]
    145147      elif typefile in ['meso','mesoapi']:
     148          if vlon is not None or vlat is not None:   indices = bidimfind(lon2d,lat2d,vlon,vlat) ; print '********** INDICES: ', indices
     149          if slon is not None: slon[0][0] = indices[0] ; slon[0][1] = indices[0]
     150          if slat is not None: slat[0][0] = indices[1] ; slat[0][1] = indices[1]
    146151          if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
    147152          else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
    … …  
    173178         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
    174179         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom)
    175          ### ... NAME FOR TITLE and FILE (have to do something for multiple vars...)
    176          basename = getname(var=varname,winds=winds,anomaly=anomaly)
    177          basename = basename + getstralt(nc,level)  ## can be moved elsewhere for a more generic routine
    178180
    179181      ##### SPECIFIC
    … …  
    183185          else:                                 tinput=tt
    184186          all_var[k]=get_tsat(vert,tinput,zlon=lon,zlat=lat,zalt=vert,ztime=time)
    185           all_varname[k] = varname
    186187      else:
    187188      ##### GENERAL STUFF HERE
    188189          all_var[k] = getfield(nc,varname)
    189           all_varname[k] = varname
     190      all_varname[k] = varname
     191      all_namefile[k] = namefile
    190192      if var2: all_var2[k] = getfield(nc,var2)
    191193     
    … …  
    208210                numplot = numplot+2
    209211             elif ope in ["cat"]:
    210                 tab = all_var[0];k = 0
     212                tab = all_var[0];k = 1
    211213                while k != len(namefiles)-1:
    212214                    tab = np.append(tab,all_var[k],axis=0)
    … …  
    276278       #print slon, slat, svert, stime       
    277279       #print index_f
    278        #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
     280       #print "index lon,lat,vert,time", indexlon, indexlat, indexvert, indextime
    279281       #if varname != all_varname[index_f]:  indextime = first
    280282####################################################################
    … …  
    401403       ### Next subplot
    402404       basename = getname(var=varname,winds=winds,anomaly=anomaly)
     405       basename = basename + getstralt(nc,level) 
    403406       if typefile in ['mesoapi','meso']:
     407            if slon is not None: basename = basename + "_lon_" + str(int(lon2d[indices[1],indices[0]]))
     408            if slat is not None: basename = basename + "_lat_" + str(int(lat2d[indices[1],indices[0]]))
    404409            plottitle = basename
    405             if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
     410            if addchar: 
     411                [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] )
     412                plottitle = plottitle + addchar + "_LT"+str(ltst)
    406413            else:        plottitle = plottitle + "_LT"+str(ltst)
    407414       else:
Note: See TracChangeset for help on using the changeset viewer.