Changeset 489 for trunk/UTIL/PYTHON/planetoplot.py

Timestamp:

Dec 21, 2011, 3:22:19 PM (13 years ago)

Author:

aslmd

Message:

UTIL PYTHON : 1) no more negative values for --time but use --axtime lt ; this actually fix a stupid bug with negative lat lon. 2) use of -l is now possible with GCM files, try it with -i 2. 3) now correct handling of sol/ls in mesoscale files 4) for mesoscale files with chosen lat+lon now shows a topo map with a cross on it.

File:

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

trunk/UTIL/PYTHON/planetoplot.py

@@ -74 +74 @@
     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
+    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel, axis
     from matplotlib.cm import get_cmap
     import numpy as np
@@ -147 +147 @@
       if ((proj == None) and (typefile not in ['mesoideal'])):   proj = getproj(nc)                 
 
-##########################################################
+      if firstfile:
+         ##########################
+         ### Define plot boundaries
+         ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
+         if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
+         elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
+         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
+         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom)
+         elif zarea is not None: [wlon,wlat] = latinterv(area=zarea)
+
+########################################################## LOAD 4D DIMENSIONS : x, y, z, t
       if typefile == "gcm":
           lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:]
@@ -158 +168 @@
           ### the following lines are kind of dirty... not possible to ask for several lats and lons
           if vlon is not None or vlat is not None:   
-              indices = bidimfind(lon2d,lat2d,vlon,vlat) #,file=nc)  #placer un point sur carte
+              if firstfile:  iwantawhereplot = nc
+              else:          iwantawhereplot = None
+              indices = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)  #placer un point sur carte
               lonp,latp = ( lon2d[indices[0],indices[1]] , lat2d[indices[0],indices[1]] )
           if typefile in ['mesoapi','meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  ### important to do that now and not before
@@ -183 +195 @@
               elif "Time" in nc.variables:  time = count + np.arange(0,len(nc.variables["Time"]),1)
               else:                         time = count + np.arange(0,1,1)
-              count = time[-1] + 1  ## so that a cat is possible with simple subscripts
+              if nnn > 0:  count = time[-1] + 1  ## so that a cat is possible with simple subscripts
+              else:        count = 0
+          if axtime in ["lt"]:
+              for i in range(len(time)):  time[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
+              print "LOCAL TIMES.... ", time
           ###
           if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
@@ -198 +214 @@
        ## Faire d'autre checks sur les compatibilites entre fichiers!!
 ##########################################################
-
-      if firstfile:
-         ##########################
-         ### Define plot boundaries
-         ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
-         if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
-         elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
-         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
-         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom)
-         elif zarea is not None: [wlon,wlat] = latinterv(area=zarea) 
 
       all_varname[k] = varname
@@ -246 +252 @@
                     if var2: tab2 = np.append(tab2,all_var2[k],axis=0) 
                     tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
-                #all_time[0] = np.arange(0,len(tab),1) ### AS: time reference is too simplistic, should be better
                 all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
                 if var2: all_var2[0] = np.array(tab2)
@@ -287 +292 @@
        else:                            yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
        time = all_time[index_f]
-       if stime is not None:
-           if stime[0][0] < 0:
-               if typefile in ['mesoapi','meso']:
-                   for i in range(len(time)):  time[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
-                   print "OK... WORKING WITH LOCAL TIMES"
-               else: errormess("local times not supported. not too hard to modify the code though.")
        if mrate is not None:                 indextime = None 
        else:                                 indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
@@ -298 +297 @@
        if typefile in ['mesoapi','meso'] and indextime is not None:  ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) )  
        print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
+       #var = nc.variables["phisinit"][:,:]
+       #contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)
+       #show()
+       #exit()
        ####################################################################
        ########## REDUCE FIELDS
@@ -337 +340 @@
                if colorb in ["def","nobar"]:                           palette = get_cmap(name=defcolorb(fvar.upper()))
                else:                                                   palette = get_cmap(name=colorb)
-               ##### 1. ELIMINATE >3D CASES
-               if len(what_I_plot.shape) >= 4:
+               ##### 1. ELIMINATE 0D or >3D CASES
+               if len(what_I_plot.shape) == 0:   
+                 print "VALUE VALUE VALUE VALUE ::: ",what_I_plot
+                 save = 'donothing'
+               elif len(what_I_plot.shape) >= 4:
                  print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
                  errormess("Are you sure you did not forget to prescribe a dimension ?")
@@ -542 +548 @@
         elif save in ['eps','svg','pdf']:     makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)
         elif save == 'gui':                   show()
+        elif save == 'donothing':             pass
         elif save == 'txt':                   print "Saved results in txt file." 
         else: