Changeset 640

Timestamp:

Apr 30, 2012, 1:09:03 PM (14 years ago)

Author:

tnavarro

Message:

Ls axtime for gcm files. Plot files without time axis (surface,etc ...). Histogram plot for difference between two files (--ope -)

File:

• trunk/UTIL/PYTHON/planetoplot.py (modified) (4 diffs)

trunk/UTIL/PYTHON/planetoplot.py

@@ -81 +81 @@
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
-    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, subplots_adjust, axes, clabel
+    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
     from matplotlib.cm import get_cmap
     #from mpl_toolkits.basemap import cm
@@ -87 +87 @@
     from numpy.core.defchararray import find
     from videosink import VideoSink
+    from timestuff import sol2ls
     import subprocess
     #from singlet import singlet
@@ -197 +198 @@
           elif "time" in nc.variables:          time = nc.variables["time"][:]
           elif "time_counter" in nc.variables:  time = nc.variables["time_counter"][:]/86400. #### time counter cinverstion from s-> days
-          else:                                 errormess("no time axis found.")
-          if axtime in ["ls","sol"]:   errormess("not supported. should not be too difficult though.")
+          else:                                 time = [0.] #errormess("no time axis found.")
+          if axtime in ["ls"]:
+              print "converting to Ls ..."
+              for iii in range(len(time)):
+                time[iii] = sol2ls(time[iii]) # to use with timestuff
+                if iii > 0:
+                  while abs(time[iii]-time[iii-1]) > 300:
+                    time[iii] = time[iii]+360
           # for 1D plots (no need for longitude computation):
           if axtime in ["lt"]:
@@ -572 +579 @@
                                              scale=20., factor=250., color=colorvec, key=key)
                                                               #200.         ## or csmooth=stride
+                        if ope == '-' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
+                            subplot(subv,subh,nplot+1)
+                            latmin = -50.; latmax = 50. # latitude range for histogram of difference
+                            if indexlat is None: # a bit dirty, if field is not reduced along latitude, we assume lat is along the y axis
+                                zeindexlat = (lat<latmax)*(lat>latmin)
+                                what_I_plot_frame = what_I_plot_frame[zeindexlat,:]
+                            toplot = np.ravel(what_I_plot_frame[np.isnan(what_I_plot_frame)==False])
+                            zebins = np.linspace(minop,maxop,num=30)
+                            hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)
+                            zestd = np.std(toplot)
+                            zemean = np.mean(toplot)
+                            zebins = np.linspace(minop,maxop,num=300)
+                            zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)**2 / (2 * zestd**2) ) )
+                            text(zebins[2],np.max(zegauss),'mean: '+str(zemean)[0:5]+'\nstd: '+str(zestd)[0:5],fontsize = 12)
+                            plot(zebins, zegauss, linewidth=1, color='r')
+                            title("Histogram fig(1) "+ope+" fig(2)")
+                            subplot(subv,subh,nplot) # go back to last plot for title of contour difference
+
                     elif which == "contour":
                         zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))