#!/usr/bin/python -tt #======================================================================= # General Documentation """Define/calculate atmospheric constants and quantities. The atmqty package has the following classes and functions: (1) The AtmQty class of all "atmospheric quantities" (defined below). All such quantities in a given instantiation of the class share the same latitude, longitude, level locations. (2) The stand-alone AtmConst class of atmosphere-related con- stants. (3) Stand-alone functions to calculate atmospheric quantities that are directly derivative from standard state variables. Most (though not all) of these functions are referenced by methods in AtmQty, but are called "stand-alone" because they can also be used on objects that are not AtmQty attributes. Used in a stand-alone fashion, these functions can operate on arrays of any size and dimension and are not confined to a spherical or vertical grid. (4) Stand-alone functions to manipulate objects of the AtmQty class (e.g. a regridding function that moves all variables from one set of latitude, longitude, level locations to another). "Atmospheric quantities" are defined as quantities that are stan- dard state variables or are directly derivative (i.e. not requi- ring time integration or modeling) from standard state variables. The AtmQty class is defined in __init__.py. The AtmConst class and each of the stand-alone functions are defined individually in separate modules in this package; AtmConst is defined in module atmconst. For the stand-alone functions, the function name is the same as the module name (e.g. function "theta" is defined in module "theta"). The stand-alone functions may not be "truly" stand-alone since they may require methods in other modules in this package; this package is designed to be used as a unit. Some modules, like the atmconst module, however, are truly stand-alone as they requires no other modules besides themselves; see the notes in the beginning of the module source code for details. The package is written such that a single command "import atmqty" to import the package will also import all the submodules in the package. Thus, one does not need to type in individual "import atmqty.M" commands for each M submodule but rather can use the name atmqty.M as a prefix to the included functions immediately after an "import atmqty" command. Some useful online help commands for the package: * help(atmqty): Help for the package (including the AtmQty class). A list of all modules in this package is found in the "Package Contents" section of the help output. * help(atmqty.M): Details of each module "M", where "M" is the module's name. Example 1: Atmospheric constants: >>> from atmconst import AtmConst >>> const = AtmConst() >>> '%.6f' % const.R_d '287.050000' Example 2: Stand-alone function: >>> from theta import theta >>> p = N.array([1000., 850., 500.]) >>> T = N.array([273.1, 257.3, 233.1]) >>> var = theta(p,T) >>> ['%.6f' % var[i] for i in range(3)] ['273.100000', '269.537605', '284.189919'] Example 3: Single-level x-y domain AtmQty object: (a) Initialize the AtmQty object: >>> import atmqty >>> import Numeric as N >>> lat = N.arange(13) * 15.0 - 90.0 >>> lon = N.arange(25) * 15.0 - 180.0 >>> lev = N.array(850.) >>> x = atmqty.AtmQty(lat, lon, lev, missing=1e29) (b) The pattern "base" is two "distorted bulls-eyes" of opposite polarity, centered over the poles (and add extra shallow dimension to represent the single level): >>> base = N.outerproduct(N.sin(lat*N.pi/360.), N.cos(lon*N.pi/360.)) >>> base = N.reshape(base, (base.shape[0], base.shape[1], 1)) (c) Create temperature and moisture "data" using the "base" pattern and add to AtmQty object (a second, redundant call to add_qty is provided to show the method accepts both types of argument syntax shown): >>> T = (base * 10.) + 257.3 >>> q = N.absolute(base) * 5e-3 >>> x.add_qty( {'q':q, 'T':T} ) >>> x.add_qty( {'q':q}, {'T':T} ) (d) Calculate potential temperature and output temperature and potential temperature: >>> x.theta() >>> ['%.6f' % x.get_qty('T')[i,2,0] for i in range(3)] ['255.469873', '255.724409', '256.005905'] >>> ['%.6f' % x.get_qty('theta')[i,2,0] for i in range(3)] ['267.620434', '267.887077', '268.181961'] (e) Set T[1,2,0] to missing and calculate eice (note that the other quantities, such as theta, that are dependent on T do not recalculate to take into account the new missing value): >>> x.get_qty('T')[1,2,0] = 1e29 >>> x.eice() >>> ['%.7g' % x.get_qty('eice')[i,2,0] for i in range(3)] ['1.285009', '1e+29', '1.351447'] >>> ['%.6f' % x.get_qty('theta')[i,2,0] for i in range(3)] ['267.620434', '267.887077', '268.181961'] """ #----------------------------------------------------------------------- # Additional Documentation # # Package Name: # atmqty # # RCS Revision Code: # $Id: __init__.py,v 1.44 2004/05/27 20:55:21 jlin Exp $ # # Modification History: # - 19 Nov 2003: Original by Johnny Lin, Computation Institute, # University of Chicago. Passed passably reasonable tests. Class # can handle missing values. # - 16 Mar 2004: Added geopotential altitude, static stability. # Passed passable reasonable tests. # - 29 Mar 2004: Added isnew, has_qty, list_all_qty. Passed # minimally passably reasonable tests. # - 20 Apr 2004: Fix bug with AtmQty method alt. Passed minimally # passably reasonable tests. # - 26 May 2004: Add methods get_qty_asMA, ipv, vort_rel, attribute # fc, and other methods/attributes. Passed minimally passably # reasonable tests. # # Notes: # - Written for Python 2.1.3. # - Code for data array in Example 3 came from Dean Williams. # - Non-"built-in" packages and modules required: Numeric, MA. See # import statements throughout for more information. # # Copyright (c) 2003-2004 by Johnny Lin. For licensing, distribution # conditions, contact information, and additional documentation see # the URL http://www.johnny-lin.com/py_pkgs/atmqty/doc/. #======================================================================= #---------------- Module General Import and Declarations --------------- #- General imports: import copy import MA import Numeric as N #- List of modules in package and import all these modules so they # are all available to the user with a single "import atmqty" # command: __all__ = [ "atmconst" \ , "eice" \ , "esat" \ , "is_numeric_float" \ , "mix_ratio" \ , "package_version" \ , "press" \ , "press2alt" \ , "rho" \ , "transform_AQobj" \ , "static_stability" \ , "theta" \ , "vapor_press" \ , "virt_temp" ] import atmconst import eice import esat import is_numeric_float import mix_ratio import press import press2alt import rho import static_stability import theta import transform_AQobj import vapor_press import virt_temp #- Set package version number: import package_version __version__ = package_version.version __author__ = package_version.author __date__ = package_version.date __credits__ = package_version.credits #------------------------ Class AtmQty: Header ------------------------ class AtmQty: """Atmospheric quantities in a latitude, longitude, level domain. All attributes in this class that are "atmospheric quantities" are assumed to be on the same latitude, longitude, level slab on a spherical grid. Latitude and longitude are always given in decimal degrees. The level coordinate can be different units depending on the instantiation of the object. Methods operating on this object are consistent with each other (i.e. equations for all quantities are physically consistent) and thus all quantities for a given instantiation of the object, for a given set of initial state variables (added to the object by the class method add_qty) are consistent with each other (assuming there is no recalculation of already calculated values during that instantiation of the object). Usually latitude and longitude span the entire sphere, but this is not required. Vertically the domain can encompass any number of levels, with a minimum of one level. (Of course, certain methods, such as centered differencing, that require a minimum number of points, will not work if the number of points in the domain is not large enough. The method will fail appropriately in such cases.) However, the domain must be contiguous in all three spatial dimensions. This class allows one to define a missing value (self.missing). All quantities that have this value are considered missing and calculations done using those values will return a missing value. Note that the same location may have missing value in some quantities and not others; missing values do not have to be co-located. This type of tracking of missing values is used instead of masks to save memory; however, some use of the MA module is made to help us temporarily manage missing values. Note that since missing values are stored as self.missing, all long-term storage of quantities and argument passing (in and out) of quantities is through plain Numeric arrays, not MA arrays. Object Attributes (besides self.qty): * self.const: Atmospheric constants. An instance of class AtmConst in module atmconst. * self.lat: Latitude [deg] vector for all "quantities". Mono- tonically increasing. * self.lon: Longitude [deg] vector for all "quantities". Mono- tonically increasing. * self.lev: Vertical level for all quantities. Level "type" is given in self.lev_type. The units of lev depends on lev_type (see description of self.lev_type for details). Monotonic, either increasing or decreasing. * self.lev_type: The level "type" given in dimension lev, and which can have the following values: + "pressure": Vertical levels are in pressure coordinates and lev is in units hPa (default value of lev_type). + "isentropic": Vertical levels are in isentropic (constant potential temperature levels) coordinates and lev is in units K. * self.missing: Missing value value. Floating point scalar. If undefined at object instantiation, is set to 1e+20. * self.P_surf: Surface pressure [hPa]. Used to calculate altitude. Either floating point scalar or a Numeric array of shape (Numeric.size(self.lat), Numeric.size(self.lon)). Has no missing values. Value is set to None on object instantiation. The attribute value can be changed directly or via the P_surf keyword in class method alt. * self.p0: Reference pressure used to calculate potential temperature [hPa]. Floating point scalar. If lev_type is "isentropic," this is the reference pressure associated with the isentropic levels and is set on object instantiation. If lev_type is "pressure", this attribute is the reference pres- sure used in calculating potential temperature; it is thus undefined until potential temperature is calculated. In both cases, this attribute should not be changed directly. * self.qty_shape: The shape of all self.qty quantities. 3- integer tuple (Numeric.size(self.lat), Numeric.size(self.lon), Numeric.size(self.lev)). Defined on object instantiation. * self.T_surf: Surface temperature [K]. Used to calculate altitude. Either floating point scalar or a Numeric array of shape (Numeric.size(self.lat), Numeric.size(self.lon)). Has no missing values. Value is set to None on object instantiation. The attribute value can be changed directly or via the T_surf keyword in class method alt. Quantities Data Structure Object Attribute: * self.qty: Dictionary of quantities variables. All quan- tities (but not boundary conditions), whether they are input parameters or are calculated from input parameters, are stored in this dictionary. Possible key values include: + 'alt': Geopotential altitude, based on 1976 U.S. Standard Atmosphere [m] + 'e': Vapor pressure (with respect to over water) [hPa]. + 'eice': Saturation vapor pressure over ice [hPa]. + 'esat': Saturation vapor pressure over water [hPa]. + 'fc': Coriolis parameter [1/s]. + 'gpot': Geopotential height [m]. + 'N2': Static stability (a.k.a. square of the buoyancy fre- quency) [1/s**2]. + 'p': Pressure [hPa]. + 'pv': Potential vorticity [m^2 s^-1 K kg^-1]. + 'q': Specific humidity [kg/kg] + 'r': Mixing ratio [kg/kg]. + 'RH': Relative humidity (based on esat) [unitless fraction]. + 'rho': Air density (including moisture effects). + 'T': Temperature [K]. + 'theta': Potential temperature [K]. + 'theta_e': Equivalent potential temperature [K]. + 'T_v': Virtual temperature [K]. + 'u': Zonal velocity [m/s]. + 'v': Meridional velocity [m/s]. + 'vort_rel': Relative vorticity [1/s]. Values corresponding to these keys are Numeric floating point arrays, with a shape specified by self.qty_shape. For details regarding attributes set on object instantiation, particularly regarding which ones must be fixed for the life of the instance, see the manual section on the AtmQty class (sub- topic "Object Instantiation") at: http://www.johnny-lin.com/py_pkgs/atmqty/doc/manual.html. The class method add_qty adds variable(s) to the self.qty dictionary. If a quantity is an input parameter state variable, you have to use this method first to add the variable to self.qty in order for it to be available for calculations. This is used instead of requiring variables be in argument lists in class methods, to simplify bookkeeping of what variables are defined and what are not. Quantities are deleted from self.qty by the class del_qty method. One can obtain a quantity from the self.qty data structure with the get_qty method. Note that depending on what lev_type is the quantities pressure and potential temperature may be an input state variable or a calculated derived variables. Also, because all state and derived quantities are in self.qty, when an object method is called to calculate a derived quantity, the only inputs that are supplied in the method argument list are optional ones. References Quoted in this Class Definition: * Holton, J. R. (1992): An Introduction to Dynamic Meteorology, 3rd edition. San Diego, CA: Academic Press, 511 pp. """ #------------- AtmQty Function __init__: Initialize Class ------------- def __init__( self, lat, lon, lev \ , lev_type="pressure" \ , missing=1e+20 \ , p0=None ): """Initialize basic class attributes. Method Arguments: * lat: Latitude [deg]. Numeric floating point vector. * lon: Longitude [deg]. Numeric floating point vector. * lev: Vertical level. Numeric floating point vector. * lev_type: The level "type" given in dimension lev. Scalar string. * missing: Missing value value. Numeric floating point scalar, if defined. Default is 1e+20. * p0: Reference pressure [hPa]. Floating point scalar. Details regarding arguments are in the docstring for the class. Note initializing an object of this class requires the gemath package. """ #- Check domain inputs are floating, monotonic, and increasing # for lat and lon: from is_numeric_float import is_numeric_float from gemath.is_monotonic_decr import is_monotonic_decr from gemath.is_monotonic_incr import is_monotonic_incr if is_numeric_float(lat, lon, lev) != 1: raise TypeError, "AtmQty: Args not Numeric floating" if (N.size(lat) > 1) and (N.size(lon) > 1) and (N.size(lev) > 1): if (not is_monotonic_incr(lat)) or (not is_monotonic_incr(lon)): raise ValueError, "AtmQty: Input(s) not monotonic increasing" if (not is_monotonic_incr(lev)) and (not is_monotonic_decr(lev)): raise ValueError, "AtmQty: Levels not monotonic" #- Initialize most domain describing attributes: from atmconst import AtmConst self.const = AtmConst() self.lat = lat self.lon = lon self.lev = lev self.lev_type = lev_type self.missing = missing #- Initialize self.qty to None to ensure the first call of # add_qty will create the dictionary. Initialize P_surf and # T_surf to None. Initialize qty_shape: self.qty = None self.P_surf = None self.T_surf = None self.qty_shape = \ ( N.size(self.lat), N.size(self.lon), N.size(self.lev) ) #- Initialize self.qty['p'] or self.qty['theta'] (and # self.p0 if applicable): tmp = N.zeros( self.qty_shape, typecode=N.Float ) for ilev in range(len(self.lev)): tmp[:,:,ilev] = self.lev[ilev] if self.lev_type == 'pressure': self.add_qty({'p':tmp}) elif self.lev_type == 'isentropic': self.add_qty({'theta':tmp}) if p0 != None: self.p0 = p0 else: self.p0 = 1000. else: raise ValueError, "AtmQty: Bad lev_type" del tmp #----- AtmQty Method: Add Variable(s) To Quantities Data Structure ---- def add_qty(self, *quantities): """Add/replace variable(s) to quantities data structure. Method Arguments: * *quantities: Each argument is a dictionary with at least one key:value pair where the key is the name of the variable and the value is the variable. The key is a string scalar and the value is a Numeric floating point array. See the class docstring for details on possible values of the key and additional conditions the value must meet. There can be as many arguments as desired in the add_qty method call. The method adds all the key:value pairs in all the arguments to the self.qty quantities data structure of the object. If self.qty already contains a term, this method replaces that term's entry with the input argument value. Method checks that each added quantity has the same shape as attribute self.qty_shape. Examples of syntax: x = atmqty.AtmQty(lat, lon, lev, missing=1e29) x.add_qty( {'e':e} ) x.add_qty( {'q':q}, {'T':T} ) x.add_qty( {'RH':RH, 'esat':esat, 'r':mix_ratio} ) """ for quantity in quantities: if type(quantity) != type({}): raise ValueError, "AtmQty: Bad input type" for aqty in quantity.keys(): if self.qty == None: self.qty = quantity else: self.qty[aqty] = quantity[aqty] if self.qty[aqty].shape != self.qty_shape: raise ValueError, "AtmQty: Bad add_qty shape" #--- AtmQty Method: Delete Variable(s) To Quantities Data Structure --- def del_qty(self, *quantities): """Delete variable(s) from quantities data structure. Method Arguments: * *quantities: Each argument is a string scalar of the key name of the variable in the quantities data structure. See the class docstring for details on values of the variable name. Note if this method is used to remove all entries in self.qty, self.qty is not set to None but becomes an empty dictionary. """ for quantity in quantities: del self.qty[quantity] #--- AtmQty Method: Obtain Variable From Quantities Data Structure ---- def get_qty(self, quantity): """Get variable from quantities data structure. Method Argument: * quantity: String scalar of the key name of the variable in the quantities data structure to obtain. See the class docstring for details on values of the variable name. Note this method does not remove the variable from the quantities data structure. """ return self.qty[quantity] #----- AtmQty Method: Obtain Quantities Variable as a Masked Array ---- def get_qty_asMA(self, quantity): """Get variable from quantities data structure in MA form. Method Argument: * quantity: String scalar of the key name of the variable in the quantities data structure to obtain. The returned MA array is formed using the copy=0 flag in the MA.masked_values function, and thus follows the rules of that package regarding whether the return value is a reference or a memory copy. See the class docstring for details on values of the variable name. Note this method does not remove the variable from the quantities data structure, nor change the entry in self.qty to MA form. """ return MA.masked_values(self.get_qty(quantity), self.missing, copy=0) #---- AtmQty Method: Does Quantities Data Structure Have A Quantity --- def has_qty(self, quantity): """Check if quantities data structure has quantity. Method Argument: * quantity: String scalar of the key name of the variable to check the quantities data structure for. See the class docstring for details on values of the quantities names. Returns 1 if quantities data structure self.qty has quantity; 0 if not. """ return self.qty.has_key(quantity) #-------- AtmQty Method: Check If Object is Newly Instantiated -------- def isnew(self): """Check if object is newly instantiated. Method Arguments: None. Checks whether the object is in a state consistent with being a newly instantiated object (returns 1 if is in such a state and 0 if not). "Newly instantiated" in this case means an AtmQty object that is initialized solely from an argument list, with no quantities added/deleted and attributes added/ deleted. Obviously, exhaustive tests cannot be done, but this method checks a number of conditions consistent with such a state. """ if len(self.qty) != 1: return 0 init_keys = [ 'lev_type', 'const', 'lev', 'missing', 'lat' \ , 'T_surf', 'qty_shape', 'lon', 'P_surf', 'qty'] if self.lev_type == "pressure": if len(self.__dict__) != 10: return 0 for attrib_key in init_keys: if not self.__dict__.has_key(attrib_key): return 0 if self.list_all_qty()[0] != 'p': return 0 if self.lev_type == "isentropic": if len(self.__dict__) != 11: return 0 for attrib_key in init_keys: if not self.__dict__.has_key(attrib_key): return 0 if not self.__dict__.has_key('p0'): return 0 if self.list_all_qty()[0] != 'theta': return 0 if self.P_surf != None: return 0 if self.T_surf != None: return 0 return 1 #----- AtmQty Method: List All Names In Quantities Data Structure ----- def list_all_qty(self): """List all quantities names in quantities data structure. Method Arguments: None. See the class docstring for details on values of the quantities names. """ return self.qty.keys() #-------------- AtmQty: Calculate Geopotential Altitude --------------- def alt(self, P_surf=None, T_surf=None): """Calculate geopotential altitude. Method Arguments: * P_surf: Surface pressure at latitude-longitude locations. Either a Numeric array of shape tuple (Numeric.size(self.lat), Numeric.size(self.lon)) or a scalar. Can have no missing values. Default is None, which means the value of self.P_surf is used. If self.P_surf is None, standard sea level pressure is used at all locations. * T_surf: Surface temperature at latitude-longitude locations. Either a Numeric array of shape tuple (Numeric.size(self.lat), Numeric.size(self.lon)) or a scalar. Can have no missing values. Default is None, which means the value of self.T_surf is used. If self.T_surf is None, standard sea level tempera- ture is used at all locations. Required Quantities in self.qty: p. Method adds/recalculates geopotential altitude to the self.qty quantities data structure (with key 'alt'), overwriting if the term previously exists. If 'p' does not already exist, that quantity is first calculated using the press method. Method overwrites self.P_surf and self.T_surf with whatever values are used in calculating altitude. """ from press2alt import press2alt if not self.has_qty('p'): self.press() if P_surf != None: self.P_surf = P_surf if T_surf != None: self.T_surf = T_surf if self.P_surf == None: self.P_surf = self.const.sea_level_press / 100. if self.T_surf == None: self.T_surf = self.const.sea_level_temp if N.size(self.P_surf) != 1: tmp_P_surf = N.zeros(self.qty_shape, typecode=N.Float) for ilev in xrange(self.qty_shape[2]): tmp_P_surf[:,:,ilev] = self.P_surf[:,:] else: tmp_P_surf = self.P_surf if N.size(self.T_surf) != 1: tmp_T_surf = N.zeros(self.qty_shape, typecode=N.Float) for ilev in xrange(self.qty_shape[2]): tmp_T_surf[:,:,ilev] = self.T_surf[:,:] else: tmp_T_surf = self.T_surf tmp = press2alt( self.get_qty('p'), P0=tmp_P_surf, T0=tmp_T_surf \ , missing=self.missing ) self.add_qty({'alt':tmp}) del tmp, tmp_P_surf, tmp_T_surf #------- AtmQty: Calculate Saturation Vapor Pressure (over ice) ------- def eice(self): """Calculate saturation vapor pressure over ice. Method Arguments: None. Required Quantities in self.qty: T. Method adds/recalculates saturation vapor pressure over ice to the self.qty quantities data structure (with key 'eice'), overwriting the term if it previously exists. """ from eice import eice self.add_qty({'eice':eice(self.get_qty('T'), missing=self.missing)}) #----- AtmQty: Calculate Saturation Vapor Pressure (over liquid) ------ def esat(self): """Calculate saturation vapor pressure over liquid water. Method Arguments: None. Required Quantities in self.qty: T. Method adds/recalculates saturation vapor pressure over liquid water to the self.qty quantities data structure (with key 'esat'), overwriting the term if it previously exists. """ from esat import esat self.add_qty({'esat':esat(self.get_qty('T'), missing=self.missing)}) #---------------- AtmQty: Calculate Coriolis Parameter ---------------- def fc(self): """Calculate Coriolis parameter. Method Arguments: None. Required Quantities in self.qty: None. Reference: Holton, p. 40. Method adds/recalculates Coriolis parameter to the self.qty quantities data structure (with key 'fc'), overwriting if the term if it previously exists. """ tmp = N.zeros(self.qty_shape, typecode=N.Float) deg2rad = N.pi / 180.0 for ilat in xrange(self.qty_shape[0]): tmp[ilat,:,:] = 2.0 * self.const.Omega \ * N.sin( self.lat[ilat] * deg2rad ) self.add_qty({'fc':tmp}) del tmp #---------- AtmQty: Calculate Isentropic Potential Vorticity ---------- def ipv(self): """Calculate isentropic potential vorticity. Method Arguments: None. Required Quantities in self.qty: fc, N2, rho, vort_rel Method adds/recalculates isentropic potential vorticity to self.qty (with key 'pv'), overwriting the term if it previ- ously exists. If 'fc', 'N2', 'rho', and 'vort_rel' do not already exist, this method first calculates these quantities using the applicable methods. Algorithm used to calculate pv is based on that used in the NCEP/NCAR reanalysis, with major differences. If lev_type is not isentropic level(s) and this method is called an error is produced. Reference: * Iredell, Mark: "How does NCEP/NCAR reanalaysis compute the potential vorticity on isentropic surfaces?" NCEP/NCAR Reanal- ysis FAQ: http://dss.ucar.edu/pub/reanalysis/FAQ.html. """ if self.lev_type != 'isentropic': raise ValueError, "AtmQty: Method not available" if not self.has_qty('fc'): self.fc() if not self.has_qty('N2'): self.static_stability() if not self.has_qty('rho'): self.rho() if not self.has_qty('vort_rel'): self.vort_rel() tmp = copy.deepcopy(self.get_qty_asMA('fc')) tmp += self.get_qty_asMA('vort_rel') tmp *= self.get_qty_asMA('N2') tmp /= self.const.g tmp /= self.get_qty_asMA('rho') tmp *= self.get_qty_asMA('theta') self.add_qty({'pv':MA.filled(tmp)}) del tmp #------------------- AtmQty: Calculate Mixing Ratio ------------------- def mix_ratio(self): """Calculate mixing ratio. Method Arguments: None. Required Quantities in self.qty: p and e, or q. Method adds/recalculates mixing ratio to the self.qty quan- tities dictionary (with key 'r'), overwriting the term if it previously exists. If 'q' does not already exist, then 'p' and 'e' are checked to see if either already exist; if either do not, this method first calculates each of them that do not exist using the applicable method. """ from mix_ratio import mix_ratio if not self.has_qty('q'): if not self.has_qty('p'): self.press() if not self.has_qty('e'): self.vapor_press() self.add_qty({'r':mix_ratio(self.qty, missing=self.missing)}) #--------------------- AtmQty: Calculate Pressure --------------------- def press(self): """Calculate pressure (from temperature, potential temperature). Method Arguments: None Required Quantities in self.qty: theta, T. Method adds/recalculates pressure to the self.qty quantities dictionary (with key 'p'), overwriting the term if it previous- ly exists. If lev_type is pressure level(s) and this method is called, an error is produced. """ if self.lev_type == 'pressure': raise ValueError, "AtmQty: Method not available" from press import press self.add_qty({'p':press(self.qty, p0=self.p0, missing=self.missing)}) #------------------- AtmQty: Calculate Air Density -------------------- def rho(self): """Calculate air density. Method Arguments: None. Required Quantities in self.qty: T_v, p. Method adds/recalculates air density (including moisture effects) to the self.qty quantities data structure (with key 'rho'), overwriting if the term previously exists. If 'T_v' does not already exist, this method first calculates it using the virt_temp method. """ from rho import rho if not self.has_qty('T_v'): self.virt_temp() if not self.has_qty('p'): self.press() tmp = rho( self.get_qty('p'), self.get_qty('T_v') \ , missing=self.missing ) self.add_qty({'rho':tmp}) del tmp #----------------- AtmQty: Calculate Static Stability ----------------- def static_stability(self): """Calculate static stability. Method Arguments: None. Required Quantities in self.qty: T, alt. Method adds/recalculates static stability to the self.qty quantities dictionary (with key 'N2'), overwriting if the term previously exists. If 'alt' does not already exist in qty, this method first calculates 'alt' using the method alt (with no keywords). """ from static_stability import static_stability if not self.has_qty('alt'): self.alt() tmp = N.zeros(self.qty_shape, typecode=N.Float) for ilat in range(N.size(self.lat)): for ilon in range(N.size(self.lon)): tmp[ilat,ilon,:] = \ static_stability( (self.get_qty('T'))[ilat,ilon,:] \ , (self.get_qty('alt'))[ilat,ilon,:] \ , missing=self.missing ) self.add_qty({'N2':tmp}) del tmp #--------------- AtmQty: Calculate Potential Temperature -------------- def theta(self, p0=1000.0): """Calculate potential temperature. Method Arguments: * p0: Reference pressure [hPa]. Floating point scalar. Optional (default set to 1000). Class attribute self.p0 is set to the value of this argument. Required Quantities in self.qty: p, T. Method adds/recalculates potential temperature to the self.qty quantities data structure (with key 'theta'), over- writing if the term if it previously exists. self.p0 is set to the value of argument p0. If lev_type is isentropic level(s) and this method is called, an error is produced. """ if self.lev_type == 'isentropic': raise ValueError, "AtmQty: Method not available" from theta import theta self.p0 = p0 tmp = theta( self.get_qty('p'), self.get_qty('T'), p0=self.p0 \ , missing=self.missing ) self.add_qty({'theta':tmp}) del tmp #------------- AtmQty: Calculate Vapor Pressure Over Water ------------ def vapor_press(self): """Calculate vapor pressure over liquid water. Method Arguments: None. Required Quantities in self.qty: p and r, q and r, or RH and esat. Method adds/recalculates vapor pressure to the self.qty quan- tities data structure (with key 'e'), overwriting the term if it previously exists. If 'esat' does not already exist in qty (and 'RH' exists, and 'r' and 'q' both do not exist), this method first calculates 'esat' using the method esat. If 'r' does not exist in qty, and 'RH' does not exist, this method first calculates 'r' using the method mix_ratio. """ from vapor_press import vapor_press if (not self.has_qty('r')) and (not self.has_qty('q')) and \ (not self.has_qty('esat')) and self.has_qty('RH'): self.esat() if (not self.has_qty('r')) and (not self.has_qty('RH')): self.mix_ratio() self.add_qty({'e':vapor_press(self.qty, missing=self.missing)}) #---------------- AtmQty: Calculate Virtual Temperature --------------- def virt_temp(self): """Calculate vapor pressure over liquid water. Method Arguments: None. Required Quantities in self.qty: r, T. Method adds/recalculates virtual temperature to the self.qty quantities data structure (with key 'T_v'), overwriting if the term previously exists. If 'r' does not already exist in self.qty, this method first calculates it using the method mix_ratio. """ from virt_temp import virt_temp if not self.has_qty('r'): self.mix_ratio() tmp = virt_temp( self.get_qty('r'), self.get_qty('T') \ , missing=self.missing ) self.add_qty({'T_v':tmp}) del tmp #---------------- AtmQty: Calculate Relative Vorticity ---------------- def vort_rel(self): """Calculate relative vorticity. Method Arguments: None. Required Quantities in self.qty: alt, u, v. Method adds/recalculates relative vorticity to the self.qty quantities data structure (with key 'vort_rel'), overwriting if the term previously exists. The 'default_spherical' al- gorithm from gemath.curl_2d is used. If 'alt' does not al- ready exist in qty, this method first calculates 'alt' using the method alt (with no keywords). Method assumes we're on the Earth and that the mean radius of the earth coincides with an atmospheric altitude of 0. Note that if quantity 'alt' contains any missing values, for the purposes of this current method only, the altitude at those missing value locations is set to 0. """ from gemath.curl_2d import curl_2d from gemath.has_close import has_close from gemath.where_close import where_close if not self.has_qty('alt'): self.alt() if has_close(self.get_qty('alt'), self.missing): missloc = where_close(self.get_qty('alt'), self.missing).flat thisalt = self.get_qty('alt').copy().flat N.putmask(thisalt, missloc, 0.0) thisalt = N.reshape(thisalt, self.qty_shape) del missloc else: thisalt = self.get_qty('alt') tmp = N.zeros( self.qty_shape, typecode=N.Float ) for ilev in xrange(len(self.lev)): tmpalt = self.const.Rad_earth + thisalt[:,:,ilev] tmp[:,:,ilev] = curl_2d( self.lon, self.lat \ , self.get_qty('u')[:,:,ilev] \ , self.get_qty('v')[:,:,ilev] \ , missing=self.missing \ , algorithm="default_spherical" \ , R_sphere=tmpalt ) self.add_qty({'vort_rel':tmp}) del thisalt, tmp, tmpalt #-------------------------- Main: Test Module ------------------------- #- Define additional examples for doctest to use: __test__ = { 'Additional Example 1: Constants': """ >>> from atmconst import AtmConst >>> const = AtmConst() >>> const.kappa 0.28590637450199202 >>> const.epsilon 0.62195306913960091 """, 'Additional Example 2: Homogenous x-y domain AtmQty object': """ >>> import atmqty >>> import Numeric as N >>> lat = N.array([-45., 0., 45.]) >>> lon = N.array([0., 90., 180., 270.]) >>> lev = N.array([1000., 850., 500.]) >>> arrayshape = ( (len(lat), len(lon), len(lev)) ) >>> T = N.zeros(arrayshape, typecode=N.Float) >>> q = N.zeros(arrayshape, typecode=N.Float) >>> for j in range(len(lon)): ... for i in range(len(lat)): ... T[i,j,:] = [273.1, 257.3, 233.1] ... q[i,j,:] = [5.5e-3, 3.4e-3, 2.3e-3] >>> x = atmqty.AtmQty(lat, lon, lev) >>> ['%.6f' % x.lat[i] for i in range(3)] ['-45.000000', '0.000000', '45.000000'] >>> x.add_qty({'q':q}) >>> x.add_qty({'T':T}) >>> x.theta() >>> ['%.6f' % x.get_qty('theta')[1,2,i] for i in range(3)] ['273.100000', '269.537605', '284.189919'] >>> x.vapor_press() >>> ['%.6f' % x.get_qty('e')[0,0,i] for i in range(3)] ['8.813646', '4.637070', '1.846433'] """, 'Additional Example 3: More Example 3 methods, with missing values': """ >>> import atmqty >>> import Numeric as N >>> lat = N.arange(13) * 15.0 - 90.0 >>> lon = N.arange(25) * 15.0 - 180.0 >>> lev = N.array([850.]) >>> x = atmqty.AtmQty(lat, lon, lev, missing=1e29) >>> base = N.outerproduct(N.sin(lat*N.pi/360.), N.cos(lon*N.pi/360.)) >>> base = N.reshape(base, (base.shape[0], base.shape[1], 1)) >>> T = (base * 10.) + 257.3 >>> q = N.absolute(base) * 5e-3 >>> x.add_qty( {'q':q}, {'T':T} ) >>> x.get_qty('T')[1,2:8,0] = 1e29 >>> x.get_qty('q')[2,2:8,0] = 1e29 >>> ['%.7g' % x.get_qty('T')[i,3,0] for i in range(5)] ['254.594', '1e+29', '255.3866', '255.8355', '256.3095'] >>> print x.P_surf, x.T_surf None None >>> x.alt() >>> ['%.7g' % x.get_qty('alt')[i,3,0] for i in range(5)] ['1457.285', '1457.285', '1457.285', '1457.285', '1457.285'] >>> print x.P_surf, x.T_surf 1013.25 288.15 >>> x.alt(P_surf=1000.) >>> ['%.7g' % x.get_qty('alt')[i,3,0] for i in range(5)] ['1349.778', '1349.778', '1349.778', '1349.778', '1349.778'] >>> print x.P_surf, x.T_surf 1000.0 288.15 >>> x.P_surf = 1013.25 >>> x.alt() >>> ['%.7g' % x.get_qty('alt')[i,3,0] for i in range(5)] ['1457.285', '1457.285', '1457.285', '1457.285', '1457.285'] >>> print x.P_surf, x.T_surf 1013.25 288.15 >>> tmp_Ps = N.zeros(x.qty_shape[0:2], typecode=N.Float) + 1000. >>> tmp_Ts = N.zeros(x.qty_shape[0:2], typecode=N.Float) + 288.15 >>> x.alt(P_surf=tmp_Ps, T_surf=tmp_Ts) >>> ['%.7g' % x.get_qty('alt')[i,3,0] for i in range(5)] ['1349.778', '1349.778', '1349.778', '1349.778', '1349.778'] >>> x.esat() >>> ['%.7g' % x.get_qty('esat')[i,3,0] for i in range(5)] ['1.420472', '1e+29', '1.519123', '1.577691', '1.641716'] >>> x.mix_ratio() >>> ['%.7g' % x.get_qty('r')[i,3,0] for i in range(5)] ['0.001354823', '0.001166173', '1e+29', '0.0007327696', '0.0004954742'] >>> x.theta() >>> ['%.9g' % x.get_qty('theta')[i,2,0] for i in range(3)] ['267.620434', '1e+29', '268.181961'] >>> x.vapor_press() >>> ['%.8g' % x.get_qty('e')[i,3,0] for i in range(5)] ['1.8475617', '1.5907822', '1e+29', '1.0002703', '0.67660701'] >>> x.virt_temp() >>> ['%.7g' % x.get_qty('T_v')[i,3,0] for i in range(5)] ['254.8034', '1e+29', '1e+29', '255.9494', '256.3867'] """, 'Additional Example 4: 3-D pressure domain, with missing values': """ >>> import atmqty >>> import Numeric as N >>> lat = N.arange(13) * 15.0 - 90.0 >>> lon = N.arange(25) * 15.0 - 180.0 >>> lev = N.array([850., 500., 250.]) #- Test error check domain on AtmQty object initialization: >>> xb = atmqty.AtmQty(lat, lon, N.array([850., 250., 500.])) Traceback (most recent call last): ... ValueError: AtmQty: Levels not monotonic >>> xb = atmqty.AtmQty(lat[::-1], lon, lev, missing=1e29) Traceback (most recent call last): ... ValueError: AtmQty: Input(s) not monotonic increasing >>> xb = atmqty.AtmQty(lat, lon[::-1], lev, missing=1e29) Traceback (most recent call last): ... ValueError: AtmQty: Input(s) not monotonic increasing #- Initialize AtmQty object and test method isnew: >>> x = atmqty.AtmQty(lat, lon, lev, missing=1e29) >>> x.isnew() 1 >>> ['%.7g' % x.qty_shape[i] for i in range(len(x.qty_shape))] ['13', '25', '3'] >>> [x.list_all_qty()[i] for i in range(len(x.list_all_qty()))] ['p'] >>> base = N.outerproduct(N.sin(lat*N.pi/360.), N.cos(lon*N.pi/360.)) >>> T = N.zeros(x.qty_shape, typecode=N.Float) >>> q = N.zeros(x.qty_shape, typecode=N.Float) >>> T[:,:,0] = (base * 10.) + 289.3 >>> T[:,:,1] = T[:,:,0] - 20.2 >>> T[:,:,2] = T[:,:,1] - 40.8 >>> q[:,:,0] = (N.absolute(base) * 5e-3)[:,:] >>> q[:,:,1] = (N.absolute(base) * 2e-3)[:,:] >>> q[:,:,2] = (N.absolute(base) * 1e-5)[:,:] >>> x.add_qty( {'q':q}, {'T':T} ) >>> x.get_qty('T')[1,2:8,0] = 1e29 >>> x.get_qty('q')[2,2:8,0] = 1e29 >>> x.isnew() 0 #- Get values of T and q at certain points for reference: >>> ['%.7g' % x.get_qty('T')[1,5,i] for i in range(len(lev))] ['1e+29', '265.3941', '224.5941'] >>> ['%.7g' % x.get_qty('T')[3,15,i] for i in range(len(lev))] ['285.7645', '265.5645', '224.7645'] >>> ['%.7g' % x.get_qty('q')[2,6,i] for i in range(len(lev))] ['1e+29', '0.0007071068', '3.535534e-06'] #- Test method alt: >>> x.alt() >>> ['%.7g' % x.get_qty('alt')[1,5,i] for i in range(len(lev))] ['1457.285', '5574.382', '10362.85'] >>> ['%.7g' % x.get_qty('alt')[3,15,i] for i in range(len(lev))] ['1457.285', '5574.382', '10362.85'] #- Test method static_stability: >>> x.static_stability() >>> ['%.7g' % x.get_qty('N2')[1,5,i] for i in range(len(lev))] ['1e+29', '1e+29', '5.445355e-05'] >>> ['%.7g' % x.get_qty('N2')[3,15,i] for i in range(len(lev))] ['0.0001668233', '0.000107752', '5.441227e-05'] #- Test method rho and automatic calculation of dependency in rho: >>> 'r' in x.list_all_qty() 0 >>> 'T_v' in x.list_all_qty() 0 >>> x.rho() >>> ['%.7g' % x.get_qty('rho')[1,5,i] for i in range(len(lev))] ['1e+29', '0.6560328', '0.387778'] >>> ['%.7g' % x.get_qty('rho')[3,15,i] for i in range(len(lev))] ['1.035111', '0.6556255', '0.3874841'] >>> 'r' in x.list_all_qty() 1 >>> 'T_v' in x.list_all_qty() 1 >>> ['%.7g' % x.get_qty('p')[3,15,i] for i in range(len(lev))] ['850', '500', '250'] >>> ['%.7g' % x.get_qty('T_v')[3,15,i] for i in range(len(lev))] ['286.0715', '265.6786', '224.7649'] #- Test method get_qty_asMA: >>> ['%.7g' % x.get_qty('rho')[1,5,i] for i in range(len(lev))] ['1e+29', '0.6560328', '0.387778'] >>> tmp = x.get_qty_asMA('rho') >>> ['%.7g' % MA.filled(tmp)[1,5,i] for i in range(len(lev))] ['1e+29', '0.6560328', '0.387778'] >>> ['%.7g' % tmp.mask()[1,5,i] for i in range(len(lev))] ['1', '0', '0'] #- Test method fc: >>> x.fc() >>> ['%.7g' % x.get_qty('fc')[1,5,i] for i in range(len(lev))] ['-0.0001408726', '-0.0001408726', '-0.0001408726'] >>> ['%.7g' % x.get_qty('fc')[3,15,i] for i in range(len(lev))] ['-0.0001031259', '-0.0001031259', '-0.0001031259'] """, 'Additional Example 5: 3-D isentropic domains, with missing values': """ >>> import atmqty >>> import Numeric as N >>> lat = N.arange(13) * 15.0 - 90.0 >>> lon = N.arange(25) * 15.0 - 180.0 >>> lev = N.array([300., 315., 330.]) >>> x = atmqty.AtmQty(lat, lon, lev, lev_type="isentropic", missing=1e29) >>> y = atmqty.AtmQty(lat, lon, lev, lev_type="isentropic", p0=1010.) >>> x.isnew() 1 >>> x.has_qty('theta') 1 >>> x.has_qty('p') 0 >>> y.isnew() 1 >>> x.p0 1000.0 >>> y.p0 1010.0 >>> x.missing 9.9999999999999991e+28 >>> y.missing 1e+20 >>> base = N.outerproduct(N.sin(lat*N.pi/360.), N.cos(lon*N.pi/360.)) >>> T = N.zeros(x.qty_shape, typecode=N.Float) >>> q = N.zeros(x.qty_shape, typecode=N.Float) >>> T[:,:,0] = (base * 10.) + 289.3 >>> T[:,:,1] = T[:,:,0] - 20.2 >>> T[:,:,2] = T[:,:,1] - 40.8 >>> q[:,:,0] = (N.absolute(base) * 5e-3)[:,:] >>> q[:,:,1] = (N.absolute(base) * 2e-3)[:,:] >>> q[:,:,2] = (N.absolute(base) * 1e-5)[:,:] >>> x.add_qty( {'q':q}, {'T':T} ) >>> x.get_qty('T')[1,2:8,0] = 1e29 >>> x.get_qty('q')[2,2:8,0] = 1e29 >>> x.isnew() 0 >>> y.isnew() 1 >>> ['%.9g' % x.get_qty('T')[3,15,i] for i in range(len(lev))] ['285.764466', '265.564466', '224.764466'] >>> ['%.7g' % x.get_qty('theta')[3,15,i] for i in range(len(lev))] ['300', '315', '330'] >>> x.press() >>> ['%.7g' % x.get_qty('p')[3,15,i] for i in range(len(lev))] ['843.6339', '550.4049', '260.9998'] """, 'Additional Example 6: 3-D isentropic domains, calculate IPV, vorticity': """ #- Imports and establish AtmQty object: >>> import atmqty >>> import Numeric as N >>> nx = 72 >>> ny = 37 >>> lon = N.arange(nx) * 5.0 - 180.0 >>> lat = N.arange(ny) * 5.0 - 90.0 >>> lev = N.array([300., 315., 330.]) >>> x = atmqty.AtmQty( lat, lon, lev, lev_type="isentropic", missing=1e+29 ) #- Make-up data (T, u, v have values; q is zero everywhere): >>> lat_2d = N.reshape( N.repeat(lat, nx), (ny, nx) ) >>> lon_2d = N.reshape( N.repeat(N.reshape(lon,(1,nx)), ny), (ny, nx) ) >>> u_2d = N.sin(lon_2d*N.pi/180.)*N.sin((lat_2d-90.)*N.pi/30.) >>> v_2d = N.sin((lon_2d-30.)*N.pi/180.)*N.sin((lat_2d-90.)*N.pi/30.) >>> N.putmask( u_2d, N.where(N.absolute(u_2d) < 1e-10, 1, 0), 0. ) >>> N.putmask( v_2d, N.where(N.absolute(v_2d) < 1e-10, 1, 0), 0. ) >>> T_2d = N.outerproduct( N.sin(lat*N.pi/360.) \ , N.cos(lon*N.pi/360.) ) * 10. + 257.3 >>> T = N.zeros(x.qty_shape, typecode=N.Float) >>> u = N.zeros(x.qty_shape, typecode=N.Float) >>> v = N.zeros(x.qty_shape, typecode=N.Float) >>> q = N.zeros(x.qty_shape, typecode=N.Float) >>> T[:,:,0] = T_2d >>> u[:,:,0] = u_2d >>> v[:,:,0] = v_2d >>> for ilev in range(1,len(lev)): ... T[:,:,ilev] = T[:,:,ilev-1] - N.sin(ilev*4.0/N.pi)*5.0 ... u[:,:,ilev] = u[:,:,ilev-1] + N.sin(ilev*4.0/N.pi)*4.0 + 3.0 ... v[:,:,ilev] = v[:,:,ilev-1] + N.sin(ilev*4.0/N.pi)*2.0 #- Set missing values for u and v (cannot have missing values in T # otherwise there will be missing values in altitude) and load data: >>> u[10:18,23:32,0] = 1e+29 >>> v[17:24,62:70,2] = 1e+29 >>> x.add_qty({'T':T, 'u':u, 'v':v, 'q':q}) #- Calculate relative vorticity: >>> x.vort_rel() >>> ['%.9g' % x.get_qty('vort_rel')[3,15,i] for i in range(len(lev))] ['6.10689002e-08', '-4.34871759e-06', '-6.98941959e-06'] >>> ['%.9g' % x.get_qty('vort_rel')[13,24,i] for i in range(len(lev))] ['1e+29', '2.5439683e-07', '-1.82497184e-07'] #- Calculate IPV: >>> x.ipv() >>> ['%.9g' % x.get_qty('pv')[3,15,i] for i in range(len(lev))] ['-1.5339659e-06', '-2.21039439e-06', '-3.05348599e-06'] >>> ['%.9g' % x.get_qty('pv')[13,24,i] for i in range(len(lev))] ['1e+29', '-9.11526793e-07', '-1.24498967e-06'] """} #- Execute doctest if module is run from command line: if __name__ == "__main__": """Test the module. Tests the examples in all the module documentation strings, plus __test__. The parent directory to the current directory is added to sys.path. """ import doctest, sys, os sys.path.append(os.pardir) doctest.testmod(sys.modules[__name__]) # ===== end file =====