Skip to content

Loading…

Meshgrid enhancements (>2-D, sparse grids, matrix indexing) #192

Merged
merged 5 commits into from

5 participants

@charris charris commented on an outdated diff
numpy/lib/function_base.py
((15 lines not shown))
Parameters
----------
- x, y : ndarray
- Two 1-D arrays representing the x and y coordinates of a grid.
+ x1, x2,..., xn : array_like
+ 1-D arrays representing the coordinates of a grid.
+ indexing : {'xy', 'ij'}, optional
+ Cartesian ('xy', default) or matrix ('ij') indexing of output.
+ See Notes for more details.
+ sparse : bool, optional
+ If True a sparse grid is returned in order to conserve memory.
+ Default is False.
+ copy : bool, optional
+ If False, a view into the original arrays are returned in
+ order to conserve memory. Default is True. Please note that
+ ``sparse=False, copy=False`` will likely return non-contiguous arrays.
+ Furthermore, more than one element of a broadcasted array may refer to
@charris NumPy member
charris added a note

broadcast, not broadcasted. Although languages do evolve ;)

@rgommers NumPy member
rgommers added a note

fixed

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@charris charris commented on the diff
numpy/lib/function_base.py
((40 lines not shown))
- return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays
- with the elements of `x` and y repeated to fill the matrix along
- the first dimension for `x`, the second for `y`.
+ X1, X2,..., XN : ndarray
+ For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
+ return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
+ or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
+ with the elements of `xi` repeated to fill the matrix along
+ the first dimension for `x1`, the second for `x2` and so on.
+
+ Notes
+ -----
+ This function supports both indexing conventions through the indexing keyword
+ argument. Giving the string 'ij' returns a meshgrid with matrix indexing,
+ while 'xy' returns a meshgrid with Cartesian indexing. In the 2-D case
+ with inputs of length M and N, the outputs are of shape (N, M) for 'xy'
@charris NumPy member
charris added a note

Maybe explain this with a connection to image plots, where the 'x' typically runs left to right and is the 'fast' index. Or is that what this is about. Where would one use 'ij' ?

@teoliphant NumPy member

I agree. I find the "matrix indexing" versus "Cartesian indexing" description hard to understand. I think you are talking about "image" conventions versus "plotting" conventions. Is that correct?

@rgommers NumPy member
rgommers added a note

Correct, image convention vs. just standard array/matrix indexing. Since this is the most confusing part apparently (see also discussion on this point at http://projects.scipy.org/numpy/ticket/966), I'm tempted to just leave out the indexing keyword after all. This seems to be confusing for Matlab users too, see for example the comments at http://blogs.mathworks.com/loren/2007/06/21/indexing-terminology/.

There really isn't an ideal solution here.

Note that mgrid and ogrid use 'ij', while currently meshgrid uses 'xy':

In [35]: mgrid[0:2, 3:5]
Out[35]: 
array([[[0, 0],
        [1, 1]],

       [[3, 4],
        [3, 4]]])

In [36]: meshgrid([0, 1], [3, 4])
Out[36]: 
[array([[0, 1],
       [0, 1]]), array([[3, 3],
       [4, 4]])]
@charris NumPy member
charris added a note

Maybe call ir 'image' indexing? I find in practice it is easier to settle on a single convention ('ij'), always put x in the first index, and just transpose for showing images. Mixing the conventions leads to all sorts of confusion.

@rgommers NumPy member
rgommers added a note

Only using 'ij' would be best probably, but since currently meshgrid uses 'xy' we can't do that.

@charris NumPy member
charris added a note

Oops, see, I already got confused ;) I think a different choice of names would clarify things, 'xy' seems ok, 'standard' would be another option, and then maybe 'image' for the second choice, with a short explanation (images are displayed left-right, top-bottom) and a translation for Matlab users.

@rgommers NumPy member

ij --> image sounds OK to me. Although I forgot the details by now. With added explanation it should work.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@charris charris commented on the diff
numpy/lib/function_base.py
((134 lines not shown))
+ sparse = kwargs.get('sparse', False)
+ indexing = kwargs.get('indexing', 'xy')
+ if not indexing in ['xy', 'ij']:
+ raise ValueError("Valid values for `indexing` are 'xy' and 'ij'.")
+
+ s0 = (1,) * ndim
+ output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)]
+
+ shape = [x.size for x in output]
+
+ if indexing == 'xy':
+ # switch first and second axis
+ output[0].shape = (1, -1) + (1,)*(ndim - 2)
+ output[1].shape = (-1, 1) + (1,)*(ndim - 2)
+ shape[0], shape[1] = shape[1], shape[0]
+
@charris NumPy member
charris added a note

So 'xy' indexing only affects the first two axis?

@rgommers NumPy member
rgommers added a note

yes

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@charris charris commented on an outdated diff
numpy/lib/tests/test_function_base.py
((19 lines not shown))
+ z = [8, 9]
+ assert_(meshgrid(x, y)[0].shape == (4, 3))
+ assert_(meshgrid(x, y, indexing='ij')[0].shape == (3, 4))
+ assert_(meshgrid(x, y, z)[0].shape == (4, 3, 2))
+ assert_(meshgrid(x, y, z, indexing='ij')[0].shape == (3, 4, 2))
+
+ assert_raises(ValueError, meshgrid, x, y, indexing='notvalid')
+
+ def test_sparse(self):
+ [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], sparse=True)
+ assert_(all(X == array([[1, 2, 3]])))
+ assert_(all(Y == array([[4],
+ [5],
+ [6],
+ [7]])))
+
@charris NumPy member
charris added a note

Any reason not to put that on a single line?

@rgommers NumPy member
rgommers added a note

not really, changed that now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@teoliphant
NumPy member

It seems that this PR can be merged. Confusion in the naming can be improved later (and already exists between the default mgrid and meshgrid conventions).

@rgommers
NumPy member

Would be better to fix it before/when merging imho, otherwise it will simply be forgotten. Sorry I haven't made the time to finish this off yet.

@lukecampbell

Looks really good, can't wait for this to be merged.

@teoliphant teoliphant merged commit 0b2bfa9 into numpy:master
@teoliphant teoliphant added a commit to teoliphant/numpy that referenced this pull request
@teoliphant teoliphant BUG: Fix some tests in PR #192 5b4e61b
@rgommers rgommers deleted the rgommers:meshgrid3d branch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Showing with 134 additions and 33 deletions.
  1. +106 −33 numpy/lib/function_base.py
  2. +28 −0 numpy/lib/tests/test_function_base.py
View
139 numpy/lib/function_base.py
@@ -3200,62 +3200,135 @@ def add_newdoc(place, obj, doc):
pass
-# From matplotlib
-def meshgrid(x,y):
+# Based on scitools meshgrid
+def meshgrid(*xi, **kwargs):
"""
- Return coordinate matrices from two coordinate vectors.
+ Return coordinate matrices from two or more coordinate vectors.
+
+ Make N-D coordinate arrays for vectorized evaluations of
+ N-D scalar/vector fields over N-D grids, given
+ one-dimensional coordinate arrays x1, x2,..., xn.
Parameters
----------
- x, y : ndarray
- Two 1-D arrays representing the x and y coordinates of a grid.
+ x1, x2,..., xn : array_like
+ 1-D arrays representing the coordinates of a grid.
+ indexing : {'xy', 'ij'}, optional
+ Cartesian ('xy', default) or matrix ('ij') indexing of output.
+ See Notes for more details.
+ sparse : bool, optional
+ If True a sparse grid is returned in order to conserve memory.
+ Default is False.
+ copy : bool, optional
+ If False, a view into the original arrays are returned in
+ order to conserve memory. Default is True. Please note that
+ ``sparse=False, copy=False`` will likely return non-contiguous arrays.
+ Furthermore, more than one element of a broadcast array may refer to
+ a single memory location. If you need to write to the arrays, make
+ copies first.
Returns
-------
- X, Y : ndarray
- For vectors `x`, `y` with lengths ``Nx=len(x)`` and ``Ny=len(y)``,
- return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays
- with the elements of `x` and y repeated to fill the matrix along
- the first dimension for `x`, the second for `y`.
+ X1, X2,..., XN : ndarray
+ For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
+ return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
+ or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
+ with the elements of `xi` repeated to fill the matrix along
+ the first dimension for `x1`, the second for `x2` and so on.
+
+ Notes
+ -----
+ This function supports both indexing conventions through the indexing keyword
+ argument. Giving the string 'ij' returns a meshgrid with matrix indexing,
+ while 'xy' returns a meshgrid with Cartesian indexing. In the 2-D case
+ with inputs of length M and N, the outputs are of shape (N, M) for 'xy'
@charris NumPy member
charris added a note

Maybe explain this with a connection to image plots, where the 'x' typically runs left to right and is the 'fast' index. Or is that what this is about. Where would one use 'ij' ?

@teoliphant NumPy member

I agree. I find the "matrix indexing" versus "Cartesian indexing" description hard to understand. I think you are talking about "image" conventions versus "plotting" conventions. Is that correct?

@rgommers NumPy member
rgommers added a note

Correct, image convention vs. just standard array/matrix indexing. Since this is the most confusing part apparently (see also discussion on this point at http://projects.scipy.org/numpy/ticket/966), I'm tempted to just leave out the indexing keyword after all. This seems to be confusing for Matlab users too, see for example the comments at http://blogs.mathworks.com/loren/2007/06/21/indexing-terminology/.

There really isn't an ideal solution here.

Note that mgrid and ogrid use 'ij', while currently meshgrid uses 'xy':

In [35]: mgrid[0:2, 3:5]
Out[35]: 
array([[[0, 0],
        [1, 1]],

       [[3, 4],
        [3, 4]]])

In [36]: meshgrid([0, 1], [3, 4])
Out[36]: 
[array([[0, 1],
       [0, 1]]), array([[3, 3],
       [4, 4]])]
@charris NumPy member
charris added a note

Maybe call ir 'image' indexing? I find in practice it is easier to settle on a single convention ('ij'), always put x in the first index, and just transpose for showing images. Mixing the conventions leads to all sorts of confusion.

@rgommers NumPy member
rgommers added a note

Only using 'ij' would be best probably, but since currently meshgrid uses 'xy' we can't do that.

@charris NumPy member
charris added a note

Oops, see, I already got confused ;) I think a different choice of names would clarify things, 'xy' seems ok, 'standard' would be another option, and then maybe 'image' for the second choice, with a short explanation (images are displayed left-right, top-bottom) and a translation for Matlab users.

@rgommers NumPy member

ij --> image sounds OK to me. Although I forgot the details by now. With added explanation it should work.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ indexing and (M, N) for 'ij' indexing. In the 3-D case with inputs of
+ length M, N and P, outputs are of shape (N, M, P) for 'xy' indexing and (M,
+ N, P) for 'ij' indexing. The difference is illustrated by the following
+ code snippet::
+
+ xv, yv = meshgrid(x, y, sparse=False, indexing='ij')
+ for i in range(nx):
+ for j in range(ny):
+ # treat xv[i,j], yv[i,j]
+
+ xv, yv = meshgrid(x, y, sparse=False, indexing='xy')
+ for i in range(nx):
+ for j in range(ny):
+ # treat xv[j,i], yv[j,i]
See Also
--------
index_tricks.mgrid : Construct a multi-dimensional "meshgrid"
- using indexing notation.
+ using indexing notation.
index_tricks.ogrid : Construct an open multi-dimensional "meshgrid"
- using indexing notation.
+ using indexing notation.
Examples
--------
- >>> X, Y = np.meshgrid([1,2,3], [4,5,6,7])
- >>> X
- array([[1, 2, 3],
- [1, 2, 3],
- [1, 2, 3],
- [1, 2, 3]])
- >>> Y
- array([[4, 4, 4],
- [5, 5, 5],
- [6, 6, 6],
- [7, 7, 7]])
+ >>> nx, ny = (3, 2)
+ >>> x = np.linspace(0, 1, nx)
+ >>> y = np.linspace(0, 1, ny)
+ >>> xv, yv = meshgrid(x, y)
+ >>> xv
+ array([[ 0. , 0.5, 1. ],
+ [ 0. , 0.5, 1. ]])
+ >>> yv
+ array([[ 0., 0., 0.],
+ [ 1., 1., 1.]])
+ >>> xv, yv = meshgrid(x, y, sparse=True) # make sparse output arrays
+ >>> xv
+ array([[ 0. , 0.5, 1. ]])
+ >>> yv
+ array([[ 0.],
+ [ 1.]])
`meshgrid` is very useful to evaluate functions on a grid.
>>> x = np.arange(-5, 5, 0.1)
>>> y = np.arange(-5, 5, 0.1)
- >>> xx, yy = np.meshgrid(x, y)
- >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2)
+ >>> xx, yy = meshgrid(x, y, sparse=True)
+ >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
+ >>> h = plt.contourf(x,y,z)
"""
- x = asarray(x)
- y = asarray(y)
- numRows, numCols = len(y), len(x) # yes, reversed
- x = x.reshape(1,numCols)
- X = x.repeat(numRows, axis=0)
+ if len(xi) < 2:
+ msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(len(xi) > 0)
+ raise ValueError(msg)
+
+ args = np.atleast_1d(*xi)
+ ndim = len(args)
+
+ copy_ = kwargs.get('copy', True)
+ sparse = kwargs.get('sparse', False)
+ indexing = kwargs.get('indexing', 'xy')
+ if not indexing in ['xy', 'ij']:
+ raise ValueError("Valid values for `indexing` are 'xy' and 'ij'.")
+
+ s0 = (1,) * ndim
+ output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)]
+
+ shape = [x.size for x in output]
+
+ if indexing == 'xy':
+ # switch first and second axis
+ output[0].shape = (1, -1) + (1,)*(ndim - 2)
+ output[1].shape = (-1, 1) + (1,)*(ndim - 2)
+ shape[0], shape[1] = shape[1], shape[0]
+
@charris NumPy member
charris added a note

