| | |
- curl_2d(x, y, Fx, Fy, missing=1e+20, algorithm='default', R_sphere=6371220.0)
- Curl of a vector F on a 2-D "rectangular" grid.
The 2-D grid F is defined on is rectangular, meaning that while
the grid spacings do not have to be even, each grid box is
rectangular and grid boundaries are rectilinear and parallel to
each other. If the 2-D grid is on a sphere, we assume that lines
of constant longitude and latitude form grids that are rectangular
(though in reality this is not strictly true).
Since F does not have a z-component, the curl of F is positive
pointing as a normal out of the plane defined by the basis for
x and y.
Positional Input Arguments:
* x: x-coordinate of each F. 1-D Numeric array with number of
elements equal to the number of columns in array F. Typically
x is in km, or on a sphere the longitude in deg. Floating or
integer type.
* y: y-coordinate of each F. 1-D Numeric array with number of
elements equal to the number of rows in array F. Typically
y is in km, or on a sphere the latitude in deg. Floating or
integer type.
* Fx: x-component of vector quantity F. 2-D Numeric array of
shape (len(y), len(x)). Floating or integer type.
* Fy: y-component of vector quantity F. 2-D Numeric array of
same shape and size as Fx. Floating or integer type.
Keyword Input Arguments:
* missing: If Fx and/or Fy has missing values, this is the
missing value value. Scalar. Default is 1e+20.
* algorithm: Name of the algorithm to use. String scalar.
Default is 'default'. Possible values include:
+ 'default': Default method. Algorithm 'order1_cartesian'
is used.
+ 'default_spherical': If 'spherepack' can be used, chooses
that method. Otherwise 'order1_spherical' is used. Either
way this option assumes the domain is on a sphere of radius
equal to keyword R_sphere.
+ 'order1_cartesian': First-order finite-differencing (back-
ward and forward differencing used at the endpoints, and
centered differencing used everywhere else). If abscissa
intervals are irregular, differencing will be correspon-
dingly asymmetric. Grid is assumed Cartesian. The curl
returned is in units of F divided by the units of x and y
(x and y must be in the same units).
+ 'order1_spherical': First-order finite-differencing (back-
ward and forward differencing used at the endpoints, and
centered differencing used everywhere else). If abscissa
intervals are irregular, differencing will be correspon-
dingly asymmetric. Grid is assumed spherical and x and y
are longitude and latitude (in deg). The curl returned is
in units of F divided by units of R_sphere. Note because
of singularities in the algorithm near the poles, the
curl for all points north of 88N or south of 88S latitude
are set to missing.
+ 'spherepack': The NCAR package SPHEREPACK 3.0 is used for
calculating the curl. The spatial domain is the global
sphere; x and y are longitude and latitude (in deg),
respectively; x must be evenly spaced and y must be gauss-
ian or evenly spaced. The domain must cover the entire
sphere and x and y must be longitude and latitude, respec-
tively, in deg. There can be no missing values. The
algorithm also assumes that x and y are monotonically in-
creasing from index 0. The curl returned is in units of F
divided by m. Thus, if F is velocity in m/s, the curl of
F is in units 1/s (and is the vorticity). Note that in
this function's implementation, the final curl output using
algorithm 'spherepack' is adjusted for R_sphere not equal
to the default mean earth radius, if appropriate, so this
algorithm can be used on other planets.
Note that SPHEREPACK operates on single precision floating
point (Float32) values. This function silently converts
input to Float32 for the computations, as needed (the input
parameters are not altered, however).
* R_sphere: If the grid is the surface of a sphere, this is
the radius of the sphere. This keyword is only used when
the algorithm is 'order1_spherical' or 'spherepack' (either
chosen explicitly by the algorithm keyword or automatically
when 'default_spherical' is used). Default value is the
"mean" radius of the Earth, in meters. "Mean" is in quotes
because this value changes depending on what you consider the
mean; the default in this function is the same Earth radius
used in module sphere. Note that if the level is a distance
z above the "mean" radius of the earth, R_sphere should be
set to the "mean" radius plus z.
Value can be a scalar or a Numeric or MA array of the same
size and shape as Fx. If keyword is an array, the value of
R_sphere for each element corresponds to the same location
as in the Fx and Fy arrays; there should also be no missing
values in R_sphere (function doesn't check for this, however).
Output:
* Curl of F. Units depend on algorithm used to calculate curl
(see above discussion of keyword algorithm). Numeric array of
same shape as Fx and Fy inputs.
If there are missing values, those elements in the output that
used missing values in the calculation of the curl are set to
the value in keyword missing. If there are missing values in
the output due to math errors and keyword missing is set to
None, output will fill those missing values with the MA default
value of 1e+20.
References:
* Glickman, T. S. (Ed.) (2000): "Curl," Glossary of Meteorology,
Boston, MA: American Meteorological Society, ISBN 1-878220-
34-9.
* Haltiner, G. J, and R. T. Williams (1980): Numerical Prediction
and Dynamic Meteorology, New York: John Wiley & Sons, pp. 8-9.
* Holton, J. R. (1992): An Introduction to Dynamic Meteorology,
San Diego, CA: Academic Press, p. 482.
* Kauffman, B. G., and W. G. Large (2002): The CCSM Coupler,
Version 5.0, User's Guide, Source Code Reference, and Scientific
Description. Boulder, CO: NCAR. Value for earth radius is
given at URL:
http://www.ccsm.ucar.edu/models/ccsm2.0/cpl5/users_guide/node10.html
* Overview of the CDAT Interface to the NCAR SPHEREPACK 3.0.
SPHEREPACK is by J. C. Adams and P. N. Swarztrauber. More
information on SPHEREPACK can be found at:
http://www.scd.ucar.edu/css/software/spherepack/
while information on the CDAT implementation used here is at:
http://esg.llnl.gov/cdat/getting_started/cdat.htm
Example of a synthetic dataset of winds. The curl is the vorti-
city. Images of the wind field and the vorticity are online at
http://www.johnny-lin.com/py_pkgs/gemath/doc/test_curl_2d_add.html.
Velocity at y greater than 60 and less than -60 are set to 0, to
prevent spectral algorithms from giving errors at the poles. Tiny
values less than 1e-10 are also set to 0:
(a) Import statements and create data:
>>> import Numeric as N
>>> from curl_2d import curl_2d
>>> nx = 181
>>> ny = 91
>>> x = N.arange(nx) * 2.0 - 180.0
>>> y = N.arange(ny) * 2.0 - 90.0
>>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )
>>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )
>>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)
>>> N.putmask( u, N.where(N.absolute(u) < 1e-10, 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(v) < 1e-10, 1, 0), 0. )
>>> N.putmask( u, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
>>> N.putmask( v, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )
(b) Use order 1 Cartesian algorithm: If x and y are in meters,
and u and v are in m/s, the curl is in 1/s:
>>> curl = curl_2d(x, y, u, v, algorithm='order1_cartesian')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-0.007251593', '-0.003628007', '0', '0.003628007']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['0', '0', '0', '0']
(c) Use spectral method from SPHEREPACK (need to first remove the
repeating dateline point from the domain). Here we take x and
y in deg, while u and v are in m/s. If so the curl is in 1/s.
These results as much less than for example (b) above since
the domain is much larger:
>>> x = x[0:-1]
>>> u = u[:,0:-1]
>>> v = v[:,0:-1]
>>> curl = curl_2d(x, y, u, v, algorithm='spherepack')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.436402e-08', '-3.220152e-08', '1.33852e-13', '3.220165e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['-2.299736e-20', '-4.806512e-21', '-9.283145e-22', '-1.368445e-20']
(d) Use order 1 spherical algorithm: x and y are in deg and u and
v are in m/s. The curl is in 1/s. Note that these results are
nearly identical to those from SPHEREPACK, except near the poles
the values are set to missing:
>>> curl = curl_2d(x, y, u, v, algorithm='order1_spherical')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-6.517317e-08', '-3.260645e-08', '0', '3.260645e-08']
>>> ['%.7g' % curl[0,i] for i in range(8,12)]
['1e+20', '1e+20', '1e+20', '1e+20']
(e) Use "default" algorithm, which for this domain chooses the
order 1 cartesian algorithm:
>>> curl = curl_2d(x, y, u, v, algorithm='default')
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-0.007251593', '-0.003628007', '0', '0.003628007']
(f) Use "default_spherical" algorithm with missing data (equal to
1e+20), which because of the missing data will choose the order
1 spherical algorithm. We also do the calculation on Mars, and
illustrate passing in the radius as a constant as well as an
array of the same size and shape as u and v:
>>> u[30:46,90:114] = 1e+20
>>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' , R_sphere=3390e+3)
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.224875e-07', '-6.128107e-08', '1e+20', '1e+20']
>>> R_mars = N.zeros(u.shape, typecode=N.Float) + 3390e+3
>>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' , R_sphere=R_mars )
>>> ['%.7g' % curl[45,i] for i in range(88,92)]
['-1.224875e-07', '-6.128107e-08', '1e+20', '1e+20']
|