Mailman 3 June 2010 - NumPy-Discussion

issue #1431 accessing multiple fields in a recarray
by sam tygier June 28, 2010

June 28, 2010

Hello A while ago there was a discussion about field order when accessing recarrays. eg: >>> x = np.array([(1.5,2.5,(1.0,2.0)),(3.,4.,(4.,5.)),(1.,3.,(2.,6.))], dtype=[('x','f4'),('y',np.float32),('value','f4',(2,2))]) >>> x[['x','y']] array([(1.5, 2.5), (3.0, 4.0), (1.0, 3.0)], dtype=[('x', '<f4'), ('y', '<f4')]) >>> x[['y','x']] array([(1.5, 2.5), (3.0, 4.0), (1.0, 3.0)], dtype=[('x', '<f4'), ('y', '<f4')]) http://thread.gmane.org/gmane.comp.python.numeric.general/36933 There is a bug at http://projects.scipy.org/numpy/ticket/1431 The fix is very simple, but there was some worry that some existing code may require the broken behavior. eg some code may ask for: x[['y','x']] and expect the result to have 'x' as the first column in the result. I suspect that code like this is rare (it could be easily changed to x[['x','y']] if thats what is really wanted). I think the correct way to proceed is probably to raise a DeprecationWarning for code that it requiring the broken behavior. The bug has a patch to implement this. Then in a future version the patch to fix it could be applied. Thanks Sam

2 1

printing of array and leading white space ?
by Vincent Davis June 27, 2010

June 27, 2010

This is a little strange and I am not sure what is going on. Look at the number of spaced before the first number in the array. >>> x = np.array([629.54440098249688162, 26186.5470310494529258 ]) >>> x array([ 629.54440098249688162, 26186.5470310494529258 ]) 3 spaces before 629.544... >>> x = np.array([629.54440098249688162, 2634186.5470310494529258 ]) >>> x array([ 6.29544400982496882e+02, 2.63418654703104962e+06]) 2 spaces There are many different ways experiment, and I assume that the spaces are there to help the printed array look better. The 3 spaces are not there is the second digit is shorter. The thing is I can't think of an example where the 3 spaces are needed. Is there an example? Vincent

2 1

arr.copy(order='F') doesn't agree with docstring: what is intended behavior?
by Kurt Smith June 27, 2010

June 27, 2010

I'd really like arr.copy(order='F') to work -- is it supposed to as its docstring says, or is it supposed to raise a TypeError as it does now? This is on numpy 1.4 >>> import numpy as np >>> a = np.arange(10).reshape(5,2) >>> a array([[0, 1], [2, 3], [4, 5], [6, 7], [8, 9]]) >>> print a.copy.__doc__ a.copy(order='C') Return a copy of the array. Parameters ---------- order : {'C', 'F', 'A'}, optional By default, the result is stored in C-contiguous (row-major) order in memory. If `order` is `F`, the result has 'Fortran' (column-major) order. If order is 'A' ('Any'), then the result has the same order as the input. Examples -------- >>> x = np.array([[1,2,3],[4,5,6]], order='F') >>> y = x.copy() >>> x.fill(0) >>> x array([[0, 0, 0], [0, 0, 0]]) >>> y array([[1, 2, 3], [4, 5, 6]]) >>> y.flags['C_CONTIGUOUS'] True >>> a.copy(order='C') Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: copy() takes no keyword arguments >>> a.copy(order='F') Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: copy() takes no keyword arguments >>>

4 5

numpy.random.poisson docs missing "Returns"
by Vincent Davis June 27, 2010

June 27, 2010

numpy.random.poisson docs missing "Returns" http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.poisson.ht…

4 11

Strange behavior of np.sinc
by David Goldsmith June 27, 2010

June 27, 2010

