Changeset 724

Timestamp:

Jul 16, 2012, 11:15:30 PM (12 years ago)

Author:

aslmd

Message:

UTIL PYTHON. erase a duplicate in planetoplot about defining lat and lon. added the capability to guess dimension arrays if not present. this allows to plot virtually anything in any NETCDF file, in particular file with 4D or 3D variables included but no clear reference to lat lon vert time, just for instance x y z t. convenient for extracted fields.

Location:

trunk/UTIL/PYTHON

Files:

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

trunk/UTIL/PYTHON/myplot.py

@@ -645 +645 @@
             print "WARNING: domain plot artificially centered on lat,lon 0,0"
     elif typefile in ['gcm','earthgcm','ecmwf']: 
+        #### n est ce pas nc.variables ?
         if "longitude" in nc.dimensions:  dalon = "longitude"
         elif "lon" in nc.dimensions:      dalon = "lon"
+        else:                             dalon = "nothing"
         if "latitude" in nc.dimensions:   dalat = "latitude"
         elif "lat" in nc.dimensions:      dalat = "lat"
+        else:                             dalat = "nothing"
         [lon2d,lat2d] = getcoord2d(nc,nlat=dalat,nlon=dalon,is1d=True)
     elif typefile in ['geo']:
@@ -657 +660 @@
 def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
     import numpy as np
-    if is1d:
-        lat = nc.variables[nlat][:]
-        lon = nc.variables[nlon][:]
+    if nlon == "nothing" or nlat == "nothing":
+        print "NO LAT LON FIELDS. I AM TRYING MY BEST. I ASSUME GLOBAL FIELD."
+        lon = np.linspace(-180.,180.,getdimfromvar(nc)[-1])
+        lat = np.linspace(-90.,90.,getdimfromvar(nc)[-2])
         [lon2d,lat2d] = np.meshgrid(lon,lat)
     else:
-        lat = nc.variables[nlat][0,:,:]
-        lon = nc.variables[nlon][0,:,:]
-        [lon2d,lat2d] = [lon,lat]
+        if is1d:
+            lat = nc.variables[nlat][:]
+            lon = nc.variables[nlon][:]
+            [lon2d,lat2d] = np.meshgrid(lon,lat)
+        else:
+            lat = nc.variables[nlat][0,:,:]
+            lon = nc.variables[nlon][0,:,:]
+            [lon2d,lat2d] = [lon,lat]
     return lon2d,lat2d
+
+## Author: AS
+def getdimfromvar (nc):
+    varinfile = nc.variables.keys()
+    dim = nc.variables[varinfile[-1]].shape ## usually the last variable is 4D or 3D
+    return dim
 
 ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth

trunk/UTIL/PYTHON/planetoplot.py

@@ -80 +80 @@
                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
                        getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
-                       windamplitude,slopeamplitude,deltat0t1,enrichment_factor
+                       windamplitude,slopeamplitude,deltat0t1,enrichment_factor,getdimfromvar
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
@@ -184 +184 @@
           if slon is not None: sslon = slon  
           if slat is not None: sslat = slat  
-          ### LAT
-          if "latitude" in nc.variables:  lat = nc.variables["latitude"][:]
-          elif "lat" in nc.variables:     lat = nc.variables["lat"][:]
-          ### LON
-          if "longitude" in nc.variables: lon = nc.variables["longitude"][:]
-          elif "lon" in nc.variables:     lon = nc.variables["lon"][:]
+          ### LAT : done in getcoorddef
+          lat = lat2d[0,:] 
+          ### LON : done in getcoorddef
+          lon = lon2d[:,0]
           ### ALT
           if "altitude" in nc.variables:   vert = nc.variables["altitude"][:]
@@ -195 +193 @@
           elif "lev" in nc.variables:      vert = nc.variables["lev"][:]
           elif "presnivs" in nc.variables: vert = nc.variables["presnivs"][:]
-          else:                            vert = [0.]
+          else: 
+              print "No altitude found. Try to build a simple axis." 
+              dadim = getdimfromvar(nc) 
+              if   len(dadim) == 4:  print "-- 4D field. Assume z is dim 2." ; vert = np.arange(dadim[-3])
+              elif len(dadim) == 3:  print "-- 3D field. Assume z is dim 1." ; vert = [0.]
+              else:                  errormess("-- 2D or 1D field. Unclear. Stop.")
+          #else:                            vert = [0.]
           ### AIRE (to weight means with the area)
           if "aire" in nc.variables:      area = nc.variables["aire"][:,:]  
@@ -203 +207 @@
           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:                                 time = [0.] #errormess("no time axis found.")
+          else: 
+              print "No time found. Try to build a simple axis. Assume t is dim 1." 
+              dadim = getdimfromvar(nc)
+              if   len(dadim) == 4:  time = np.arange(dadim[-4])
+              elif len(dadim) == 3:  time = np.arange(dadim[-3])
+              else:                  errormess("-- I am sorry this is 2D or 1D field. I failed.")
+          #else:                                time = [0.] #errormess("no time axis found.")
           if axtime in ["ls"]:
               print "converting to Ls ..."