Changeset 637

Timestamp:

Apr 30, 2012, 11:09:59 AM (13 years ago)

Author:

aslmd

Message:

Working whether dimension is latitude/longitude or lat/lon (so, OK for UK assimilation files now). Better centering of lcc projection. Labels for contour plots.

Location:

trunk/UTIL/PYTHON

Files:

• myplot.py (modified) (9 diffs)
• myscript.py (modified) (2 diffs)
• planetoplot.py (modified) (5 diffs)

trunk/UTIL/PYTHON/myplot.py

@@ -290 +290 @@
     from mymath import max,mean
     from scipy import integrate
+    import numpy as np
     if yint and vert is not None and indice is not None:
        if type(input).__name__=='MaskedArray':
@@ -320 +321 @@
         fig.subplots_adjust(hspace = 0.5)
         rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
-    elif   numplot == 4:
+    elif numplot == 4:
         subv = 2
         subh = 2
@@ -534 +535 @@
             [lon2d,lat2d] = np.meshgrid(falselon,falselat) ## dummy coordinates
             print "WARNING: domain plot artificially centered on lat,lon 0,0"
-    elif typefile in ['gcm']: 
-        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
-    elif typefile in ['earthgcm','ecmwf']: 
-        [lon2d,lat2d] = getcoord2d(nc,nlat="lat",nlon="lon",is1d=True)
+    elif typefile in ['gcm','earthgcm','ecmwf']: 
+        if "longitude" in nc.dimensions:  dalon = "longitude"
+        elif "lon" in nc.dimensions:      dalon = "lon"
+        if "latitude" in nc.dimensions:   dalat = "latitude"
+        elif "lat" in nc.dimensions:      dalat = "lat"
+        [lon2d,lat2d] = getcoord2d(nc,nlat=dalat,nlon=dalon,is1d=True)
     elif typefile in ['geo']:
         [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
@@ -692 +695 @@
     meanlon = 0.5*(wlon[0]+wlon[1])
     meanlat = 0.5*(wlat[0]+wlat[1])
+    zewidth = np.abs(wlon[0]-wlon[1])*60000.*np.cos(3.14*meanlat/180.)
+    zeheight = np.abs(wlat[0]-wlat[1])*60000.
     if blat is None:
         ortholat=meanlat
@@ -703 +708 @@
     h = 50.  ## en km
     radius = 3397200.
+    #print meanlat, meanlon
     if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
                               llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
@@ -708 +714 @@
     elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=ortholon,lat_0=ortholat)
     elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
-                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
+                              width=zewidth,height=zeheight)
+                              #llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
     elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
     elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=180.)
     elif char == "nplaea":  m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
     elif char == "laea":    m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
-                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
+                              width=zewidth,height=zeheight)
+                              #llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
     elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
     elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
@@ -724 +732 @@
     zelinewidth = 1
     zelatmax = 80
+    if meanlat > 75.: zelatmax = 90. ; step = step/2.
     # to show gcm grid:
     #zecolor = 'r'
@@ -833 +842 @@
              "DOWNDRAFT":    "%.0f",\
              "TK":           "%.0f",\
+             "T":            "%.0f",\
              #"ZMAX_TH":      "%.0f",\
              #"WSTAR":        "%.0f",\
@@ -1014 +1024 @@
         "South_Pole":            [[-90.,-50.],[-180., 180.]],\
         "Close_South_Pole":      [[-90.,-75.],[-180., 180.]],\
-        "Sirenum_Crater_large":  [[-46.,-34.],[-166., -151.]],\
-        "Sirenum_Crater_small":  [[-36.,-26.],[-168., -156.]],\
-
+        "Sirenum_Crater_large":  [[-46.,-34.],[-166.,-151.]],\
+        "Sirenum_Crater_small":  [[-36.,-26.],[-168.,-156.]],\
+        "Rupes":                 [[ 72., 90.],[-120.,- 20.]],\
               }
     if area not in list:   area = "Whole"