Hi! The docstring for numpy.lib.function_base.sinc indicates that the parameter has to be an ndarray, and that it will return the limiting value 1 for sinc(0). Checking to see if it should actually say array_like, I found the following (Python 2.6): >>> np.sinc(np.array((0,0.5))) array([ 1. , 0.63661977]) >>> np.sinc((0,0.5)) array([ NaN, 0.63661977]) >>> np.sinc([0,0.5]) array([ NaN, 0.63661977]) >>> np.version.version '1.4.1' So, it doesn't choke on non-array sequences, and appears to return values consistent w/ array input, except at 0. Bug in code (failure at 0 if in a sequence) and in the doc (ndarray should be array_like)? DG -- Mathematician: noun, someone who disavows certainty when their uncertainty set is non-empty, even if that set has measure zero. Hope: noun, that delusive spirit which escaped Pandora's jar and, with her lies, prevents mankind from committing a general suicide. (As interpreted by Robert Graves)

2 3

Solving a NLLSQ Problem by Pieces?
by Wayne Watson June 27, 2010

June 27, 2010

The context here is astronomy and optics. The real point is in the last paragraph. I'm looking at a paper that deals with 5 NL (nonlinear) equations and 8 unknown parameters. A. a=a0+arctan((y-y0)/(x-x0) B. z=V*r+S*e**(D*r) r=sqrt((x-x0)**2+(y-y0)**2) and C. cos(z)=cos(u)*cos(z)-sin(u)*sin(ep)*cos(b) sin(a-E) = sin(b)*sin(u)/sin(z) He's trying to estimate parameters of a fisheye lens which has taken star images on a photo plate. For example, x0,y0 is the offset of the center of projection from the zenith (camera not pointing straight up in the sky.) Eq. 2 expresses some nonlinearity in the lens. a0, xo, y0, V, S, D, ep, and E are the parameters. It looks like he uses gradient descent (NLLSQ is nonlinear least squares in Subject.), and takes each equation in turn using the parameter values from the preceding one in the next, B. He provides reasonable initial estimates. A final step uses all eight parameters. He re-examines ep and E, and assigns new estimates. For all (star positions) on the photo plate, he minimizes SUM (Fi**2*Gi) using values from the step for A and B, except for x0,y0. He then does some more dithering, which I'll skip. What I've presented is probably a bit difficult to understand without a good optics understanding, but my question is something like is this commonly done to solve a system of NLLSQ? It looks a bit wild. I guess if one knows his subject well, then bringing some "extra" knowledge to the process helps. As I understand it, he solves parameters in A, then uses them in B, and so on. I guess that's a reasonable way to do it. -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7" N, 121° 2' 32" W, 2700 feet Air travel safety Plus Three/Minus 8 rule. Eighty % of crashes take place in the first 3 minutes and last 8 minutes. Survival chances are good in you're paying attention. No hard drinks prior to those periods. No sleeping pills either. Sit within 5 rows of an exit door. --The Survivors Club, Ben Sherwood Web Page:<www.speckledwithstars.net/>

2 3

Try this, new tests for numpy.random
by Vincent Davis June 26, 2010

June 26, 2010

