[5390] | 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: |
---|
[6304] | 77 | * y: Ordinate values of data. Rank 1 numeric vector. Required. |
---|
[5390] | 78 | Can have missing values. Floating or integer type. |
---|
| 79 | |
---|
[6304] | 80 | * x: Abscissa values of data. Rank 1 numeric vector. Required. |
---|
[5390] | 81 | Can have no missing values. Must be monotonically ascending. |
---|
| 82 | Floating or integer type. |
---|
| 83 | |
---|
| 84 | * xinterp: Abscissa values to calculate interpolated ordinate |
---|
[6304] | 85 | values at. Rank 1 numeric vector or numeric scalar. Required. |
---|
[5390] | 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: |
---|
[6304] | 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 |
---|
[5390] | 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 |
---|
[6304] | 113 | >>> import numpy as N |
---|
[5390] | 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 |
---|
[6304] | 140 | import numpy.ma MA |
---|
| 141 | import numpy as N |
---|
[5390] | 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 |
---|
[6304] | 218 | >>> import numpy as N |
---|
[5390] | 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'] |
---|
[6304] | 286 | >>> yint.dtype.char |
---|
[5390] | 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'] |
---|
[6304] | 294 | >>> yint.dtype.char |
---|
[5390] | 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 | |
---|