trunk/UTIL/PYTHON/myscript.py

@@ -39 +39 @@
     parser.add_option('--ylabel',       action='store',dest='ylab',       type="string",  default=None, help='customize the y-axis label')
     parser.add_option('--labels',       action='store',dest='labels',     type="string",  default=None, help='customize 1D curve labels. Str comma-separated. [None]')
-    parser.add_option('--lstyle',       action='store',dest='linestyle',  type="string",  default=None, help='customize 1D curve linestyles. String separated by commas. [None]')
+    parser.add_option('--lstyle',       action='store',dest='linestyle',  type="string",  default=None, help='customize 1D curve linestyles. Str comma-separ. [None]')
 
     ### SPECIFIC FOR MAPPING [MAPMODE 1]
@@ -49 +49 @@
     parser.add_option('--blat',         action='store',dest='blat',      type="int",     default=None,  help='reference lat (or bounding lat for stere) [computed]')
     parser.add_option('--blon',         action='store',dest='blon',      type="int",     default=None,  help='reference lon [computed]')
-    parser.add_option('--mark',         action='append',dest='mark',  type="string",  default=None, help='superimpose a crossmark at given coords. lon,lat float separated by commas. [None]')
+    parser.add_option('--mark',         action='append',dest='mark',  type="string",  default=None, help='superimpose a crossmark at given lon,lat [None]')
 
     ### SPECIFIC FOR SLICING [MAPMODE 0]

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
+    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.cm import get_cmap
     #from mpl_toolkits.basemap import cm
@@ -141 +141 @@
       elif typefile != typefile0:   errormess("Not the same kind of files !", [typefile0, typefile])
     ### we care for input file being 1D simulations 
-      if typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1:       mapmode=0 ; winds=False
-      elif typefile in ['earthgcm'] and len(nc.variables["lon"][:]) is 1 and len(nc.variables["lat"][:]) is 1:           mapmode=0 ; winds=False
+      if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc.variables["longitude"][:])*len(nc.variables["latitude"][:])
+      elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc.variables["lon"][:])*len(nc.variables["lat"][:])
+      if typefile in ['gcm','earthgcm'] and is1d == 1:       mapmode=0 ; winds=False
     ### we create default vert and time prescriptions if not here in case mapping mode is on (lat/lon maps)
       if redope is None and mapmode == 1:
@@ -282 +283 @@
       all_time[k] = time
       if var2: 
+         all_var2[k] = getfield(nc,var2)
          if ((var2 in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (var2 not in nc.variables)): all_var2[k] = slopeamplitude(nc)
-         else: all_var2[k] = getfield(nc,var2)
       if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
       #####################
@@ -561 +562 @@
                                 [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
                                 key = True
+                                if fvar in ['UV','uv','uvmet']: key = False
                             elif typefile in ['gcm']:
                                 key = False
@@ -573 +575 @@
                         zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
                         zelevels = np.linspace(zevminc,zevmaxc,ticks/2) #20)
-                        if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,1000.)
+                        if var2 == 'HGT':     zelevels = np.arange(-10000.,30000.,1000.)
+                        elif var2 == 'tpot':  zelevels = np.arange(270,370,5)
+                        elif var2 == 'tk':    zelevels = np.arange(150,250,5)
                         if mapmode == 0:   
-                            what_I_plot_frame, x, y = define_axis( lon,lat,vert,time,indexlon,indexlat,indexvert,\
+                            what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
                                                               itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode )
-                            cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33, alpha=0.5, linestyles='solid')
-                        elif mapmode == 1:  cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33, alpha=0.5, linestyles='solid')
+                            ## is this needed? only if len(all_var2[index_f].shape) != len(all_var[index_f].shape)
+                            cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
+                            clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper()))
+                        elif mapmode == 1:  cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
                     if which in ["regular","unidim"]: