source: branches/numpy/anuga/utilities/interp.py @ 7208

Last change on this file since 7208 was 7208, checked in by rwilson, 14 years ago

Fixed problems with back-merge.

File size: 10.4 KB
Line 
1#!/usr/bin/python -tt
2#=======================================================================
3#                        General Documentation
4
5"""Single-function module.
6
7   See function docstring for description.
8"""
9
10#-----------------------------------------------------------------------
11#                       Additional Documentation
12#
13# RCS Revision Code:
14#   $Id: interp.py,v 1.2 2004/03/23 04:28:16 jlin Exp $
15#
16# Modification History:
17# - 22 Mar 2004:  Original by Johnny Lin, Computation Institute,
18#   University of Chicago.  Passed passably reasonable tests.
19#
20# Notes:
21# - Written for Python 2.2.
22# - Module docstrings can be tested using the doctest module.  To
23#   test, execute "python interp.py".
24# - See import statements throughout for non-"built-in" packages and
25#   modules required.
26#
27# Copyright (c) 2004 by Johnny Lin.  For licensing, distribution
28# conditions, contact information, and additional documentation see
29# the URL http://www.johnny-lin.com/py_pkgs/gemath/doc/;
30
31# This library is free software; you can redistribute it and/or modify
32# it under the terms of the GNU Lesser General Public License as
33# published by the Free Software Foundation; either version 2.1 of the
34# License, or (at your option) any later version.
35
36# This library is distributed in the hope that it will be useful, but
37# WITHOUT ANY WARRANTY; without even the implied warranty of
38# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
39# Lesser General Public License for more details.
40
41# You should have received a copy of the GNU Lesser General Public
42# License along with this library; if not, write to the Free Software
43# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
44# USA.
45
46# You can contact Johnny Lin at his email address or at the University
47# of Chicago, Department of the Geophysical Sciences, 5734 S. Ellis
48# Ave., Chicago, IL 60637, USA.
49
50
51#=======================================================================
52
53
54
55
56#----------------------- Overall Module Imports ------------------------
57
58#- Set module version to package version:
59
60
61
62
63#------------------------ Non-Private Function -------------------------
64
65def interp(y, x, xinterp, missing=1e+20):
66    """Simple linear interpolation for ordinate with missing values.
67
68
69    Vectors x and y are the data describing a piecewise linear function.
70    Function returns the interpolated values of the ordinate function
71    at abscissa values in xinterp.  Values of xinterp outside the range
72    of x are returned as missing.  Any elements in the output that uses
73    missing values in y for the interpolation are set to missing.
74
75
76    Positional Input Arguments:
77    * y:  Ordinate values of data.  Rank 1 numeric vector.  Required.
78      Can have missing values.  Floating or integer type.
79
80    * x:  Abscissa values of data.  Rank 1 numeric vector.  Required.
81      Can have no missing values.  Must be monotonically ascending. 
82      Floating or integer type.
83
84    * xinterp:  Abscissa values to calculate interpolated ordinate
85      values at.  Rank 1 numeric vector or numeric scalar.  Required. 
86      Can have no missing values.  Can be in any order.  Floating or
87      integer type.
88
89
90    Keyword Input Arguments:
91    * missing:  If input has missing values, this is the missing value
92      value.  Scalar.  Floating or integer type.  Default is 1e+20.
93
94
95    Output Result:
96    * Interpolated ordinate values at xinterp.  Rank 1 numeric vector
97      of same length as xinterp (if xinterp is a numeric scalar,
98      output is also a numeric scalar).  Missing values are set to the
99      value of argument missing.  Type is Float, even if argument
100      missing and inputs are all integer.
101
102
103    References:
104    * Lin, J. W.-B.:  "Simple Interpolation."
105      Python/CDAT for Earth Scientists: Tips and Examples.
106      http://www.johnny-lin.com/cdat_tips/tips_math/interp.html
107
108
109    Example with no missing values (gives same output as function
110    arrayfns.interp):
111
112    >>> from interp import interp
113    >>> import numpy as N
114    >>> x = N.array([1., 2., 3., 4., 5.])
115    >>> y = N.array([3., 6., 2.,-5.,-3.])
116    >>> xint = N.array([3.4, 2.3])
117    >>> yint = interp(y, x, xint, missing=1e+20)
118    >>> ['%.7g' % yint[i] for i in range(len(yint))]
119    ['-0.8', '4.8']
120
121    Example with missing values:
122
123    >>> x = N.array([1.,    2., 3.,  4.,  5.])
124    >>> y = N.array([3., 1e+20, 2., -5., -3.])
125    >>> xint = N.array([3.4, 2.3])
126    >>> yint = interp(y, x, xint, missing=1e+20)
127    >>> ['%.7g' % yint[i] for i in range(len(yint))]
128    ['-0.8', '1e+20']
129
130    Example with values out of range of the data:
131
132    >>> x = N.array([1.,   2.1, 3.,  4., 5.1])
133    >>> y = N.array([3., 1e+20, 2., -5., -3.])
134    >>> xint = N.array([3.4, -2.3, 6.])
135    >>> yint = interp(y, x, xint, missing=1e+20)
136    >>> ['%.7g' % yint[i] for i in range(len(yint))]
137    ['-0.8', '1e+20', '1e+20']
138    """
139    import arrayfns
140    import numpy.ma MA
141    import numpy as N
142    from where_close import where_close
143
144
145    #- Check inputs for possible errors:
146
147    if (N.rank(y) != 1) or (N.rank(x) != 1):
148        raise ValueError, "interp:  Input(s) not a vector"
149    if N.rank(xinterp) > 1:
150        raise ValueError, "interp:  xinterp not a vector or scalar"
151    if x[-1] <= x[0]:
152        raise ValueError, "interp:  x not monotonically increasing"
153
154
155    #- Establish constants and define xint, a rank 1 version of
156    #  xinterp to be used for the rest of the function:
157
158    if N.rank(xinterp) == 0:
159        xint = N.reshape(xinterp, (1,))
160    else:
161        xint = xinterp
162
163    num_xint = N.size(xint)
164
165
166    #- Mask as missing values of xint that are outside of the range
167    #  of x:
168
169    yint_outrange_mask = N.logical_or( N.less(xint, x[0]) \
170                                     , N.greater(xint, x[-1]) )
171
172
173    #- Mask of elements with missing values in y, if there are any
174    #  missing values in y.  If xint equals a value in x, missing
175    #  values mask for that xint is the same as the corresponding
176    #  value in x; and mask elements in xint which fall in an interval
177    #  (whose upper bound index is top_idx) where one of the endpoints
178    #  is missing:
179
180    y_miss_mask    = where_close(y, missing)
181    yint_miss_mask = N.zeros(num_xint)
182
183    if MA.maximum(y_miss_mask) == 1:
184
185        for i in xrange(num_xint):
186            if yint_outrange_mask[i] == 0:
187                x_eq_xint = where_close(x, xint[i])
188                if MA.maximum(x_eq_xint) == 1:
189                    yint_miss_mask[i] = y_miss_mask[N.nonzero(x_eq_xint)]
190                else:
191                    top_idx = N.nonzero(N.greater(x, xint[i]))[0]
192                    yint_miss_mask[i] = y_miss_mask[top_idx] or \
193                                        y_miss_mask[top_idx-1]
194
195
196    #- Return interpolated values, set to missing values as
197    #  appropriate, and making a scalar if xinterp is a scalar:
198
199    yint = arrayfns.interp(y, x, xint)
200    N.putmask( yint, N.logical_or(yint_miss_mask, yint_outrange_mask) \
201             , missing)
202    if N.rank(xinterp) == 0:  yint = yint[0]
203
204    return yint
205
206
207
208
209#-------------------------- Main:  Test Module -------------------------
210
211#- Define additional examples for doctest to use:
212
213__test__ = {'Additional Examples':
214    """
215    (1) General error catching:
216
217    >>> from interp import interp
218    >>> import numpy as N
219    >>> x = N.array([1.,    2., 3.,  4.,  5.,  6.])
220    >>> y = N.array([3., 1e+20, 2., -5., -3., -4.])
221    >>> x = N.reshape(x, (2,3))
222    >>> y = N.reshape(y, (2,3))
223    >>> xint = N.array([3.4, 2.3])
224    >>> yint = interp(y, x, xint, missing=1e+20)
225    Traceback (most recent call last):
226        ...
227    ValueError: interp:  Input(s) not a vector
228
229    >>> x = N.array([1.,    2., 3.,  4.,  5.,  6.])
230    >>> y = N.array([3., 1e+20, 2., -5., -3., -4.])
231    >>> xint = N.array([[3.4, 2.3],[3.4, 2.3]])
232    >>> yint = interp(y, x, xint, missing=1e+20)
233    Traceback (most recent call last):
234        ...
235    ValueError: interp:  xinterp not a vector or scalar
236
237    >>> x = N.array([1.,    2., 3.,  4.,  5.,  0.])
238    >>> y = N.array([3., 1e+20, 2., -5., -3., -4.])
239    >>> xint = N.array([3.4, 2.3])
240    >>> yint = interp(y, x, xint, missing=1e+20)
241    Traceback (most recent call last):
242        ...
243    ValueError: interp:  x not monotonically increasing
244
245    >>> x = N.array([1.,    2., 3.,  4.,  5.,  6.])
246    >>> y = N.array([3., None, 2., -5., -3., -4.])
247    >>> xint = N.array([3.4, 2.3, 2., 5., 3., 1.])
248    >>> yint = interp(y, x, xint, missing=None)
249    Traceback (most recent call last):
250        ...
251    ValueError: where_close:  Inputs must be Float or Integer
252
253    (2) Values right on the border of intervals:
254
255    >>> x = N.array([1.,    2., 3.,  4.,  5.,  6.])
256    >>> y = N.array([3., 1e+20, 2., -5., -3., -4.])
257    >>> xint = N.array([3.4, 2.3, 2., 5., 3., 1.])
258    >>> yint = interp(y, x, xint, missing=1e+20)
259    >>> ['%.7g' % yint[i] for i in range(len(yint))]
260    ['-0.8', '1e+20', '1e+20', '-3', '2', '3']
261
262    (3) Single element vector input:
263
264    >>> yint = interp(y, x, N.array([6.]), missing=1e+20)
265    >>> ['%.7g' % yint[i] for i in range(len(yint))]
266    ['-4']
267
268    (4) Scalar xint:
269
270    >>> x = N.array([1.,    2., 3.,  4.,  5.,  6.])
271    >>> y = N.array([3., 1e+20, 2., -5., -3., -4.])
272    >>> yint = interp(y, x, N.array(6.), missing=1e+20)
273    >>> yint
274    -4.0
275    >>> N.rank(yint)
276    0
277
278    (5) Integer values:
279
280    >>> x = N.arange(6)
281    >>> y = N.arange(6)
282    >>> xint = N.array([3.4, 2.3])
283    >>> yint = interp(y, x, xint, missing=-9999999)
284    >>> ['%.7g' % yint[i] for i in range(len(yint))]
285    ['3.4', '2.3']
286    >>> yint.dtype.char
287    'd'
288    >>> x = N.arange(6)
289    >>> y = N.arange(6)
290    >>> xint = N.array([3, 2])
291    >>> yint = interp(y, x, xint, missing=-9999999)
292    >>> ['%.7g' % yint[i] for i in range(len(yint))]
293    ['3', '2']
294    >>> yint.dtype.char
295    'd'
296    """}
297
298
299#- Execute doctest if module is run from command line:
300
301if __name__ == "__main__":
302    """Test the module.
303
304    Tests the examples in all the module documentation strings, plus
305    __test__.
306
307    Note:  To help ensure that module testing of this file works, the
308    parent directory to the current directory is added to sys.path.
309    """
310    import doctest, sys, os
311    sys.path.append(os.pardir)
312    doctest.testmod(sys.modules[__name__])
313
314
315
316
317# ===== end file =====
318
Note: See TracBrowser for help on using the repository browser.