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 | |
---|
65 | def 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 | |
---|
301 | if __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 | |
---|