Ignore:
Timestamp:
Nov 17, 2011, 2:17:48 AM (13 years ago)
Author:
aslmd
Message:

PYTHON: quite a lot of modifs and tests. extended multivar subplot capabilities + cosmetic changes + general cleaning. corrected some stuff that were not working with mesoscale. and there is now a common script for meso and gcm svn status -uq! it is named pp.py [stands for: planetoplot]. put some examples in README file. testing is probably a bit necessary for most complex options.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/myplot.py

    r391 r392  
    4343    elif 'HGT_M' in nc.variables:    typefile = 'geo'
    4444    #else:                            errormess("whatkindfile: typefile not supported.")
    45     else: typefile = 'gcm' # for lslin-ed files from gcm
     45    else:                            typefile = 'gcm' # for lslin-ed files from gcm
    4646    return typefile
    4747
    … …  
    5050    ## this allows to get much faster (than simply referring to nc.variables[var])
    5151    dimension = len(nc.variables[var].dimensions)
    52     print "   Opening variable with", dimension, "dimensions ..."
     52    #print "   Opening variable with", dimension, "dimensions ..."
    5353    if dimension == 2:    field = nc.variables[var][:,:]
    5454    elif dimension == 3:  field = nc.variables[var][:,:,:]
    … …  
    6565    shape = np.array(input).shape
    6666    #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
    67     print 'dim,shape: ',dimension,shape
     67    print 'IN  REDUCEFIELD dim,shape: ',dimension,shape
    6868    output = input
    6969    error = False
    … …  
    9898        elif max(d1) >= shape[3]: error = True
    9999        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
    100              output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
     100             output = mean(input[d4,:,:,:],axis=0)
     101             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     102             output = mean(output[d2,:],axis=0)
     103             output = mean(output[d1],axis=0)
    101104        elif d4 is not None and d3 is not None and d2 is not None:
    102             output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[d2,:],axis=0)
     105             output = mean(input[d4,:,:,:],axis=0)
     106             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     107             output = mean(output[d2,:],axis=0)
    103108        elif d4 is not None and d3 is not None and d1 is not None:
    104             output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[:,d1],axis=1)
     109             output = mean(input[d4,:,:,:],axis=0)
     110             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     111             output = mean(output[:,d1],axis=1)
    105112        elif d4 is not None and d2 is not None and d1 is not None:
    106             output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
     113             output = mean(input[d4,:,:,:],axis=0)
     114             output = mean(output[:,d2,:],axis=1)
     115             output = mean(output[:,d1],axis=1)
    107116        elif d3 is not None and d2 is not None and d1 is not None:
    108             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
    109         elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3])
    110         elif d4 is not None and d2 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1)
    111         elif d4 is not None and d1 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,:,d1],axis=2)
    112         elif d3 is not None and d2 is not None:  output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,d2,:],axis=1)
    113         elif d3 is not None and d1 is not None:  output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,:,d1],axis=0)
    114         elif d2 is not None and d1 is not None:  output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)
    115         elif d1 is not None:                     output = mean(input[:,:,:,d1],axis=3)
    116         elif d2 is not None:                     output = mean(input[:,:,d2,:],axis=2)
    117         elif d3 is not None:                     output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt[d3])
    118         elif d4 is not None:                     output = mean(input[d4,:,:,:],axis=0)
     117             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
     118             output = mean(output[:,d2,:],axis=1)
     119             output = mean(output[:,d1],axis=1)
     120        elif d4 is not None and d3 is not None: 
     121             output = mean(input[d4,:,:,:],axis=0)
     122             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     123        elif d4 is not None and d2 is not None: 
     124             output = mean(input[d4,:,:,:],axis=0)
     125             output = mean(output[:,d2,:],axis=1)
     126        elif d4 is not None and d1 is not None: 
     127             output = mean(input[d4,:,:,:],axis=0)
     128             output = mean(output[:,:,d1],axis=2)
     129        elif d3 is not None and d2 is not None: 
     130             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
     131             output = mean(output[:,d2,:],axis=1)
     132        elif d3 is not None and d1 is not None: 
     133             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
     134             output = mean(output[:,:,d1],axis=0)
     135        elif d2 is not None and d1 is not None: 
     136             output = mean(input[:,:,d2,:],axis=2)
     137             output = mean(output[:,:,d1],axis=2)
     138        elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
     139        elif d2 is not None:        output = mean(input[:,:,d2,:],axis=2)
     140        elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     141        elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
    119142    dimension = np.array(output).ndim
    120143    shape = np.array(output).shape
    121     print 'dim,shape: ',dimension,shape
     144    print 'OUT REDUCEFIELD dim,shape: ',dimension,shape
    122145    return output, error
    123146
    124 ## Author: AC
    125 
    126 def reduce_zaxis (input,ax=None,yint=False,vert=None):
     147## Author: AC + AS
     148def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
    127149    from mymath import max,mean
    128150    from scipy import integrate
    129     if yint:
     151    if yint and vert is not None and indice is not None:
    130152       if type(input).__name__=='MaskedArray':
    131153          input.set_fill_value([np.NaN])
    132           output = integrate.trapz(input.filled(),x=vert,axis=ax)
     154          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
    133155       else:
    134           output = integrate.trapz(input.filled(),x=vert,axis=ax)
     156          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
    135157    else:
    136158       output = mean(input,axis=ax)
    … …  
    257279            if cen_lat > 10.:   
    258280                proj="npstere"
    259                 print "NP stereographic polar domain"
     281                #print "NP stereographic polar domain"
    260282            else:           
    261283                proj="spstere"
    262                 print "SP stereographic polar domain"
     284                #print "SP stereographic polar domain"
    263285        elif map_proj == 1:
    264             print "lambert projection domain"
     286            #print "lambert projection domain"
    265287            proj="lcc"
    266288        elif map_proj == 3:
    267             print "mercator projection"
     289            #print "mercator projection"
    268290            proj="merc"
    269291        else:
    … …  
    460482        elif wlat[1] >= 0.:    blat = wlat[0]
    461483        elif wlat[0] <= 0.:    blat = wlat[1]
    462     print "blat ", blat
     484    #print "blat ", blat
    463485    h = 50.  ## en km
    464486    radius = 3397200.
    … …  
    484506    #    step = np.min([5.,step])
    485507    #    steplon = step
    486     print step
    487508    m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
    488509    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
    … …  
    534555    ###
    535556    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
    536     print "field ", min(fieldcalc), max(fieldcalc)
    537     print "bounds ", zevmin, zevmax
     557    print "BOUNDS field ", min(fieldcalc), max(fieldcalc)
     558    print "BOUNDS adopted ", zevmin, zevmax
    538559    return zevmin, zevmax
    539560
    … …  
    545566    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
    546567    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
    547     print "new min ", min(what_I_plot)
     568    print "NEW MIN ", min(what_I_plot)
    548569    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
    549570    what_I_plot[ what_I_plot > zevmax ] = zevmax
    550     print "new max ", max(what_I_plot)
    551    
     571    print "NEW MAX ", max(what_I_plot)
    552572    return what_I_plot
    553573
    … …  
    556576    from mymath import max,min
    557577    lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
    558     print "on vire en dessous de ", lim
     578    print "NO PLOT BELOW VALUE ", lim
    559579    what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
    560580    return what_I_plot
    … …  
    566586    [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
    567587                    [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
    568     print "zoom %",zoom,wlon,wlat
     588    print "ZOOM %",zoom,wlon,wlat
    569589    return wlon,wlat
    570590
    … …  
    642662#W --> spectral ou jet
    643663#spectral BrBG RdBu_r
    644     print "predefined colorbars"
     664    #print "predefined colorbars"
    645665    if whichone not in whichcolorb:
    646666        whichone = "def"
    … …  
    813833      if axis[imin] == -180 and axis[imax] == 180:
    814834         zeindex = zeindex[0:len(zeindex)-1]
    815          print "whole longitude averaging asked, so last point is not taken into account."
     835         print "INFO: whole longitude averaging asked, so last point is not taken into account."
    816836  return zeindex
    817837     
    … …  
    828848   shape = what_I_plot.shape
    829849   if indextime is None:
    830       print "axis is time"
     850      print "AXIS is time"
    831851      x = time
    832852      count = count+1
    833853   if indexlon is None:
    834       print "axis is lon"
     854      print "AXIS is lon"
    835855      if count == 0: x = lon
    836856      else: y = lon
    837857      count = count+1
    838858   if indexlat is None:
    839       print "axis is lat"
     859      print "AXIS is lat"
    840860      if count == 0: x = lat
    841861      else: y = lat
    842862      count = count+1
    843863   if indexvert is None and dim0 is 4:
    844       print "axis is vert"
     864      print "AXIS is vert"
    845865      if vertmode == 0: # vertical axis is as is (GCM grid)
    846866         if count == 0: x=range(len(vert))
    … …  
    853873   x = array(x)
    854874   y = array(y)
    855    print "what_I_plot.shape", what_I_plot.shape
    856    print "x.shape, y.shape", x.shape, y.shape
     875   print "CHECK: what_I_plot.shape", what_I_plot.shape
     876   print "CHECK: x.shape, y.shape", x.shape, y.shape
    857877   if len(shape) == 1:
    858878      print shape[0]
    … …  
    864884      if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
    865885         what_I_plot = swapaxes(what_I_plot,0,1)
    866          print "swapaxes", what_I_plot.shape, shape
     886         print "INFO: swapaxes", what_I_plot.shape, shape
    867887         #x0 = x
    868888         #x = y
    … …  
    896916    if slon is None and slat is None:
    897917       mapmode = 1 # in this case we plot a map, with the given projection
    898     #elif proj is not None:
    899     #   print "WARNING: you specified a", proj,\
    900     #   "projection but asked for slices", slon,"in longitude and",slat,"in latitude"
    901     print "mapmode: ", mapmode
    902918
    903919    return nlon, nlat, nvert, ntime, mapmode, nslices
Note: See TracChangeset for help on using the changeset viewer.