Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py (revision 350) @@ -26,5 +26,5 @@ parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') - parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -*numplot*)') + parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=None, help='plot number (def=2)(<0: plot LT -*numplot*)') parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)') parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') @@ -138,5 +138,5 @@ ############# ### Main call - name = planetoplot (zenamefiles,int(zelevel),\ + name = planetoplot (zenamefiles,nvert=int(zelevel),vertmode=opt.interp,\ proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\ numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py (revision 350) @@ -26,5 +26,5 @@ parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') - parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -*numplot*)') + parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=None, help='plot number (def=2)(<0: plot LT -*numplot*)') parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)') parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') @@ -133,5 +133,5 @@ ############# ### Main call - name = planetoplot (zefile,int(zelevel),\ + name = planetoplot (zefile,nvert=int(zelevel),vertmode=opt.interp,\ proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\ numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 350) @@ -68,5 +68,5 @@ output = input error = False - ### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays + #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4] if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3] @@ -77,12 +77,13 @@ if d2 >= shape[0]: error = True elif d1 >= shape[1]: error = True - elif d1 is not None and d2 is not None: output = input[d2,d1] - elif d1 is not None: output = input[:,d1] - elif d2 is not None: output = input[d2,:] + elif d1 is not None and d2 is not None: output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0) + elif d1 is not None: output = mean(input[:,d1],axis=1) + elif d2 is not None: output = mean(input[d2,:],axis=0) elif dimension == 3: if max(d4) >= shape[0]: error = True elif max(d2) >= shape[1]: error = True elif max(d1) >= shape[2]: error = True - elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,d2,d1] + elif d4 is not None and d2 is not None and d1 is not None: + output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0) elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0) elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1) @@ -96,9 +97,14 @@ elif max(d2) >= shape[2]: error = True elif max(d1) >= shape[3]: error = True - 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] - elif d4 is not None and d3 is not None and d2 is not None: output = input[d4,d3,d2,:] - elif d4 is not None and d3 is not None and d1 is not None: output = input[d4,d3,:,d1] - elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,:,d2,d1] - elif d3 is not None and d2 is not None and d1 is not None: output = input[:,d3,d2,d1] + elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: + output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0) + elif d4 is not None and d3 is not None and d2 is not None: + output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0) + elif d4 is not None and d3 is not None and d1 is not None: + output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[:,d1],axis=1) + elif d4 is not None and d2 is not None and d1 is not None: + output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1) + elif d3 is not None and d2 is not None and d1 is not None: + output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1) elif d4 is not None and d3 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0) elif d4 is not None and d2 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1) @@ -107,8 +113,8 @@ elif d3 is not None and d1 is not None: output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0) elif d2 is not None and d1 is not None: output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2) - elif d1 is not None: output = input[:,:,:,d1] - elif d2 is not None: output = input[:,:,d2,:] - elif d3 is not None: output = input[:,d3,:,:] - elif d4 is not None: output = input[d4,:,:,:] + elif d1 is not None: output = mean(input[:,:,:,d1],axis=3) + elif d2 is not None: output = mean(input[:,:,d2,:],axis=2) + elif d3 is not None: output = mean(input[:,d3,:,:],axis=1) + elif d4 is not None: output = mean(input[d4,:,:,:],axis=0) dimension = np.array(output).ndim shape = np.array(output).shape @@ -801,15 +807,19 @@ shape = what_I_plot.shape if indextime is None: + print "axis is time" x = time count = count+1 if indexlon is None: + print "axis is lon" if count == 0: x = lon else: y = lon count = count+1 if indexlat is None: + print "axis is lat" if count == 0: x = lat else: y = lat count = count+1 if indexvert is None and dim0 is 4: + print "axis is vert" if vertmode == 0: # vertical axis is as is (GCM grid) if count == 0: x=range(len(vert)) @@ -822,6 +832,9 @@ x = array(x) y = array(y) + print "what_I_plot.shape", what_I_plot.shape + print "x.shape, y.shape", x.shape, y.shape if len(shape) == 1: - if shape != len(x): + print shape[0] + if shape[0] != len(x): print "WARNING HERE !!!" x = y Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py (revision 350) @@ -7,10 +7,11 @@ def planetoplot (namefiles,\ - vertmode,\ + nvert=0,\ + vertmode=0,\ proj=None,\ back=None,\ target=None, stride=3,\ - numplot=2,\ + numplot=None,\ var=None,\ colorb="def",\ @@ -59,5 +60,5 @@ if not isinstance(namefiles, np.ndarray): namefiles = [namefiles] zelen = len(namefiles) - numplot = zelen*nslices + if numplot == None: numplot = zelen*nslices print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot all_var = [[]]*zelen @@ -97,15 +98,27 @@ ########################################################## if typefile == "gcm": - lat = nc.variables["latitude"][:] + lat = nc.variables["latitude"][:] lon = nc.variables["longitude"][:] time = nc.variables["Time"][:] vert = nc.variables["altitude"][:] - #if firstfile: - # lat0 = lat - #elif len(lat0) != len(lat): - # errormess("Not the same latitude lengths !", [len(lat0), len(lat)]) - #elif sum((lat == lat0) == False) != 0: - # errormess("Not the same latitudes !", [lat,lat0]) - ## Faire d'autre checks sur les compatibilites entre fichiers!! + elif typefile in ['meso','mesoapi']: + if mapmode == 0: + if var in ['PHTOT','W']: vertdim='BOTTOM-TOP_PATCH_END_STAG' + else: vertdim='BOTTOM-TOP_PATCH_END_UNSTAG' + if var in ['V']: latdim='SOUTH-NORTH_PATCH_END_STAG' + else: latdim='SOUTH-NORTH_PATCH_END_UNSTAG' + if var in ['U']: londim='WEST-EAST_PATCH_END_STAG' + else: londim='WEST-EAST_PATCH_END_UNSTAG' + lon = np.arange(0,getattr(nc,londim),1) + lat = np.arange(0,getattr(nc,latdim),1) + vert = np.arange(0,getattr(nc,vertdim),1) + time = np.arange(0,len(nc.variables["Times"]),1) + #if firstfile: + # lat0 = lat + #elif len(lat0) != len(lat): + # errormess("Not the same latitude lengths !", [len(lat0), len(lat)]) + #elif sum((lat == lat0) == False) != 0: + # errormess("Not the same latitudes !", [lat,lat0]) + ## Faire d'autre checks sur les compatibilites entre fichiers!! ########################################################## @@ -122,5 +135,5 @@ ### Name for title and graphics save file basename = getname(var=var,winds=winds,anomaly=anomaly) - basename = basename + getstralt(nc,vertmode) ## can be moved elsewhere for a more generic routine + basename = basename + getstralt(nc,nvert) ## can be moved elsewhere for a more generic routine print "var, var2: ", var, var2 @@ -179,26 +192,25 @@ #################################################################### - if typefile in ['mesoapi','meso']: - indexlon = None - indexlat = None - indexvert = vertmode - indextime = itime - nlon = 1 - nlat = 1 - nvert = 1 - ntime = 1 - index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles) - print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime) - elif typefile in ['gcm']: + if typefile in ['meso','mesoapi'] and mapmode == 1: + indextime = itime + indexlon = None + indexlat = None + indexvert = nvert ## ou svert ??? + nlon = 1 + nlat = 1 + nvert = 1 + ntime = 1 + else: ## get all indexes to be taken into account for this subplot and then reduce field ## 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 - #print "(nplot-1)%nlon", (nplot-1)%nlon indexlon = getsindex(slon,(nplot-1)%nlon,lon) - print "indexlon", indexlon indexlat = getsindex(slat,((nplot-1)//nlon)%nlat,lat) indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) - index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles) - #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime) + index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles) + print nlon, nlat, nvert, ntime + print slon, slat, svert, stime + print index_f + print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime) #################################################################### @@ -208,5 +220,6 @@ #what_I_plot = what_I_plot*mult if not error: - if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) + if mapmode == 1: + if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) zevmin, zevmax = calculate_bounds(what_I_plot) zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20) @@ -214,5 +227,5 @@ if mapmode == 0: x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ - indextime,what_I_plot.shape, len(all_var[index_f].shape),vertmode) + indextime,what_I_plot.shape, len(all_var2[index_f].shape),vertmode) ### If we plot a 2-D field if len(what_I_plot.shape) is 2: @@ -221,5 +234,5 @@ ### If we plot a 1-D field elif len(what_I_plot.shape) is 1: - plot(x,y,what_I_plot, zelevels, colors='g', linewidths = 0.33 ) #colors='w' )# , alpha=0.5) + plot(what_I_plot,x) else: errormess("There is an error in reducing field !") @@ -239,43 +252,32 @@ # fvar = str(mult) + "*" + var ### - if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) + if mapmode == 1: + if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) + elif mapmode == 0: + what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ + indextime,what_I_plot, len(all_var[index_f].shape),vertmode) zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) if colorb in ["def","nobar"]: palette = get_cmap(name=defcolorb(fvar.upper())) else: palette = get_cmap(name=colorb) - ### If we plot a 2-D field + ##### 2D field if len(what_I_plot.shape) is 2: - if mapmode == 0: - what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ - indextime,what_I_plot, len(all_var[index_f].shape),vertmode) if not tile: if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax) #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20) - print "what_I_plot.shape", what_I_plot.shape - print "x.shape, y.shape", x.shape, y.shape - zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) + zelevels = np.linspace(zevmin,zevmax)#,num=ndiv+1) #contourf(what_I_plot, zelevels, cmap = palette ) - contourf(x,y,what_I_plot, zelevels, cmap = palette ) + contourf( x, y, what_I_plot, zelevels, cmap = palette ) else: if hole: what_I_plot = nolow(what_I_plot) - pcolor(x,y,what_I_plot, cmap = palette, \ - vmin=zevmin, vmax=zevmax ) - if colorb != 'nobar' and var != 'HGT' : - #subplot(111) - #colorbar() + pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax ) + if colorb != 'nobar' and var != 'HGT' : colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\ ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\ extend='neither',spacing='proportional') # both min max neither - ### If we plot a 1-D field + ##### 1D field elif len(what_I_plot.shape) is 1: - what_I_plot,x,y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ - indextime,what_I_plot,len(all_var[index_f].shape),vertmode) - # x = np.array(x) - # y = np.array(y) - print "what_I_plot.shape", what_I_plot.shape - print "x.shape", x.shape plot(x,what_I_plot) - ### If we plot something that is not a 2-D or 1-D field - ### (maybe plot 3-D field one day or movie ??) + #### Other cases: (maybe plot 3-D field one day or movie ??) else: print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!!"