I added a set of simple testes for numpy.random. These test only test that you get results and that they are stable. (see attachment) That said I get several errors (see below) as a result I assume of how the resent enthought 6.2 beta does the calculation (The tests results are from numpy running of 3.1) I think it would be a good idea if several others ran the test to see if there are other problems or if it is actually stable across versions. Any other feedback Vincent SUMMARY OF UNIT TEST RESULTS ============================ FAILED TESTS ------------ * TestRandomDist.test_lognormal -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 1.81865312, 0.07476378], [ 0.51747188, 0.61901346], [ 44.09839408, 0.35777415]]) y: array([[ 16.50698632, 36.54846706], [ 22.678866 , 0.71617561], [ 65.72798502, 86.84341601]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 272 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_multivariate_normal -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 50.0%) x: array([[[ -0.15876548, 10. ], [ 0.58858952, 10. ]], ... y: array([[[ -1.47027513, 10. ], [ -1.65915082, 10. ]], ... unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 305 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_negative_binomial -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not equal (mismatch 100.0%) x: array([[702, 777], [692, 646], [630, 806]]) y: array([[848, 841], [892, 611], [779, 647]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 313 np.testing.assert_array_equal(actual, desired) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 677 verbose=verbose, header='Arrays are not equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_noncentral_chisquare -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 5.42550628, 3.68327961], [ 7.8305713 , 5.45961594], [ 10.68988949, 5.17932387]]) y: array([[ 23.91905354, 13.35324693], [ 31.22452661, 16.60047399], [ 5.03461598, 17.94973089]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 321 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_noncentral_f -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 1.41010542, 106.52208522], [ 1.35761411, 45.00197918], [ 9.97943281, 13.86047453]]) y: array([[ 1.405981 , 0.34207973], [ 3.57715069, 7.92632663], [ 0.43741599, 1.17742088]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 330 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_normal -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 0.59809618, -2.59342171], [-0.6588001 , -0.47962827], [ 3.78642337, -1.02785335]]) y: array([[ 2.8037837 , 3.59863924], [ 3.12143348, -0.33382988], [ 4.18552479, 4.46410668]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 338 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_rand -- test_random.py (numpy-vmd-dev/numpy/random/tests) NameError: global name 'array' is not defined unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 113 desired1 = array([[ 0.61879477158567997, 0.59162362775974664], * TestRandomDist.test_randint -- test_random.py (numpy-vmd-dev/numpy/random/tests) AttributeError: 'module' object has no attribute 'randint' unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 129 actual = np.randint(-99, 99, size=(3,2)) * TestRandomDist.test_randn -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 0.2373197 , -1.35843925], [-0.39112844, -0.30154253], [ 1.83148329, -0.57565507]]) y: array([[ 1.34016346, 1.73759123], [ 1.49898834, -0.22864333], [ 2.031034 , 2.17032495]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 125 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_standard_cauchy -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[-0.17470026, 1.29709215], [-3.18156373, -0.45005445], [-0.46997791, -1.57823221]]) y: array([[ 0.7712766 , -6.55601162], [ 0.93582023, -2.07479293], [-4.74601644, 0.18338989]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 378 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_standard_gamma -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 3.0732848 , 1.30782909], [ 3.57713935, 0.9544672 ], [ 1.09782087, 3.73036312]]) y: array([[ 5.50841531, 6.6295347 ], [ 5.93988485, 2.31044849], [ 7.54838614, 8.01275609]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 394 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_standard_normal -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 0.2373197 , -1.35843925], [-0.39112844, -0.30154253], [ 1.83148329, -0.57565507]]) y: array([[ 1.34016346, 1.73759123], [ 1.49898834, -0.22864333], [ 2.031034 , 2.17032495]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 402 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_standard_t -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 0.34958651, -0.07515709], [-0.36478493, -0.16467516], [ 0.65905206, -1.33550263]]) y: array([[ 0.97140612, -0.08830487], [ 1.36311144, -0.55317464], [-0.18473749, 0.61181537]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 410 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) * TestRandomDist.test_wald -- test_random.py (numpy-vmd-dev/numpy/random/tests) AssertionError: Arrays are not almost equal (mismatch 100.0%) x: array([[ 0.99532782, 0.5024266 ], [ 0.78691891, 0.37161609], [ 0.42210803, 0.73660839]]) y: array([[ 3.82935266, 5.13125249], [ 0.35045404, 1.50832397], [ 0.2412432 , 0.22031101]]) unittest.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6): 279 testMethod() test_random.py (numpy-vmd-dev/numpy/random/tests): 444 np.testing.assert_array_almost_equal(actual, desired, decimal=15) utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 765 header='Arrays are not almost equal') utils.py (/Library/Frameworks/Python.framework/Versions/6.2/lib/python2.6/site-packages/numpy/testing): 609 raise AssertionError(msg) PASSED TESTS ------------ * TestMultinomial.test_basic -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestMultinomial.test_int_negative_interval -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestMultinomial.test_zero_probability -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_beta -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_binomial -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_bytes -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_chisquare -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_dirichlet -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_exponential -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_f -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_gamma -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_geometric -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_gumbel -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_hypergeometric -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_laplace -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_logistic -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_logseries -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_multinomial -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_pareto -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_poisson -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_power -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_random_integers -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_random_sample -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_rayleigh -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_shuffle -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_standard_exponential -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_triangular -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_uniform -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_vonmises -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_weibull -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRandomDist.test_zipf -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRegression.test_VonMises_range -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRegression.test_hypergeometric_range -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestRegression.test_logseries_convergence -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestSetState.test_backwards_compatibility -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestSetState.test_basic -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestSetState.test_gaussian_reset -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestSetState.test_gaussian_reset_in_media_res -- test_random.py (numpy-vmd-dev/numpy/random/tests) * TestSetState.test_negative_binomial -- test_random.py (numpy-vmd-dev/numpy/random/tests)

1 0

numpy.random.pareto, missing "returns" in the docs
by Vincent Davis June 26, 2010

June 26, 2010

numpy.random.pareto, missing "returns" in the docs http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.pareto.htm… Vincent

1 0

Documentation error in numpy.random.logseries
by Vincent Davis June 26, 2010

June 26, 2010

numpy.random.logseries(p, size=None) but the parameters section, Parameters: loc : float scale : float > 0. size : {tuple, int} Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. Notice that p <> loc and what about scale. I'll file a ticket unless I am mission something, Which should it be loc or p What about scale. There is no numpy-dev list right? Should this list be used or the scipy-dev list> Vincent

3 7

2d binning and linear regression
by Tom Durrant June 23, 2010

June 23, 2010

Hi All, I have a problem involving lat/lon data. Basically, I am evaluating numerical weather model data against satellite data, and trying to produce gridded plots of various statistics. There are various steps involved with this, but basically, I get to the point where I have four arrays of the same length, containing lat, lon, model and observed values respectively along the path of the sattelite. eg: lat = [ 50.32 50.78 51.24 51.82 52.55 53.15 53.75 54.28 54.79 55.16 ... ] lon = [ 192.83 193.38 193.94 194.67 195.65 196.49 197.35 198.15 198.94 199.53 ... ] obs = [ 1.42 1.35 1.55 1.50 1.59 1.76 2.15 1.90 1.55 0.73 ... ] model = [ 1.67 1.68 1.70 1.79 2.04 2.36 2.53 2.38 2.149 1.57 ... ] I then want to calculate statistics based on bins of say 2 X 2 degree boxes. These arrays are very large, on the order of 10^6. For ease of explanation, I will focus initially on bias. My first approach was to use loops through lat/lon boxes and use np.where statements to extract all the values of the model and observations within the given box, and calculate the bias as the mean of the difference. This was obviously very slow. histogram2d provided a FAR better way to do this. i.e. import numpy as np latedges=np.arange(-90,90,2) lonedges=np.arange(0,360,2) diff = model-obs grid_N, xedges, yedges = np.histogram2d(lat, lon], bins=(latedges,lonedges)) grid_bias_sum, xedges, yedges = np.histogram2d(lat, lon, weights=diff, bins=(latedges,lonedges)) grid_bias = grid_bias_sum/grid_N I now want to determine the the linear regression coefficients for each each box after fitting a least squares linear regression to the data in each bin. I have been looking at using np.digitize to extract the bin indeces, but haven't had much success trying to do this in two dimensions. I am back to looping through the lat and lon box values, using np.where to extract the observations and model data within that box, and using np.polyfit to obtain the regression coefficients. This is, of course, very slow. Can anyone advise a smarter, vectorized way to do this? Any thoughts would be greatly appreciated. Thanks in advance Tom

6 13