So 'xy' indexing only affects the first two axis?

@rgommers NumPy member
rgommers added a note

yes

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ if sparse:
+ if copy_:
+ return [x.copy() for x in output]
+ else:
+ return output
+ else:
+ # Return the full N-D matrix (not only the 1-D vector)
+ if copy_:
+ mult_fact = np.ones(shape, dtype=int)
+ return [x * mult_fact for x in output]
+ else:
+ return np.broadcast_arrays(*output)
- y = y.reshape(numRows,1)
- Y = y.repeat(numCols, axis=1)
- return X, Y
def delete(arr, obj, axis=None):
"""
View
28 numpy/lib/tests/test_function_base.py
@@ -1026,6 +1026,34 @@ def test_simple(self):
[6, 6, 6],
[7, 7, 7]])))
+ def test_single_input(self):
+ assert_raises(ValueError, meshgrid, np.arange(5))
+
+ def test_indexing(self):
+ x = [1, 2, 3]
+ y = [4, 5, 6, 7]
+ [X, Y] = meshgrid(x, y, indexing='ij')
+ assert_(all(X == array([[1, 1, 1, 1],
+ [2, 2, 2, 2],
+ [3, 3, 3, 3]])))
+ assert_(all(Y == array([[4, 5, 6, 7],
+ [4, 5, 6, 7],
+ [4, 5, 6, 7]])))
+
+ # Test expected shapes:
+ z = [8, 9]
+ assert_(meshgrid(x, y)[0].shape == (4, 3))
+ assert_(meshgrid(x, y, indexing='ij')[0].shape == (3, 4))
+ assert_(meshgrid(x, y, z)[0].shape == (4, 3, 2))
+ assert_(meshgrid(x, y, z, indexing='ij')[0].shape == (3, 4, 2))
+
+ assert_raises(ValueError, meshgrid, x, y, indexing='notvalid')
+
+ def test_sparse(self):
+ [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], sparse=True)
+ assert_(all(X == array([[1, 2, 3]])))
+ assert_(all(Y == array([[4], [5], [6], [7]])))
+
class TestPiecewise(TestCase):
def test_simple(self):
Something went wrong with that request. Please try again.