Ignore:
Timestamp:
Nov 6, 2011, 9:28:12 PM (13 years ago)
Author:
aslmd
Message:

PYTHON: improvements on handling 1D 2D slices for gcm and meso. a few minor corrections added for the scripts to work better. it might be possible at some point to merge meso.py and gcm.py, wait and see.

Location:
trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py

    r349 r350  
    2626    parser.add_option('-s', '--stride', action='store',dest='stride',    type="int",     default=3,     help='stride vectors (def=3)')
    2727    parser.add_option('-v', '--var',    action='append',dest='var',      type="string",  default=None,  help='variable color-shaded (append)')
    28     parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=2,     help='plot number (def=2)(<0: plot LT -*numplot*)')
     28    parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=None,  help='plot number (def=2)(<0: plot LT -*numplot*)')
    2929    parser.add_option('-i', '--interp', action='store',dest='interp',    type="int",     default=None,  help='interpolation (2: p, 3: z-amr, 4:z-als)')
    3030    parser.add_option('-c', '--color',  action='store',dest='colorb',    type="string",  default="def", help='change colormap (nobar: no colorbar)')
    … …  
    138138        #############
    139139        ### Main call
    140         name = planetoplot (zenamefiles,int(zelevel),\
     140        name = planetoplot (zenamefiles,nvert=int(zelevel),vertmode=opt.interp,\
    141141                proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\
    142142                numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py

    r349 r350  
    2626    parser.add_option('-s', '--stride', action='store',dest='stride',    type="int",     default=3,     help='stride vectors (def=3)')
    2727    parser.add_option('-v', '--var',    action='append',dest='var',      type="string",  default=None,  help='variable color-shaded (append)')
    28     parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=2,     help='plot number (def=2)(<0: plot LT -*numplot*)')
     28    parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=None,  help='plot number (def=2)(<0: plot LT -*numplot*)')
    2929    parser.add_option('-i', '--interp', action='store',dest='interp',    type="int",     default=None,  help='interpolation (2: p, 3: z-amr, 4:z-als)')
    3030    parser.add_option('-c', '--color',  action='store',dest='colorb',    type="string",  default="def", help='change colormap (nobar: no colorbar)')
    … …  
    133133            #############
    134134            ### Main call
    135             name = planetoplot (zefile,int(zelevel),\
     135            name = planetoplot (zefile,nvert=int(zelevel),vertmode=opt.interp,\
    136136                proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\
    137137                numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py

    r349 r350  
    6868    output = input
    6969    error = False
    70     ### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
     70    #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
    7171    if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
    7272    if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
    … …  
    7777        if   d2 >= shape[0]: error = True
    7878        elif d1 >= shape[1]: error = True
    79         elif d1 is not None and d2 is not None:  output = input[d2,d1]
    80         elif d1 is not None:         output = input[:,d1]
    81         elif d2 is not None:         output = input[d2,:]
     79        elif d1 is not None and d2 is not None:  output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0)
     80        elif d1 is not None:         output = mean(input[:,d1],axis=1)
     81        elif d2 is not None:         output = mean(input[d2,:],axis=0)
    8282    elif dimension == 3:
    8383        if   max(d4) >= shape[0]: error = True
    8484        elif max(d2) >= shape[1]: error = True
    8585        elif max(d1) >= shape[2]: error = True
    86         elif d4 is not None and d2 is not None and d1 is not None:  output = input[d4,d2,d1]
     86        elif d4 is not None and d2 is not None and d1 is not None: 
     87            output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
    8788        elif d4 is not None and d2 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
    8889        elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
    … …  
    9697        elif max(d2) >= shape[2]: error = True
    9798        elif max(d1) >= shape[3]: error = True
    98         elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:  output = input[d4,d3,d2,d1]
    99         elif d4 is not None and d3 is not None and d2 is not None:         output = input[d4,d3,d2,:]
    100         elif d4 is not None and d3 is not None and d1 is not None:         output = input[d4,d3,:,d1]
    101         elif d4 is not None and d2 is not None and d1 is not None:         output = input[d4,:,d2,d1]
    102         elif d3 is not None and d2 is not None and d1 is not None:         output = input[:,d3,d2,d1]
     99        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 = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
     101        elif d4 is not None and d3 is not None and d2 is not None:
     102            output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0)
     103        elif d4 is not None and d3 is not None and d1 is not None:
     104            output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[:,d1],axis=1)
     105        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)
     107        elif d3 is not None and d2 is not None and d1 is not None:
     108            output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
    103109        elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0)
    104110        elif d4 is not None and d2 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1)
    … …  
    107113        elif d3 is not None and d1 is not None:  output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0)
    108114        elif d2 is not None and d1 is not None:  output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)
    109         elif d1 is not None:                       output = input[:,:,:,d1]
    110         elif d2 is not None:                       output = input[:,:,d2,:]
    111         elif d3 is not None:                       output = input[:,d3,:,:]
    112         elif d4 is not None:                       output = input[d4,:,:,:]
     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 = mean(input[:,d3,:,:],axis=1)
     118        elif d4 is not None:                     output = mean(input[d4,:,:,:],axis=0)
    113119    dimension = np.array(output).ndim
    114120    shape = np.array(output).shape
    … …  
    801807   shape = what_I_plot.shape
    802808   if indextime is None:
     809      print "axis is time"
    803810      x = time
    804811      count = count+1
    805812   if indexlon is None:
     813      print "axis is lon"
    806814      if count == 0: x = lon
    807815      else: y = lon
    808816      count = count+1
    809817   if indexlat is None:
     818      print "axis is lat"
    810819      if count == 0: x = lat
    811820      else: y = lat
    812821      count = count+1
    813822   if indexvert is None and dim0 is 4:
     823      print "axis is vert"
    814824      if vertmode == 0: # vertical axis is as is (GCM grid)
    815825         if count == 0: x=range(len(vert))
    … …  
    822832   x = array(x)
    823833   y = array(y)
     834   print "what_I_plot.shape", what_I_plot.shape
     835   print "x.shape, y.shape", x.shape, y.shape
    824836   if len(shape) == 1:
    825       if shape != len(x):
     837      print shape[0]
     838      if shape[0] != len(x):
    826839         print "WARNING HERE !!!"
    827840         x = y

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py

    r349 r350  
    77
    88def planetoplot (namefiles,\
    9            vertmode,\
     9           nvert=0,\
     10           vertmode=0,\
    1011           proj=None,\
    1112           back=None,\
    1213           target=None,
    1314           stride=3,\
    14            numplot=2,\
     15           numplot=None,\
    1516           var=None,\
    1617           colorb="def",\
    … …  
    5960    if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
    6061    zelen = len(namefiles)
    61     numplot = zelen*nslices
     62    if numplot == None:  numplot = zelen*nslices
    6263    print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot
    6364    all_var   = [[]]*zelen
    … …  
    9798##########################################################
    9899      if typefile == "gcm":
    99           lat = nc.variables["latitude"][:]
     100          lat = nc.variables["latitude"][:] 
    100101          lon = nc.variables["longitude"][:]
    101102          time = nc.variables["Time"][:]
    102103          vert = nc.variables["altitude"][:]
    103           #if firstfile:
    104           #   lat0 = lat
    105           #elif len(lat0) != len(lat):
    106           #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
    107           #elif sum((lat == lat0) == False) != 0:
    108           #   errormess("Not the same latitudes !", [lat,lat0])
    109           ## Faire d'autre checks sur les compatibilites entre fichiers!!
     104      elif typefile in ['meso','mesoapi']:
     105          if mapmode == 0:
     106              if var in ['PHTOT','W']:  vertdim='BOTTOM-TOP_PATCH_END_STAG'
     107              else:                     vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
     108              if var in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
     109              else:             latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
     110              if var in ['U']:  londim='WEST-EAST_PATCH_END_STAG'
     111              else:             londim='WEST-EAST_PATCH_END_UNSTAG'
     112              lon = np.arange(0,getattr(nc,londim),1)
     113              lat = np.arange(0,getattr(nc,latdim),1)
     114              vert = np.arange(0,getattr(nc,vertdim),1)
     115              time = np.arange(0,len(nc.variables["Times"]),1)
     116       #if firstfile:
     117       #   lat0 = lat
     118       #elif len(lat0) != len(lat):
     119       #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
     120       #elif sum((lat == lat0) == False) != 0:
     121       #   errormess("Not the same latitudes !", [lat,lat0])
     122       ## Faire d'autre checks sur les compatibilites entre fichiers!!
    110123##########################################################
    111124
    … …  
    122135         ### Name for title and graphics save file
    123136         basename = getname(var=var,winds=winds,anomaly=anomaly)
    124          basename = basename + getstralt(nc,vertmode)  ## can be moved elsewhere for a more generic routine
     137         basename = basename + getstralt(nc,nvert)  ## can be moved elsewhere for a more generic routine
    125138     
    126139      print "var, var2: ", var, var2
    … …  
    179192
    180193####################################################################
    181        if typefile in ['mesoapi','meso']:
    182            indexlon = None
    183            indexlat = None
    184            indexvert = vertmode 
    185            indextime = itime
    186            nlon = 1
    187            nlat = 1
    188            nvert = 1
    189            ntime = 1 
    190            index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
    191            print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
    192        elif typefile in ['gcm']:
     194       if typefile in ['meso','mesoapi'] and mapmode == 1:
     195               indextime = itime
     196               indexlon = None
     197               indexlat = None
     198               indexvert = nvert  ## ou svert ???
     199               nlon = 1
     200               nlat = 1
     201               nvert = 1
     202               ntime = 1
     203       else:
    193204           ## get all indexes to be taken into account for this subplot and then reduce field
    194205           ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
    195            #print "(nplot-1)%nlon", (nplot-1)%nlon
    196206           indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
    197            print "indexlon", indexlon
    198207           indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
    199208           indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) 
    200209           indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) 
    201            index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
    202            #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
     210       index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
     211       print nlon, nlat, nvert, ntime 
     212       print slon, slat, svert, stime       
     213       print index_f
     214       print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
    203215####################################################################
    204216
    … …  
    208220           #what_I_plot = what_I_plot*mult
    209221           if not error:
    210               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
     222              if mapmode == 1:
     223                  if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
    211224              zevmin, zevmax = calculate_bounds(what_I_plot)
    212225              zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20)
    … …  
    214227              if mapmode == 0:
    215228                  x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    216                         indextime,what_I_plot.shape, len(all_var[index_f].shape),vertmode)
     229                        indextime,what_I_plot.shape, len(all_var2[index_f].shape),vertmode)
    217230              ### If we plot a 2-D field
    218231              if len(what_I_plot.shape) is 2:
    … …  
    221234              ### If we plot a 1-D field
    222235              elif len(what_I_plot.shape) is 1:
    223                   plot(x,y,what_I_plot, zelevels, colors='g', linewidths = 0.33 ) #colors='w' )# , alpha=0.5)
     236                  plot(what_I_plot,x)
    224237           else:
    225238              errormess("There is an error in reducing field !")
    … …  
    239252               #    fvar = str(mult) + "*" + var
    240253               ###
    241                if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
     254               if mapmode == 1:
     255                   if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
     256               elif mapmode == 0:
     257                   what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
     258                         indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
    242259               zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
    243260               if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar.upper()))
    244261               else:                           palette = get_cmap(name=colorb)
    245                ### If we plot a 2-D field
     262               ##### 2D field
    246263               if len(what_I_plot.shape) is 2:
    247                  if mapmode == 0:
    248                      what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    249                             indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
    250264                 if not tile:
    251265                     if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
    252266                     #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
    253                      print "what_I_plot.shape", what_I_plot.shape
    254                      print "x.shape, y.shape", x.shape, y.shape
    255                      zelevels = np.linspace(zevmin,zevmax,num=ndiv+1)
     267                     zelevels = np.linspace(zevmin,zevmax)#,num=ndiv+1)
    256268                     #contourf(what_I_plot, zelevels, cmap = palette )
    257                      contourf(x,y,what_I_plot, zelevels, cmap = palette )
     269                     contourf( x, y, what_I_plot, zelevels, cmap = palette )
    258270                 else:
    259271                     if hole:  what_I_plot = nolow(what_I_plot)
    260                      pcolor(x,y,what_I_plot, cmap = palette, \
    261                            vmin=zevmin, vmax=zevmax )
    262                  if colorb != 'nobar' and var != 'HGT' :
    263                          #subplot(111)
    264                          #colorbar()         
     272                     pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
     273                 if colorb != 'nobar' and var != 'HGT' :       
    265274                         colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
    266275                                           ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
    267276                                           extend='neither',spacing='proportional')
    268277                                           # both min max neither
    269                ### If we plot a 1-D field
     278               ##### 1D field
    270279               elif len(what_I_plot.shape) is 1:
    271                  what_I_plot,x,y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    272                        indextime,what_I_plot,len(all_var[index_f].shape),vertmode)
    273                 # x = np.array(x)
    274                 # y = np.array(y)
    275                  print "what_I_plot.shape", what_I_plot.shape
    276                  print "x.shape", x.shape
    277280                 plot(x,what_I_plot)
    278                ### If we plot something that is not a 2-D or 1-D field
    279                ### (maybe plot 3-D field one day or movie ??)
     281               #### Other cases: (maybe plot 3-D field one day or movie ??)
    280282               else:
    281283                 print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!!"
Note: See TracChangeset for help on using the changeset viewer.