Changeset 490

Timestamp:

Jan 4, 2012, 1:25:44 AM (14 years ago)

Author:

aslmd

Message:

UTIL PYTHON: several lat/lon now works for mesoscale in 1D mode. also added a capability to plot wind module if one of the wind component is asked as field. some stuff was added as commented for demonstration

Location:

trunk/UTIL/PYTHON

Files:

• myplot.py (modified) (2 diffs)
• planetoplot.py (modified) (8 diffs)

trunk/UTIL/PYTHON/myplot.py

@@ -751 +751 @@
              "QSURFICE":     "%.0f",\
              "UM":           "%.0f",\
+             "WIND":         "%.0f",\
              "ALBBARE":      "%.2f",\
              "TAU":          "%.1f",\
@@ -779 +780 @@
              "QH2O":         "PuBu",\
              "USTM":         "YlOrRd",\
+             "WIND":         "YlOrRd",\
              #"T_nadir_nit":  "RdBu_r",\
              #"T_nadir_day":  "RdBu_r",\

trunk/UTIL/PYTHON/planetoplot.py

@@ -76 +76 @@
     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
+    #from mpl_toolkits.basemap import cm
     import numpy as np
     from numpy.core.defchararray import find
     from videosink import VideoSink
     import subprocess
+    #from singlet import singlet
 
     ################################
@@ -94 +96 @@
     ################################
     nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
-    vlon = None ; vlat = None
-    if slon is not None: vlon = slon[0][0]
-    if slat is not None: vlat = slat[0][0]
     if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!")
     zelen = len(namefiles)*len(var)
@@ -157 +156 @@
          elif zarea is not None: [wlon,wlat] = latinterv(area=zarea)
 
-########################################################## LOAD 4D DIMENSIONS : x, y, z, t
+##########################################################
+############ LOAD 4D DIMENSIONS : x, y, z, t #############
+##########################################################
       if typefile == "gcm":
-          lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:]
+          lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:] ; vert = nc.variables["altitude"][:]
           if "Time" in nc.variables:      time = nc.variables["Time"][:]
           elif "time" in nc.variables:    time = nc.variables["time"][:]
           else:                           errormess("no time axis found.")
-          vert = nc.variables["altitude"][:]
           if axtime in ["ls","sol"]:   errormess("not supported. should not be too difficult though.")
       elif typefile in ['meso','mesoapi','geo','mesoideal']:
-          ### 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:   
-              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]] )
+          ###### STUFF TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
+          ###### principle: calculate correct indices then repopulate slon and slat
+          if slon is not None or slat is not None:
+              if firstfile and save == 'png' and typefile == 'meso':   iwantawhereplot = nc     #show a topo map with a cross on the chosen point
+              else:                                                    iwantawhereplot = None   #do not show anything, just select indices
+              numlon = 1 ; numlat = 1 
+              if slon is not None:   numlon = slon.shape[0]    
+              if slat is not None:   numlat = slat.shape[0]
+              indices = np.ones([numlon,numlat,2]) ; vlon = None ; vlat = None
+              for iii in range(numlon):  
+               for jjj in range(numlat):
+                 if slon is not None:  vlon = slon[iii][0]  ### note: slon[:][0] does not work
+                 if slat is not None:  vlat = slat[jjj][0]  ### note: slon[:][0] does not work
+                 indices[iii,jjj,:] = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)  
+                 lonp,latp = ( lon2d[indices[iii,jjj,0],indices[iii,jjj,1]] , lat2d[indices[iii,jjj,0],indices[iii,jjj,1]] )
+                 #print vlon, lonp, vlat, latp
+              for iii in range(numlon):
+               for jjj in range(numlat):
+                 if slon is not None: slon[iii][0] = indices[iii,0,1] ; slon[iii][1] = indices[iii,0,1]  #...this is idx
+                 if slat is not None: slat[jjj][0] = indices[0,jjj,0] ; slat[jjj][1] = indices[0,jjj,0]  #...this is idy
+              lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,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
-          if slon is not None: slon[0][0] = indices[1] ; slon[0][1] = indices[1]  #...this is idx
-          if slat is not None: slat[0][0] = indices[0] ; slat[0][1] = indices[0]  #...this is idy
+          ######
           if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
           else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
-          if (var2 is not None and var2 not in ['PHTOT','W']):
-               vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
-               dumped_vert_stag=True
-          else: dumped_vert_stag=False
+          if (var2 is not None and var2 not in ['PHTOT','W']): dumped_vert_stag=True ; vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
+          else:                                                dumped_vert_stag=False
           if varname in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
           else:                 latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
@@ -213 +226 @@
        #   errormess("Not the same latitudes !", [lat,lat0])
        ## Faire d'autre checks sur les compatibilites entre fichiers!!
+##########################################################
+##########################################################
 ##########################################################
 
@@ -297 +312 @@
        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()
+       ##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()
+       #truc = True 
+       #truc = False
+       #if truc: indexvert = None
        ####################################################################
        ########## REDUCE FIELDS
@@ -316 +334 @@
            vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
            vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
+           if varname in [uchar,vchar]: what_I_plot = np.sqrt( np.square(vecx) + np.square(vecy) ) ; varname = "wind"
+
+       #####################################################################
+       #if truc:
+       #   nx = what_I_plot.shape[2] ; ny = what_I_plot.shape[1] ; nz = what_I_plot.shape[0] 
+       #   for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / smooth(what_I_plot[k,:,:],12)
+       #   for iii in range(nx):
+       #    for jjj in range(ny):
+       #     deviation = what_I_plot[:,jjj,iii] ; mx = max(deviation) ; mn = min(deviation)
+       #     if iii > 6 and iii < nx-6 and jjj > 6 and jjj < ny-6:   what_I_plot[0,jjj,iii],rel = singlet(deviation,vert/1000.)  ### z must be in km
+       #     else:                                                   what_I_plot[0,jjj,iii]     = 0.
+       #     if np.abs(what_I_plot[0,jjj,iii]) > 1.5: 
+       #         print iii,jjj,what_I_plot[0,jjj,iii],int(abs(1.-mx)*100.),int(abs(1.-mn)*100.)
+       #         plot(rel)
+       #   show()
+       #   anomaly = True ### pour avoir les bons reglages plots
+       #   what_I_plot = what_I_plot[0,:,:]  
+       #####################################################################
+
        ####################################################################
        ### General plot settings
@@ -340 +377 @@
                if colorb in ["def","nobar"]:                           palette = get_cmap(name=defcolorb(fvar.upper()))
                else:                                                   palette = get_cmap(name=colorb)
+               #palette = cm.GMT_split
                ##### 1. ELIMINATE 0D or >3D CASES
                if len(what_I_plot.shape) == 0:   
@@ -391 +429 @@
                         if indextime is not None: lbl = lbl + " it" + str(indextime[0])
                         if mrate is not None: x = y  ## because swapaxes...
+                        #what_I_plot_frame = np.diff(what_I_plot_frame, n=1) ; x = x[1:]
                         if indexvert is not None or indextime is None:    plot(x,what_I_plot_frame,label=lbl)  ## regular plot
                         else:                                             plot(what_I_plot_frame,x,label=lbl)  ## vertical profile