Enhance meshgrid to generate 3D grids + add option for sparse grids (migrated from Trac #967)

[thouis]

Original ticket http://projects.scipy.org/numpy/ticket/967 Reported 2008-12-09 by atmention:pbrod, assigned to unknown.

I would suggest to replace numpy's meshgrid with the meshgrid found in SciTools at:

http://code.google.com/p/scitools/

for the following reasons:

The meshgrid function in NumPy is limited to two dimensions only, while the SciTools version can also work with 3D grids. In addition, the NumPy version of meshgrid has no option for generating sparse grids to conserve memory, like it is in SciTools by specifying the sparse argument:

xv, yv = meshgrid(linspace(-2,2,5), linspace(-1,1,3), sparse=True) xv array([[-2., -1., 0., 1., 2.]]) yv array([[-1.], [ 0.], [ 1.]])

Actually, this is the default behavior for the meshgrid function in SciTools. In NumPy, however, we will in this case get a full 2D grid:

xv, yv = numpy.meshgrid(linspace(-2,2,5), linspace(-1,1,3)) xv array([[-2., -1., 0., 1., 2.], [-2., -1., 0., 1., 2.], [-2., -1., 0., 1., 2.]]) yv array([[-1., -1., -1., -1., -1.], [ 0., 0., 0., 0., 0.], [ 1., 1., 1., 1., 1.]])

This is the same result we get by setting sparse=False in meshgrid in SciTools.

The NumPy functions mgrid and ogrid does provide support for, respectively, full and sparse n-dimensional meshgrids, however, these functions uses slices to generate the meshgrids rather than one-dimensional coordinate arrays such as in Matlab. With slices, the user does not have the option to generate meshgrid with, e.g., irregular spacings, like

x = array([-1,-0.5,1,4,5], float) y = array([0,-2,-5], float) xv, yv = meshgrid(x, y, sparse=False) xv array([[-1. , -0.5, 1. , 4. , 5. ], [-1. , -0.5, 1. , 4. , 5. ], [-1. , -0.5, 1. , 4. , 5. ]]) yv array([[ 0., 0., 0., 0., 0.], [-2., -2., -2., -2., -2.], [-5., -5., -5., -5., -5.]])

In addition to the reasons mentioned above, the meshgrid function in NumPy supports only Cartesian indexing, i.e., x and y, not matrix indexing, i.e., rows and columns (mgrid and ogrid supports only matrix indexing). The meshgrid function in SciTools 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. The difference is illustrated by the following code snippet: nx = 10 ny = 15

x = linspace(-2,2,nx)
y = linspace(-2,2,ny)

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]

It is not entirely true that matrix indexing is not supported by the meshgrid function in NumPy because we can just switch the order of the first two input and output arguments: yv, xv = numpy.meshgrid(y, x) is the same as xv, yv = meshgrid(x, y, sparse=False, indexing='ij') However, we think it is clearer to have the logical "x, y" sequence on the left-hand side and instead adjust a keyword argument.

[thouis]

Attachment in Trac by atmention:pbrod, 2008-12-09: meshgrid.py

[thouis]

Comment in Trac by atmention:charris, 2008-12-09